anchor.diff:        
         ! --- Set MBB, MBM, and KBB from Eq. 4:
         CALL AllocAry( PhiR_T_MLL,  nR, nL, 'PhiR_T_MLL', ErrStat2, ErrMsg2); if(Failed()) return
  +      CALL AllocAry( Temp , nR, nR, 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
            
  -      PhiR_T_MLL = TRANSPOSE(PhiR)
  -      PhiR_T_MLL = MATMUL(PhiR_T_MLL, MLL)
  -      !CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, PhiR_T_MLL, MLL, 0.0_FeKi, PhiR_T_MLL  , ErrStat2, ErrMsg2); if(Failed()) return
  -      MBB = MATMUL(MRL, PhiR)
  +      ! PhiR_T_MLL = TRANSPOSE(PhiR) * MLL
  +      CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, PhiR, MLL, 0.0_FeKi, PhiR_T_MLL  , ErrStat2, ErrMsg2); if(Failed()) return
  +      ! MBB1 = MATMUL(MRL, PhiR)
         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, MRL, PhiR, 0.0_FeKi, MBB  , ErrStat2, ErrMsg2); if(Failed()) return
  -      MBB = MRR + MBB + TRANSPOSE( MBB ) + MATMUL( PhiR_T_MLL, PhiR )
  +      ! MBB2 = MATMUL( PhiR_T_MLL, PhiR )
  +      CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, PhiR_T_MLL, PhiR, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  +      MBB = MRR + MBB + TRANSPOSE( MBB ) + Temp
  +      DEALLOCATE(Temp)
            
         IF ( nM == 0) THEN
            MBM = 0.0_FEKi
         ELSE
  -         MBM = MATMUL( PhiR_T_MLL, PhiL(:,1:nM))  ! last half of operation
  -         MBM = MATMUL( MRL, PhiL(:,1:nM) ) + MBM    !This had PhiM      
  +         CALL AllocAry( Temp , nR, nM, 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
  +         !MBM = MATMUL( PhiR_T_MLL, PhiL(:,1:nM))  ! last half of operation
  +         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, PhiR_T_MLL, PhiL(:,1:nM), 0.0_FeKi, MBM  , ErrStat2, ErrMsg2); if(Failed()) return
  +         ! Temp = MATMUL( MRL, PhiL(:,1:nM) )
  +         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, MRL, PhiL(:,1:nM), 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  +         MBM = Temp + MBM    !This had PhiM      
  +         DEALLOCATE(Temp)
         ENDIF
         
         !KBB = MATMUL(KRL, PhiR)   
  
positive.code:                               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(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
        endif
     endif
     call CleanUp()
  
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
  