anchor.diff:                                                      !    states (z)
      PUBLIC :: ExtPtfm_GetOP                          !  Routine to get the operating-point values for linearization (from data structures to arrays)
   
  +
  +
  +   INTERFACE LAPACK_COPY    
  +       SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
  +           USE Precision, only: R8Ki
  +           INTEGER    :: INCX,INCY,N
  +           real(R8Ki) :: DX(*),DY(*)
  +       ENDSUBROUTINE
  +       SUBROUTINE SCOPY(N,X,INCX,Y,INCY)
  +           USE Precision, only: SiKi
  +           INTEGER    :: INCX,INCY,N
  +           real(SiKi) :: X(*),Y(*)
  +       ENDSUBROUTINE
  +   END INTERFACE
  +   INTERFACE LAPACK_GEMV    
  +       SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  +           USE Precision, only: R8Ki
  +           real(R8Ki) :: ALPHA,BETA
  +           integer    :: INCX,INCY,LDA,M,N
  +           character  :: TRANS
  +           real(R8Ki) :: A(LDA,*),X(*),Y(*)
  +       ENDSUBROUTINE
  +       SUBROUTINE SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  +           USE Precision, only: SiKi
  +           real(SiKi) :: ALPHA,BETA
  +           integer    :: INCX,INCY,LDA,M,N
  +           character  :: TRANS
  +           real(SiKi) :: A(LDA,*),X(*),Y(*)
  +       ENDSUBROUTINE
  +   END INTERFACE LAPACK_GEMV    
  +   INTERFACE LAPACK_AXPY    
  +      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
  +           USE Precision, only: R8Ki
  +           real(R8Ki) :: DA
  +           integer    :: INCX,INCY,N
  +           real(R8Ki) :: DX(*),DY(*)
  +       ENDSUBROUTINE
  +      SUBROUTINE SAXPY(N,A,X,INCX,Y,INCY)
  +           USE Precision, only: SiKi
  +           real(SiKi) :: A
  +           integer    :: INCX,INCY,N
  +           real(SiKi) :: X(*),Y(*)
  +       ENDSUBROUTINE
  +   END INTERFACE
   CONTAINS
      
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
positive.code:       enddo
        ! LinNames_x
        do I=1,p%nCB; 
            InitOut%LinNames_x(I)       = 'Mode '//trim(Num2LStr(p%ActiveCBDOF(I)))//' displacement, -'; ! TODO  TODO IDOF
            InitOut%LinNames_x(I+p%nCB) = 'Mode '//trim(Num2LStr(p%ActiveCBDOF(I)))//' velocity, -';
        enddo
        ! 
  
negative.code:    m%uFlat(10:12) = u%PtfmMesh%RotationVel   (:,1)
     m%uFlat(13:15) = u%PtfmMesh%TranslationAcc(:,1)
     m%uFlat(16:18) = u%PtfmMesh%RotationAcc   (:,1)
  
     !--- Computing output:  y = Cx + Du + Fy  
     ! 
     if (p%nCB>0) then
         ! x flat
         m%xFlat(      1:p%nCB  ) = x%qm   (1:p%nCB)
         m%xFlat(p%nCB+1:2*p%nCB) = x%qmdot(1:p%nCB)
  
         ! >>> MATMUL implementation
         !Fc = matmul(p%CMat, m%xFlat) + matmul(p%DMat, m%uFlat) + m%F_at_t(1:6) - matmul(p%M12, m%F_at_t(6+1:6+p%nCB))
  
         ! >>> LAPACK implementation
         Fc(1:6) = m%F_at_t(1:6) ! Fc = F1r + ...
         !           GEMV(TRS, M  , N      , alpha    , A     , LDA, X                    ,INCX, Beta  ,  Y, IncY)
         CALL LAPACK_GEMV('n', 6  , 2*p%nCB,  1.0_ReKi, p%CMat, 6  , m%xFlat              , 1, 1.0_ReKi, Fc, 1   ) ! = C*x + (F1r)
         CALL LAPACK_GEMV('n', 6  ,   18   ,  1.0_ReKi, p%DMat, 6  , m%uFlat              , 1, 1.0_ReKi, Fc, 1   ) ! + D*u
         CALL LAPACK_GEMV('n', 6  , p%nCB  , -1.0_ReKi, p%M12 , 6  , m%F_at_t(6+1:6+p%nCB), 1, 1.0_ReKi, Fc, 1   ) ! - M12*F2r
     else
         Fc =                           matmul(p%DMat, m%uFlat) + m%F_at_t(1:6) 
     endif
  