lu factorization: adding pivot indices as parameter

Numeric/fullMatrix.cpp

@@ -290,24 +290,22 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
 }
 
 template<>
-bool fullMatrix<double>::luFactor()
+bool fullMatrix<double>::luFactor(fullVector<int> &ipiv)
 {
   int M = size1(), N = size2(), lda = size1(), info;
-  int *ipiv = new int[std::min(M, N)];
-  F77NAME(dgetrf)(&M, &N, _data, &lda, ipiv, &info);
-  delete [] ipiv;
+  ipiv.resize(std::min(M, N));
+  F77NAME(dgetrf)(&M, &N, _data, &lda, &ipiv(0), &info);
   if(info == 0) return true;
   return false;
 }
 
 template<>
-bool fullMatrix<double>::luSubstitute(const fullVector<double> &rhs, fullVector<double> &result)
+bool fullMatrix<double>::luSubstitute(const fullVector<double> &rhs, fullVector<int> &ipiv, fullVector<double> &result)
 {
   int N = size1(), nrhs=1, lda = N, ldb = N, info;
-  int *ipiv = new int[N];
+  char trans = 'N';
   for(int i = 0; i < N; i++) result(i) = rhs(i);
-  F77NAME(dgetrs)("N", &N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
-  delete [] ipiv;
+  F77NAME(dgetrs)(&trans, &N, &nrhs, _data, &lda, &ipiv(0), result._data, &ldb, &info);
   if(info == 0) return true;
   return false;
 }

Numeric/fullMatrix.h

@@ -686,7 +686,7 @@ class fullMatrix
   }
 #endif
   ;
- bool luFactor()
+ bool luFactor(fullVector<int> &ipiv)
 #if !defined(HAVE_LAPACK)
   {
     Msg::Error("LU factorization requires LAPACK");
@@ -694,7 +694,7 @@ class fullMatrix
   }
 #endif
   ;
- bool luSubstitute(const fullVector<double> &rhs, fullVector<double> &result)
+ bool luSubstitute(const fullVector<scalar> &rhs, fullVector<int> &ipiv, fullVector<scalar> &result)
 #if !defined(HAVE_LAPACK)
   {
     Msg::Error("LU substitution requires LAPACK");