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:       if(Failed()) return;
     endif
     if (nM_out<nM) then
        ErrMsg2='Cannot request output modes than modes.'; ErrStat2=ErrID_Fatal; 
        if(Failed()) return;
     endif
     
  
negative.code: !!
  !! NOTE: generic code
  SUBROUTINE CraigBamptonReduction(MM, KK, IDR, nR, IDL, nL, nM, nM_Out, MBB, MBM, KBB, PhiL, PhiR, OmegaL, ErrStat, ErrMsg, FG, FGR, FGL, FGB, FGM, CC, CBB, CBM, CMM) 
     use NWTC_LAPACK, only: LAPACK_GEMV
     REAL(FEKi),             INTENT(IN   ) :: MM(:, :) !< Mass matrix
     REAL(FEKi),             INTENT(IN   ) :: KK(:, :) !< Stiffness matrix
     INTEGER(IntKi),         INTENT(IN   ) :: nR
  