anchor.diff:                                CBB=CBB, CBM=CBM, CMM=CMM)  !< Optional Outputs
      if(Failed()) return
   
  -
      ! --- Reduction of force if provided
      if (present(FG).and.present(FGR).and.present(FGL)) then
  +      if (present(FGB)) then
  +         !FGB = FGR + matmul( transpose(PhiR), FGL)
  +         if (nL>0) then
  +            CALL LAPACK_GEMV('t', nL  , nR,  1.0_FeKi, PhiR, nL, FGL, 1, 0.0_FeKi, FGB, 1 )
  +            FGB = FGR + FGB
  +         else
  +            FGB = FGR 
  +         endif
  +      endif
         if (present(FGM)) then
  -         FGB = FGR + matmul( transpose(PhiR), FGL)
  -      endif
  -      if (present(FGB)) then
  -         FGM = matmul( FGL, PhiL(:,1:nM) ) != matmul( transpose(PhiM), FGL ) because FGL is 1-D
  +         !FGM = matmul( FGL, PhiL(:,1:nM) ) != matmul( transpose(PhiM), FGL ) because FGL is 1-D
  +         if (nM>0) then
  +            CALL LAPACK_GEMV('t', nL  , nM,  1.0_FeKi, PhiL(:,1:nM), nL, FGL, 1, 0.0_FeKi, FGM, 1 )
  +         endif
         endif
      endif
      call CleanUp()
  
positive.code:    end do
     ! Compute Ainv = 1.0*V^t * U^t + 0.0*Ainv     V*(inv(S))*U' 
     !call DGEMM( 'T', 'T', N, M, K, 1.0, V, K, U, M, 0.0, Ainv, N)
     call LAPACK_GEMM( 'T', 'T', 1.0_FEKi, Vt, U, 0.0_FEKi, Ainv, ErrStat, ErrMsg)
     ! --- Compute rank
     !tol=maxval(shape(A))*epsilon(maxval(S))
  
negative.code: 
        if (present(CRR)) then
           ! Guyan damping CBB = CRR + (CRL*PhiR) + (CRL*PhiR)^T + PhiR^T*CLL*PhiR
           ! PhiR_T_CLL = TRANSPOSE(PhiR) * CLL
           CALL AllocAry( Temp , nR, nR, 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
           CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, PhiR, MLL, 0.0_FeKi, PhiR_T_MLL  , ErrStat2, ErrMsg2); if(Failed()) return
           ! CBB = MATMUL(CRL, PhiR)
           CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, CRL, PhiR, 0.0_FeKi, CBB  , ErrStat2, ErrMsg2); if(Failed()) return
           ! CBB2 = MATMUL( PhiR_T_CLL, PhiR )
           CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, PhiR_T_MLL, PhiR, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
           CBB = CRR + CBB + TRANSPOSE( CBB ) + Temp
           DEALLOCATE(Temp)
           ! Cross coupling CMB = PhiM^T*CLR + PhiM^T CLL PhiR
           !                CBM = CRL*PhiM + PhiR^T CLL^T PhiM (NOTE: assuming CLL symmetric)
           IF ( nM == 0) THEN
  