anchor.diff:        ! bCheckSingularity = True
         CALL EigenSolveWrap(KLL, MLL, nL, nM_out, .True., PhiL(:,1:nM_out), OmegaL(1:nM_out),  ErrStat2, ErrMsg2); if(Failed()) return
         ! --- Normalize PhiL
  -      ! bjj: break up this equation to avoid as many temporary variables on the stack
         ! MU = MATMUL ( MATMUL( TRANSPOSE(PhiL), MLL ), PhiL )
  -      CALL AllocAry( Temp , nL    , nL , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
  -      CALL AllocAry( MU   , nM_out, nL , 'Mu'   , ErrStat2 , ErrMsg2); if(Failed()) return
  -      MU   = TRANSPOSE(PhiL)
  -      !Temp = MATMUL( MU, MLL )
  -      CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, MU, MLL, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  -      !MU(1:nM_Out,1:nM_Out) = MATMUL( Temp, PhiL )
  +      CALL AllocAry( Temp , nM_out, nL     , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
  +      CALL AllocAry( MU   , nM_out, nM_out , 'Mu'   , ErrStat2 , ErrMsg2); if(Failed()) return
  +      CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, PhiL, MLL, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, Temp, PhiL, 0.0_FeKi, MU  , ErrStat2, ErrMsg2); if(Failed()) return
         DEALLOCATE(Temp)
  -      ! PhiL = MATMUL( PhiL, MU2 )  !this is the nondimensionalization (MU2 is diagonal)   
  +      ! PhiL = MATMUL( PhiL, MU ) ! this is the normalization (MU is diagonal)   
         DO I = 1, nM_out
            PhiL(:,I) = PhiL(:,I) / SQRT( MU(I, I) )
         ENDDO    
         DEALLOCATE(MU)
         if (present(CRR)) then
            ! CB damping CMM = PhiL^T  CLL PhiL
  -         CALL AllocAry( Temp , nM, nL , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
  -         CALL AllocAry( MU   , nM, nL , 'Mu'   , ErrStat2 , ErrMsg2); if(Failed()) return
  -         MU  = TRANSPOSE(PhiL(1:nL, 1:nM))
  -         !Temp = MATMUL( MU, CLL )
  -         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, MU, CLL, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  -         !   CMM = MATMUL( Temp, PhiL(1:nL, 1:nM) )
  -         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, Temp, PhiL, 0.0_FeKi, CMM  , ErrStat2, ErrMsg2); if(Failed()) return
  -         DEALLOCATE(MU)
  +         CALL AllocAry( Temp , nM, nL    , 'Temp' , ErrStat2 , ErrMsg2); if(Failed()) return
  +         CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, PhiL(1:nL, 1:nM), CLL,  0.0_FeKi, Temp , ErrStat2, ErrMsg2); if(Failed()) return
  +         CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, Temp, PhiL(1:nL, 1:nM), 0.0_FeKi, CMM  , ErrStat2, ErrMsg2); if(Failed()) return
            DEALLOCATE(Temp)
         endif
      else
  
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
  