From 1c191c73b89507ff2e400d61be78f5af5164902f Mon Sep 17 00:00:00 2001
From: Sebastien Blaise <sebastien.blaise@uclouvain.be>
Date: Tue, 5 Nov 2013 13:33:35 +0000
Subject: [PATCH] lu factorization: adding pivot indices as parameter
---
Numeric/fullMatrix.cpp | 14 ++++++--------
Numeric/fullMatrix.h | 4 ++--
2 files changed, 8 insertions(+), 10 deletions(-)
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 58abf9f24a..b4325e0ef6 100644
--- a/Numeric/fullMatrix.cpp
+++ b/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;
}
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 74cf557c33..74debb530e 100644
--- a/Numeric/fullMatrix.h
+++ b/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");
--
GitLab