anchor.diff:     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)
  -       ! y = Cx + Du + Fy  - TODO consider using lapack for efficiency
  -       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))
  +
  +       ! >>> 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
  
positive.code:                                                     !    states (z)
     PUBLIC :: ExtPtfm_GetOP                          !  Routine to get the operating-point values for linearization (from data structures to arrays)
  
  CONTAINS
     
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
negative.code:       enddo
        ! LinNames_x
        do I=1,p%nCB; 
            InitOut%LinNames_x(I)       = 'Mode '//trim(Num2LStr(p%ActiveCBDOF(I)))//' displacement, -';
            InitOut%LinNames_x(I+p%nCB) = 'Mode '//trim(Num2LStr(p%ActiveCBDOF(I)))//' velocity, -';
        enddo
        ! 
  