anchor.diff:        CALL AllocAry( Init%FG,     nDOF,             'Init%FG'  ,  ErrStat2, ErrMsg2); if(Failed()) return; ! system gravity force vector 
         CALL AllocAry( Temp   ,size(MM,1), nDOF,      'Temp'     ,  ErrStat2, ErrMsg2); if(Failed()) return; 
         CALL AllocAry( p%T_red_T,nDOF   , size(MM,1), 'T_red_T' ,  ErrStat2, ErrMsg2); if(Failed()) return; 
  -      ! Elimination
  +      ! --- Elimination (stack expensive)
         !Init%M  = matmul(transpose(p%T_red), matmul(MM, p%T_red))
         !Init%K  = matmul(transpose(p%T_red), matmul(KK, p%T_red))
  -      p%T_red_T = transpose(p%T_red)
  -      Temp    = matmul(MM, p%T_red)
  -      Init%M  = matmul(p%T_red_T, Temp)
  -      Temp    = matmul(KK, p%T_red)
  -      Init%K  = matmul(p%T_red_T, Temp)
  +      !p%T_red_T = transpose(p%T_red)
  +      do i = 1, size(p%T_red,1)
  +         do j = 1, size(p%T_red,2)
  +            p%T_red_T(j,i) = p%T_red(i,j)
  +         enddo
  +      enddo
  +      !Temp    = matmul(MM, p%T_red)
  +      CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, MM     , p%T_red, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  +      !Init%M  = matmul(p%T_red_T, Temp)
  +      CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, p%T_red, Temp   , 0.0_FeKi, Init%M, ErrStat2, ErrMsg2); if(Failed()) return
  +      !Temp    = matmul(KK, p%T_red)
  +      CALL LAPACK_gemm( 'N', 'N', 1.0_FeKi, KK     , p%T_red, 0.0_FeKi, Temp  , ErrStat2, ErrMsg2); if(Failed()) return
  +      !Init%K  = matmul(p%T_red_T, Temp)
  +      CALL LAPACK_gemm( 'T', 'N', 1.0_FeKi, p%T_red, Temp   , 0.0_FeKi, Init%K, ErrStat2, ErrMsg2); if(Failed()) return
  +      if (allocated(Temp))    deallocate(Temp)
         Init%FG = matmul(p%T_red_T, FF)
  -      if (allocated(Temp))    deallocate(Temp)
      endif
      !CALL AllocAry( Init%D,      nDOF, nDOF,  'Init%D'   ,  ErrStat2, ErrMsg2); if(Failed()) return; ! system damping matrix 
      !Init%D = 0 !< Used for additional damping 
  
positive.code:    real(FEKi), dimension(:,:), allocatable :: MM, KK
     real(FEKi), dimension(:),   allocatable :: FF
     real(FEKi), dimension(:,:), allocatable :: Temp
     integer(IntKi) :: nDOF, iDOF, nDOFPerNode, iNode, iiDOF
     ErrStat = ErrID_None
     ErrMsg  = ""
  
  
negative.code: !------------------------------------------------------------------------------------------------------
  !> Assemble stiffness and mass matrix, and gravity force vector
  SUBROUTINE DirectElimination(Init, p, ErrStat, ErrMsg)
     use NWTC_LAPACK, only: LAPACK_GEMM
     use IntegerList, only: len
     TYPE(SD_InitType),            INTENT(INOUT) :: Init
     TYPE(SD_ParameterType),target,INTENT(INOUT) :: p
  