From 05393d34de4fbdc896c321cfc2dcc4736dde3081 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 27 Mar 2009 17:21:46 +0000
Subject: [PATCH] added invert()
---
Numeric/GmshMatrix.cpp | 36 ++++++++++++++++++++++++++++++++++++
Numeric/GmshMatrix.h | 8 ++++++++
2 files changed, 44 insertions(+)
diff --git a/Numeric/GmshMatrix.cpp b/Numeric/GmshMatrix.cpp
index 16bf2c8de2..75ad9f134e 100644
--- a/Numeric/GmshMatrix.cpp
+++ b/Numeric/GmshMatrix.cpp
@@ -99,6 +99,8 @@ extern "C" {
void dgesv_(int *N, int *nrhs, double *A, int *lda, int *ipiv,
double *b, int *ldb, int *info);
void dgetrf_(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
+ void dgetri_(int *M, double *A, int *lda, int *ipiv, double *work, int *lwork,
+ int *info);
void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
double *A, int *lda, double *S, double* U, int *ldu,
double *VT, int *ldvt, double *work, int *lwork, int *info);
@@ -120,6 +122,40 @@ bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<doub
return false;
}
+template<>
+bool gmshMatrix<double>::invert(gmshMatrix<double> &result)
+{
+ int M = size1(), N = size2(), lda = size1(), info;
+ int *ipiv = new int[std::min(M, N)];
+ bool ret = true;
+ result = *this;
+ dgetrf_(&M, &N, result._data, &lda, ipiv, &info);
+ if(info == 0){
+ int lwork = M * 4;
+ double *work = new double[lwork];
+ dgetri_(&M, result._data, &lda, ipiv, work, &lwork, &info);
+ if(info > 0){
+ Msg::Error("U(%d, %d)=0 in LU decomposition", info, info);
+ ret = false;
+ }
+ else if(info < 0){
+ Msg::Error("Wrong %d-th argument in matrix inversion", -info);
+ ret = false;
+ }
+ delete [] work;
+ }
+ else if(info > 0){
+ Msg::Error("U(%d, %d)=0 in LU decomposition", info, info);
+ ret = false;
+ }
+ else{
+ Msg::Error("Wrong %d-th argument in matrix factorization", -info);
+ ret = false;
+ }
+ delete [] ipiv;
+ return ret;
+}
+
template<>
double gmshMatrix<double>::determinant() const
{
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index e78be16628..40df2676e8 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -173,6 +173,14 @@ class gmshMatrix
Msg::Error("LU factorization requires LAPACK");
return false;
}
+#endif
+ ;
+ bool invert(gmshMatrix<scalar> &result)
+#if !defined(HAVE_LAPACK)
+ {
+ Msg::Error("LU factorization requires LAPACK");
+ return false;
+ }
#endif
;
gmshMatrix<scalar> cofactor(int i, int j) const
--
GitLab