diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index eb75fa485ae4dc541ba85b3c69d07fb37ddd2521..58abf9f24acedef0ec3e4a1410a495207efbaa46 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -208,6 +208,8 @@ extern "C" {
   void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv,
                       double *b, int *ldb, int *info);
   void F77NAME(dgetrf)(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
+  void F77NAME(dgetrs)(char *trans, int *N, int *nrhs, double *A, int *lda, int *ipiv,
+                       double *b, int *ldb, int *info);
   void F77NAME(dgetri)(int *M, double *A, int *lda, int *ipiv, double *work,
                        int *lwork, int *info);
   void F77NAME(dgesvd)(const char* jobu, const char *jobvt, int *M, int *N,
@@ -284,10 +286,29 @@ bool fullMatrix<double>::luSolve(const fullVector<double> &rhs, fullVector<doubl
   F77NAME(dgesv)(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
   delete [] ipiv;
   if(info == 0) return true;
-  //  if(info > 0)
-  //    Msg::Error("U(%d,%d)=0 in LU decomposition", info, info);
-  //  else
-  //    Msg::Error("Wrong %d-th argument in LU decomposition", -info);
+  return false;
+}
+
+template<>
+bool fullMatrix<double>::luFactor()
+{
+  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;
+  if(info == 0) return true;
+  return false;
+}
+
+template<>
+bool fullMatrix<double>::luSubstitute(const fullVector<double> &rhs, fullVector<double> &result)
+{
+  int N = size1(), nrhs=1, lda = N, ldb = N, info;
+  int *ipiv = new int[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;
+  if(info == 0) return true;
   return false;
 }
 
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 03fc23a71b7ba5a6f77181d330795eb3a9df0fca..74cf557c339c2bf627cb68157b415dc9dbde60a5 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -679,11 +679,27 @@ class fullMatrix
       }
   }
   bool luSolve(const fullVector<scalar> &rhs, fullVector<scalar> &result)
+#if !defined(HAVE_LAPACK)
+  {
+    Msg::Error("LU factorization and substitution requires LAPACK");
+    return false;
+  }
+#endif
+  ;
+ bool luFactor()
 #if !defined(HAVE_LAPACK)
   {
     Msg::Error("LU factorization requires LAPACK");
     return false;
   }
+#endif
+  ;
+ bool luSubstitute(const fullVector<double> &rhs, fullVector<double> &result)
+#if !defined(HAVE_LAPACK)
+  {
+    Msg::Error("LU substitution requires LAPACK");
+    return false;
+  }
 #endif
   ;
   bool invertInPlace()