[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Clean up
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Clean up |
|
Date: |
Wed, 18 Oct 2023 03:59:12 -0400 |
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch master
in repository getfem.
The following commit(s) were added to refs/heads/master by this push:
new f50e5995 Clean up
f50e5995 is described below
commit f50e5995c0d69f0559fc9f75f8e96e39cfbe6949
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Wed Oct 18 09:58:59 2023 +0200
Clean up
- remove unused preprocessor macros
- whitespace
- include local gmm_dense_lu.h header file
---
contrib/crack_plate/crack_bilaplacian_sif.cc | 2 +-
src/dal_backtrace.cc | 17 +-
src/getfem/getfem_crack_sif.h | 2 +-
src/getfem_superlu.cc | 167 ++++++++-------
src/gmm/gmm_blas_interface.h | 3 +-
src/gmm/gmm_except.h | 2 +-
src/gmm/gmm_opt.h | 2 +-
src/gmm/gmm_superlu_interface.h | 290 +++++++++++++--------------
8 files changed, 241 insertions(+), 244 deletions(-)
diff --git a/contrib/crack_plate/crack_bilaplacian_sif.cc
b/contrib/crack_plate/crack_bilaplacian_sif.cc
index 3da6139f..6dda5e0d 100644
--- a/contrib/crack_plate/crack_bilaplacian_sif.cc
+++ b/contrib/crack_plate/crack_bilaplacian_sif.cc
@@ -59,7 +59,7 @@ namespace getfem {
void get_crack_tip_and_orientation_KL(const level_set &/* ls */,
base_node &P,
base_small_vector &T, base_small_vector
&N) {
- cerr << __PRETTY_FUNCTION__ << " IS TO BE DONE\n";
+ cerr << GMM_PRETTY_FUNCTION << " IS TO BE DONE\n";
/* too lazy to do it now */
P.resize(2); P[0] = 0.; P[1] = 0;
T.resize(2); T[0] = 1; T[1] = 0;
diff --git a/src/dal_backtrace.cc b/src/dal_backtrace.cc
index 2f711ffc..464d57ff 100644
--- a/src/dal_backtrace.cc
+++ b/src/dal_backtrace.cc
@@ -32,25 +32,24 @@ namespace dal {
If you call this function with a non-mangled name (i.e. "main"),
you will get strange results.
*/
- std::string demangle(const char *
-#ifdef GETFEM_HAVE_CXXABI_H
- s
-#endif
- ) {
#ifdef GETFEM_HAVE_CXXABI_H
+ std::string demangle(const char *s) {
int status;
/* documented in
http://gcc.gnu.org/onlinedocs/libstdc++/latest-doxygen/namespaceabi.html */
char *sd = abi::__cxa_demangle(s, NULL, NULL, &status);
if (sd == NULL || status) {
- if (sd) free(sd);
+ if (sd)
+ free(sd);
return std::string(""); // + " [could not be demangled]";
} else {
- std::string res(sd); free(sd); return res;
+ std::string res(sd);
+ free(sd);
+ return res;
}
+ }
#else
- return std::string("");
+ std::string demangle(const char *) { return std::string(""); }
#endif
- }
#ifdef GETFEM_HAVE_BACKTRACE
void dump_glibc_backtrace() {
diff --git a/src/getfem/getfem_crack_sif.h b/src/getfem/getfem_crack_sif.h
index 6d849ba9..4f33d581 100644
--- a/src/getfem/getfem_crack_sif.h
+++ b/src/getfem/getfem_crack_sif.h
@@ -73,7 +73,7 @@ namespace getfem {
inline void get_crack_tip_and_orientation(const level_set &/* ls */,
base_node &P,
base_small_vector &T, base_small_vector
&N) {
- cerr << __PRETTY_FUNCTION__ << " IS TO BE DONE\n";
+ cerr << GMM_PRETTY_FUNCTION << " IS TO BE DONE\n";
/* too lazy to do it now */
P.resize(2); P[0] = .5; P[1] = 0;
T.resize(2); T[0] = 1; T[1] = 0;
diff --git a/src/getfem_superlu.cc b/src/getfem_superlu.cc
index be113dd5..4e7c3538 100644
--- a/src/getfem_superlu.cc
+++ b/src/getfem_superlu.cc
@@ -45,7 +45,7 @@ namespace SuperLU_C {
#include "superlu/slu_cdefs.h"
}
namespace SuperLU_Z {
-#include "superlu/slu_zdefs.h"
+#include "superlu/slu_zdefs.h"
}
@@ -55,28 +55,28 @@ namespace gmm {
/* interface for Create_CompCol_Matrix */
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- float *a, int *ir, int *jc) {
+ float *a, int *ir, int *jc) {
SuperLU_S::sCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
- SLU_NC, SLU_S, SLU_GE);
+ SLU_NC, SLU_S, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- double *a, int *ir, int *jc) {
+ double *a, int *ir, int *jc) {
SuperLU_D::dCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
- SLU_NC, SLU_D, SLU_GE);
+ SLU_NC, SLU_D, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- std::complex<float> *a, int *ir, int *jc) {
+ std::complex<float> *a, int *ir, int *jc) {
SuperLU_C::cCreate_CompCol_Matrix(A, m, n, nnz, (SuperLU_C::complex *)(a),
- ir, jc, SLU_NC, SLU_C, SLU_GE);
+ ir, jc, SLU_NC, SLU_C, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- std::complex<double> *a, int *ir, int *jc) {
+ std::complex<double> *a, int *ir, int *jc) {
SuperLU_Z::zCreate_CompCol_Matrix(A, m, n, nnz,
- (SuperLU_Z::doublecomplex *)(a), ir, jc,
- SLU_NC, SLU_Z, SLU_GE);
+ (SuperLU_Z::doublecomplex *)(a), ir, jc,
+ SLU_NC, SLU_Z, SLU_GE);
}
/* interface for Create_Dense_Matrix */
@@ -86,14 +86,14 @@ namespace gmm {
inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, double *a, int
k)
{ SuperLU_D::dCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_D, SLU_GE); }
inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
- std::complex<float> *a, int k) {
+ std::complex<float> *a, int k) {
SuperLU_C::cCreate_Dense_Matrix(A, m, n, (SuperLU_C::complex *)(a),
- k, SLU_DN, SLU_C, SLU_GE);
+ k, SLU_DN, SLU_C, SLU_GE);
}
- inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
- std::complex<double> *a, int k) {
+ inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
+ std::complex<double> *a, int k) {
SuperLU_Z::zCreate_Dense_Matrix(A, m, n, (SuperLU_Z::doublecomplex *)(a),
- k, SLU_DN, SLU_Z, SLU_GE);
+ k, SLU_DN, SLU_Z, SLU_GE);
}
/* interface for gssv */
@@ -113,18 +113,18 @@ namespace gmm {
/* interface for gssvx */
#define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
- inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,
\
- int *perm_c, int *perm_r, int *etree, char *equed, \
- FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
- SuperMatrix *U, void *work, int lwork, \
- SuperMatrix *B, SuperMatrix *X, \
- FLOATTYPE *recip_pivot_growth, \
- FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
- SuperLUStat_t *stats, int *info, KEYTYPE) { \
+ inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,
\
+ int *perm_c, int *perm_r, int *etree, char
*equed, \
+ FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,
\
+ SuperMatrix *U, void *work, int lwork,
\
+ SuperMatrix *B, SuperMatrix *X,
\
+ FLOATTYPE *recip_pivot_growth,
\
+ FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE
*berr, \
+ SuperLUStat_t *stats, int *info, KEYTYPE) {
\
NAMESPACE::mem_usage_t mem_usage; \
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
- U, work, lwork, B, X, recip_pivot_growth, rcond, \
- ferr, berr, &mem_usage, stats, info); \
+ U, work, lwork, B, X, recip_pivot_growth, rcond, \
+ ferr, berr, &mem_usage, stats, info); \
return mem_usage.for_lu; /* bytes used by the factor storage */ \
}
@@ -139,10 +139,10 @@ namespace gmm {
template<typename T>
int SuperLU_solve(const gmm::csc_matrix<T> &csc_A, T *sol, T *rhs,
- double& rcond_, int permc_spec) {
+ double& rcond_, int permc_spec) {
/*
* Get column permutation vector perm_c[], according to permc_spec:
- * permc_spec = 0: use the natural ordering
+ * permc_spec = 0: use the natural ordering
* permc_spec = 1: use minimum degree ordering on structure of A'*A
* permc_spec = 2: use minimum degree ordering on structure of A'+A
* permc_spec = 3: use approximate minimum degree column ordering
@@ -157,7 +157,7 @@ namespace gmm {
if ((2 * nz / n) >= m)
GMM_WARNING2("CAUTION : it seems that SuperLU has a problem"
- " for nearly dense sparse matrices");
+ " for nearly dense sparse matrices");
superlu_options_t options;
set_default_options(&options);
@@ -174,11 +174,10 @@ namespace gmm {
SuperMatrix SA, SL, SU, SB, SX; // SuperLU format.
Create_CompCol_Matrix(&SA, m, n, nz, const_cast<T*>(&csc_A.pr[0]),
- const_cast<int *>((const int *)(&csc_A.ir[0])),
+ const_cast<int *>((const int *)(&csc_A.ir[0])),
const_cast<int *>((const int *)(&csc_A.jc[0])));
Create_Dense_Matrix(&SB, m, nrhs, &rhs[0], m);
Create_Dense_Matrix(&SX, m, nrhs, &sol[0], m);
-
memset(&SL,0,sizeof SL);
memset(&SU,0,sizeof SU);
@@ -189,21 +188,21 @@ namespace gmm {
R recip_pivot_gross, rcond;
std::vector<int> perm_r(m), perm_c(n);
- SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, equed /* output */,
- &Rscale[0] /* row scale factors (output) */,
- &Cscale[0] /* col scale factors (output) */,
- &SL /* fact L (output)*/, &SU /* fact U (output)*/,
- NULL /* work */,
- 0 /* lwork: superlu auto allocates (input) */,
- &SB /* rhs */, &SX /* solution */,
- &recip_pivot_gross /* reciprocal pivot growth */
- /* factor max_j( norm(A_j)/norm(U_j) ). */,
- &rcond /*estimate of the reciprocal condition */
- /* number of the matrix A after equilibration */,
- &ferr[0] /* estimated forward error */,
- &berr[0] /* relative backward error */,
- &stat, &info, T());
+ SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, equed /* output */,
+ &Rscale[0] /* row scale factors (output) */,
+ &Cscale[0] /* col scale factors (output) */,
+ &SL /* fact L (output)*/, &SU /* fact U (output)*/,
+ NULL /* work */,
+ 0 /* lwork: superlu auto allocates (input) */,
+ &SB /* rhs */, &SX /* solution */,
+ &recip_pivot_gross /* reciprocal pivot growth */
+ /* factor max_j( norm(A_j)/norm(U_j) ). */,
+ &rcond /*estimate of the reciprocal condition */
+ /* number of the matrix A after equilibration */,
+ &ferr[0] /* estimated forward error */,
+ &berr[0] /* relative backward error */,
+ &stat, &info, T());
rcond_ = rcond;
if (SB.Store) Destroy_SuperMatrix_Store(&SB);
if (SX.Store) Destroy_SuperMatrix_Store(&SX);
@@ -232,17 +231,17 @@ namespace gmm {
mutable char equed;
void free_supermatrix() {
if (is_init) {
- if (SB.Store) Destroy_SuperMatrix_Store(&SB);
- if (SX.Store) Destroy_SuperMatrix_Store(&SX);
- if (SA.Store) Destroy_SuperMatrix_Store(&SA);
- if (SL.Store) Destroy_SuperNode_Matrix(&SL);
- if (SU.Store) Destroy_CompCol_Matrix(&SU);
+ if (SB.Store) Destroy_SuperMatrix_Store(&SB);
+ if (SX.Store) Destroy_SuperMatrix_Store(&SX);
+ if (SA.Store) Destroy_SuperMatrix_Store(&SA);
+ if (SL.Store) Destroy_SuperNode_Matrix(&SL);
+ if (SU.Store) Destroy_CompCol_Matrix(&SU);
}
}
SuperLU_factor_impl_common() : is_init(false) {}
virtual ~SuperLU_factor_impl_common() { free_supermatrix(); }
};
-
+
template <typename T> struct SuperLU_factor_impl : public
SuperLU_factor_impl_common {
typedef typename gmm::number_traits<T>::magnitude_type R;
@@ -259,7 +258,7 @@ namespace gmm {
void SuperLU_factor_impl<T>::build_with(const gmm::csc_matrix<T> &A, int
permc_spec) {
/*
* Get column permutation vector perm_c[], according to permc_spec:
- * permc_spec = 0: use the natural ordering
+ * permc_spec = 0: use the natural ordering
* permc_spec = 1: use minimum degree ordering on structure of A'*A
* permc_spec = 2: use minimum degree ordering on structure of A'+A
* permc_spec = 3: use approximate minimum degree column ordering
@@ -273,7 +272,7 @@ namespace gmm {
GMM_ASSERT1(nz != 0, "Cannot factor a matrix full of zeros!");
GMM_ASSERT1(n == m, "Cannot factor a non-square matrix");
-
+
set_default_options(&options);
options.ColPerm = NATURAL;
options.PrintStat = NO;
@@ -284,7 +283,7 @@ namespace gmm {
case 3 : options.ColPerm = COLAMD; break;
}
StatInit(&stat);
-
+
Create_CompCol_Matrix(&SA, m, n, nz, const_cast<T*>(&A.pr[0]),
const_cast<int *>((const int *)(&A.ir[0])),
const_cast<int *>((const int *)(&A.jc[0])));
@@ -299,34 +298,34 @@ namespace gmm {
ferr.resize(1); berr.resize(1);
R recip_pivot_gross, rcond;
perm_r.resize(m); perm_c.resize(n);
- memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, &equed /* output
*/,
- &Rscale[0] /* row scale factors (output)
*/,
+ memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, &equed /* output
*/,
+ &Rscale[0] /* row scale factors (output)
*/,
&Cscale[0] /* col scale factors (output)
*/,
- &SL /* fact L (output)*/, &SU /* fact U
(output)*/,
- NULL /* work
*/,
- 0 /* lwork: superlu auto allocates (input)
*/,
+ &SL /* fact L (output)*/, &SU /* fact U
(output)*/,
+ NULL /* work
*/,
+ 0 /* lwork: superlu auto allocates (input)
*/,
&SB /* rhs */, &SX /* solution
*/,
&recip_pivot_gross /* reciprocal pivot growth
*/
- /* factor max_j( norm(A_j)/norm(U_j) ).
*/,
+ /* factor max_j( norm(A_j)/norm(U_j) ).
*/,
&rcond /*estimate of the reciprocal condition
*/
/* number of the matrix A after equilibration
*/,
&ferr[0] /* estimated forward error
*/,
&berr[0] /* relative backward error
*/,
&stat, &info, T());
-
+
Destroy_SuperMatrix_Store(&SB);
Destroy_SuperMatrix_Store(&SX);
Create_Dense_Matrix(&SB, m, 1, &rhs[0], m);
Create_Dense_Matrix(&SX, m, 1, &sol[0], m);
StatFree(&stat);
-
+
GMM_ASSERT1(info != -333333333, "SuperLU was cancelled.");
GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
is_init = true;
}
- template <typename T>
+ template <typename T>
void SuperLU_factor_impl<T>::solve(int transp) {
options.Fact = FACTORED;
options.IterRefine = NOREFINE;
@@ -339,16 +338,16 @@ namespace gmm {
StatInit(&stat);
int info = 0;
R recip_pivot_gross, rcond;
- SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, &equed /* output */,
- &Rscale[0] /* row scale factors (output) */,
+ SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, &equed /* output */,
+ &Rscale[0] /* row scale factors (output) */,
&Cscale[0] /* col scale factors (output) */,
- &SL /* fact L (output)*/, &SU /* fact U (output)*/,
- NULL /* work */,
- 0 /* lwork: superlu auto allocates (input) */,
+ &SL /* fact L (output)*/, &SU /* fact U (output)*/,
+ NULL /* work */,
+ 0 /* lwork: superlu auto allocates (input) */,
&SB /* rhs */, &SX /* solution */,
&recip_pivot_gross /* reciprocal pivot growth */
- /* factor max_j( norm(A_j)/norm(U_j) ). */,
+ /* factor max_j( norm(A_j)/norm(U_j) ). */,
&rcond /*estimate of the reciprocal condition */
/* number of the matrix A after equilibration */,
&ferr[0] /* estimated forward error */,
@@ -357,8 +356,8 @@ namespace gmm {
StatFree(&stat);
GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
}
-
- template<typename T> void
+
+ template<typename T> void
SuperLU_factor<T>::build_with(const gmm::csc_matrix<T> &A, int permc_spec) {
((SuperLU_factor_impl<T>*)impl.get())->build_with(A,permc_spec);
}
@@ -378,27 +377,27 @@ namespace gmm {
return ((SuperLU_factor_impl<T>*)impl.get())->rhs;
}
- template<typename T>
+ template<typename T>
SuperLU_factor<T>::SuperLU_factor() {
impl = std::make_shared<SuperLU_factor_impl<T>>();
}
- template<typename T>
+ template<typename T>
SuperLU_factor<T>::SuperLU_factor(const SuperLU_factor& other) {
impl = std::make_shared<SuperLU_factor_impl<T>>();
GMM_ASSERT1(!(other.impl->is_init),
- "copy of initialized SuperLU_factor is forbidden");
+ "copy of initialized SuperLU_factor is forbidden");
other.impl->is_init = false;
}
- template<typename T> SuperLU_factor<T>&
+ template<typename T> SuperLU_factor<T>&
SuperLU_factor<T>::operator=(const SuperLU_factor& other) {
GMM_ASSERT1(!(other.impl->is_init) && !(impl->is_init),
- "assignment of initialized SuperLU_factor is forbidden");
+ "assignment of initialized SuperLU_factor is forbidden");
return *this;
}
- template<typename T> float
+ template<typename T> float
SuperLU_factor<T>::memsize() const {
return impl->memory_used;
}
@@ -409,8 +408,8 @@ namespace gmm {
SuperLU_factor<std::complex<float> > c;
SuperLU_factor<std::complex<double> > d;
//a = 0; b = 0; c = 0; d = 0;
- }
- */
+ }
+ */
}
template class gmm::SuperLU_factor<float>;
diff --git a/src/gmm/gmm_blas_interface.h b/src/gmm/gmm_blas_interface.h
index df5d7158..22351ab6 100644
--- a/src/gmm/gmm_blas_interface.h
+++ b/src/gmm/gmm_blas_interface.h
@@ -35,8 +35,7 @@
@brief gmm interface for fortran BLAS.
*/
-#if defined(GETFEM_USES_BLAS) || defined(GMM_USES_BLAS) \
- || defined(GMM_USES_LAPACK) || defined(GMM_USES_ATLAS)
+#if defined(GMM_USES_BLAS) || defined(GMM_USES_LAPACK)
#ifndef GMM_BLAS_INTERFACE_H
#define GMM_BLAS_INTERFACE_H
diff --git a/src/gmm/gmm_except.h b/src/gmm/gmm_except.h
index 3b38aa16..442de139 100644
--- a/src/gmm/gmm_except.h
+++ b/src/gmm/gmm_except.h
@@ -63,7 +63,7 @@ namespace gmm {
int errorLevel_;
};
-#ifdef GETFEM_HAVE_PRETTY_FUNCTION
+#ifdef GMM_HAVE_PRETTY_FUNCTION
# define GMM_PRETTY_FUNCTION __PRETTY_FUNCTION__
#elif _MSC_VER
# define GMM_PRETTY_FUNCTION __FUNCTION__
diff --git a/src/gmm/gmm_opt.h b/src/gmm/gmm_opt.h
index 7be45c44..1a4b1f14 100644
--- a/src/gmm/gmm_opt.h
+++ b/src/gmm/gmm_opt.h
@@ -37,7 +37,7 @@
#ifndef GMM_OPT_H__
#define GMM_OPT_H__
-#include <gmm/gmm_dense_lu.h>
+#include "gmm_dense_lu.h"
namespace gmm {
diff --git a/src/gmm/gmm_superlu_interface.h b/src/gmm/gmm_superlu_interface.h
index 2a7a2a4a..9605dc65 100644
--- a/src/gmm/gmm_superlu_interface.h
+++ b/src/gmm/gmm_superlu_interface.h
@@ -65,7 +65,7 @@ namespace SuperLU_C {
#include "superlu/slu_cdefs.h"
}
namespace SuperLU_Z {
-#include "superlu/slu_zdefs.h"
+#include "superlu/slu_zdefs.h"
}
@@ -75,28 +75,28 @@ namespace gmm {
/* interface for Create_CompCol_Matrix */
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- float *a, int *ir, int *jc) {
+ float *a, int *ir, int *jc) {
SuperLU_S::sCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
- SLU_NC, SLU_S, SLU_GE);
+ SLU_NC, SLU_S, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- double *a, int *ir, int *jc) {
+ double *a, int *ir, int *jc) {
SuperLU_D::dCreate_CompCol_Matrix(A, m, n, nnz, a, ir, jc,
- SLU_NC, SLU_D, SLU_GE);
+ SLU_NC, SLU_D, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- std::complex<float> *a, int *ir, int *jc) {
+ std::complex<float> *a, int *ir, int *jc) {
SuperLU_C::cCreate_CompCol_Matrix(A, m, n, nnz, (SuperLU_C::complex *)(a),
- ir, jc, SLU_NC, SLU_C, SLU_GE);
+ ir, jc, SLU_NC, SLU_C, SLU_GE);
}
-
+
inline void Create_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
- std::complex<double> *a, int *ir, int *jc) {
+ std::complex<double> *a, int *ir, int *jc) {
SuperLU_Z::zCreate_CompCol_Matrix(A, m, n, nnz,
- (SuperLU_Z::doublecomplex *)(a), ir, jc,
- SLU_NC, SLU_Z, SLU_GE);
+ (SuperLU_Z::doublecomplex *)(a), ir, jc,
+ SLU_NC, SLU_Z, SLU_GE);
}
/* interface for Create_Dense_Matrix */
@@ -106,14 +106,14 @@ namespace gmm {
inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n, double *a, int
k)
{ SuperLU_D::dCreate_Dense_Matrix(A, m, n, a, k, SLU_DN, SLU_D, SLU_GE); }
inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
- std::complex<float> *a, int k) {
+ std::complex<float> *a, int k) {
SuperLU_C::cCreate_Dense_Matrix(A, m, n, (SuperLU_C::complex *)(a),
- k, SLU_DN, SLU_C, SLU_GE);
+ k, SLU_DN, SLU_C, SLU_GE);
}
- inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
- std::complex<double> *a, int k) {
+ inline void Create_Dense_Matrix(SuperMatrix *A, int m, int n,
+ std::complex<double> *a, int k) {
SuperLU_Z::zCreate_Dense_Matrix(A, m, n, (SuperLU_Z::doublecomplex *)(a),
- k, SLU_DN, SLU_Z, SLU_GE);
+ k, SLU_DN, SLU_Z, SLU_GE);
}
/* interface for gssv */
@@ -133,18 +133,18 @@ namespace gmm {
/* interface for gssvx */
#define DECL_GSSVX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
- inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,
\
- int *perm_c, int *perm_r, int *etree, char *equed, \
- FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
- SuperMatrix *U, void *work, int lwork, \
- SuperMatrix *B, SuperMatrix *X, \
- FLOATTYPE *recip_pivot_growth, \
- FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
- SuperLUStat_t *stats, int *info, KEYTYPE) { \
+ inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,
\
+ int *perm_c, int *perm_r, int *etree, char
*equed, \
+ FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,
\
+ SuperMatrix *U, void *work, int lwork,
\
+ SuperMatrix *B, SuperMatrix *X,
\
+ FLOATTYPE *recip_pivot_growth,
\
+ FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE
*berr, \
+ SuperLUStat_t *stats, int *info, KEYTYPE) {
\
NAMESPACE::mem_usage_t mem_usage; \
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
- U, work, lwork, B, X, recip_pivot_growth, rcond, \
- ferr, berr, &mem_usage, stats, info); \
+ U, work, lwork, B, X, recip_pivot_growth, rcond, \
+ ferr, berr, &mem_usage, stats, info); \
return mem_usage.for_lu; /* bytes used by the factor storage */ \
}
@@ -159,11 +159,11 @@ namespace gmm {
template <typename MAT, typename VECTX, typename VECTB>
int SuperLU_solve(const MAT &A, const VECTX &X_, const VECTB &B,
- double& rcond_, int permc_spec = 3) {
+ double& rcond_, int permc_spec = 3) {
VECTX &X = const_cast<VECTX &>(X_);
/*
* Get column permutation vector perm_c[], according to permc_spec:
- * permc_spec = 0: use the natural ordering
+ * permc_spec = 0: use the natural ordering
* permc_spec = 1: use minimum degree ordering on structure of A'*A
* permc_spec = 2: use minimum degree ordering on structure of A'+A
* permc_spec = 3: use approximate minimum degree column ordering
@@ -180,7 +180,7 @@ namespace gmm {
int nz = nnz(csc_A);
if ((2 * nz / n) >= m)
GMM_WARNING2("CAUTION : it seems that SuperLU has a problem"
- " for nearly dense sparse matrices");
+ " for nearly dense sparse matrices");
superlu_options_t options;
set_default_options(&options);
@@ -197,7 +197,7 @@ namespace gmm {
SuperMatrix SA, SL, SU, SB, SX; // SuperLU format.
Create_CompCol_Matrix(&SA, m, n, nz, (double *)(&(csc_A.pr[0])),
- (int *)(&(csc_A.ir[0])), (int *)(&(csc_A.jc[0])));
+ (int *)(&(csc_A.ir[0])), (int *)(&(csc_A.jc[0])));
Create_Dense_Matrix(&SB, m, nrhs, &rhs[0], m);
Create_Dense_Matrix(&SX, m, nrhs, &sol[0], m);
memset(&SL,0,sizeof SL);
@@ -210,21 +210,21 @@ namespace gmm {
R recip_pivot_gross, rcond;
std::vector<int> perm_r(m), perm_c(n);
- SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, equed /* output */,
- &Rscale[0] /* row scale factors (output) */,
- &Cscale[0] /* col scale factors (output) */,
- &SL /* fact L (output)*/, &SU /* fact U (output)*/,
- NULL /* work */,
- 0 /* lwork: superlu auto allocates (input) */,
- &SB /* rhs */, &SX /* solution */,
- &recip_pivot_gross /* reciprocal pivot growth */
- /* factor max_j( norm(A_j)/norm(U_j) ). */,
- &rcond /*estimate of the reciprocal condition */
- /* number of the matrix A after equilibration */,
- &ferr[0] /* estimated forward error */,
- &berr[0] /* relative backward error */,
- &stat, &info, T());
+ SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, equed /* output */,
+ &Rscale[0] /* row scale factors (output) */,
+ &Cscale[0] /* col scale factors (output) */,
+ &SL /* fact L (output)*/, &SU /* fact U (output)*/,
+ NULL /* work */,
+ 0 /* lwork: superlu auto allocates (input) */,
+ &SB /* rhs */, &SX /* solution */,
+ &recip_pivot_gross /* reciprocal pivot growth */
+ /* factor max_j( norm(A_j)/norm(U_j) ). */,
+ &rcond /*estimate of the reciprocal condition */
+ /* number of the matrix A after equilibration */,
+ &ferr[0] /* estimated forward error */,
+ &berr[0] /* relative backward error */,
+ &stat, &info, T());
rcond_ = rcond;
Destroy_SuperMatrix_Store(&SB);
Destroy_SuperMatrix_Store(&SX);
@@ -258,7 +258,7 @@ namespace gmm {
enum { LU_NOTRANSP, LU_TRANSP, LU_CONJUGATED };
void free_supermatrix(void);
template <class MAT> void build_with(const MAT &A, int permc_spec = 3);
- template <typename VECTX, typename VECTB>
+ template <typename VECTX, typename VECTB>
/* transp = LU_NOTRANSP -> solves Ax = B
transp = LU_TRANSP -> solves A'x = B
transp = LU_CONJUGATED -> solves conj(A)X = B */
@@ -266,12 +266,12 @@ namespace gmm {
SuperLU_factor(void) { is_init = false; }
SuperLU_factor(const SuperLU_factor& other) {
GMM_ASSERT2(!(other.is_init),
- "copy of initialized SuperLU_factor is forbidden");
+ "copy of initialized SuperLU_factor is forbidden");
is_init = false;
}
SuperLU_factor& operator=(const SuperLU_factor& other) {
GMM_ASSERT2(!(other.is_init) && !is_init,
- "assignment of initialized SuperLU_factor is forbidden");
+ "assignment of initialized SuperLU_factor is forbidden");
return *this;
}
~SuperLU_factor() { free_supermatrix(); }
@@ -280,117 +280,117 @@ namespace gmm {
template <class T> void SuperLU_factor<T>::free_supermatrix(void) {
- if (is_init) {
- if (SB.Store) Destroy_SuperMatrix_Store(&SB);
- if (SX.Store) Destroy_SuperMatrix_Store(&SX);
- if (SA.Store) Destroy_SuperMatrix_Store(&SA);
- if (SL.Store) Destroy_SuperNode_Matrix(&SL);
- if (SU.Store) Destroy_CompCol_Matrix(&SU);
- }
+ if (is_init) {
+ if (SB.Store) Destroy_SuperMatrix_Store(&SB);
+ if (SX.Store) Destroy_SuperMatrix_Store(&SX);
+ if (SA.Store) Destroy_SuperMatrix_Store(&SA);
+ if (SL.Store) Destroy_SuperNode_Matrix(&SL);
+ if (SU.Store) Destroy_CompCol_Matrix(&SU);
}
+ }
-
- template <class T> template <class MAT>
- void SuperLU_factor<T>::build_with(const MAT &A, int permc_spec) {
+
+ template <class T> template <class MAT>
+ void SuperLU_factor<T>::build_with(const MAT &A, int permc_spec) {
/*
* Get column permutation vector perm_c[], according to permc_spec:
- * permc_spec = 0: use the natural ordering
+ * permc_spec = 0: use the natural ordering
* permc_spec = 1: use minimum degree ordering on structure of A'*A
* permc_spec = 2: use minimum degree ordering on structure of A'+A
* permc_spec = 3: use approximate minimum degree column ordering
*/
- free_supermatrix();
- int n = mat_nrows(A), m = mat_ncols(A), info = 0;
- csc_A.init_with(A);
-
- rhs.resize(m); sol.resize(m);
- gmm::clear(rhs);
- int nz = nnz(csc_A);
-
- set_default_options(&options);
- options.ColPerm = NATURAL;
- options.PrintStat = NO;
- options.ConditionNumber = NO;
- switch (permc_spec) {
+ free_supermatrix();
+ int n = mat_nrows(A), m = mat_ncols(A), info = 0;
+ csc_A.init_with(A);
+
+ rhs.resize(m); sol.resize(m);
+ gmm::clear(rhs);
+ int nz = nnz(csc_A);
+
+ set_default_options(&options);
+ options.ColPerm = NATURAL;
+ options.PrintStat = NO;
+ options.ConditionNumber = NO;
+ switch (permc_spec) {
case 1 : options.ColPerm = MMD_ATA; break;
case 2 : options.ColPerm = MMD_AT_PLUS_A; break;
case 3 : options.ColPerm = COLAMD; break;
- }
- StatInit(&stat);
-
- Create_CompCol_Matrix(&SA, m, n, nz, (double *)(&(csc_A.pr[0])),
- (int *)(&(csc_A.ir[0])), (int *)(&(csc_A.jc[0])));
-
- Create_Dense_Matrix(&SB, m, 0, &rhs[0], m);
- Create_Dense_Matrix(&SX, m, 0, &sol[0], m);
- memset(&SL,0,sizeof SL);
- memset(&SU,0,sizeof SU);
- equed = 'B';
- Rscale.resize(m); Cscale.resize(n); etree.resize(n);
- ferr.resize(1); berr.resize(1);
- R recip_pivot_gross, rcond;
- perm_r.resize(m); perm_c.resize(n);
- memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, &equed /* output */,
- &Rscale[0] /* row scale factors (output) */,
- &Cscale[0] /* col scale factors (output) */,
- &SL /* fact L (output)*/, &SU /* fact U (output)*/,
- NULL /* work */,
- 0 /* lwork: superlu auto allocates (input) */,
- &SB /* rhs */, &SX /* solution */,
- &recip_pivot_gross /* reciprocal pivot growth */
- /* factor max_j( norm(A_j)/norm(U_j) ). */,
- &rcond /*estimate of the reciprocal condition */
- /* number of the matrix A after equilibration */,
- &ferr[0] /* estimated forward error */,
- &berr[0] /* relative backward error */,
- &stat, &info, T());
-
- Destroy_SuperMatrix_Store(&SB);
- Destroy_SuperMatrix_Store(&SX);
- Create_Dense_Matrix(&SB, m, 1, &rhs[0], m);
- Create_Dense_Matrix(&SX, m, 1, &sol[0], m);
- StatFree(&stat);
-
- GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
- is_init = true;
}
-
- template <class T> template <typename VECTX, typename VECTB>
- void SuperLU_factor<T>::solve(const VECTX &X_, const VECTB &B,
- int transp) const {
- VECTX &X = const_cast<VECTX &>(X_);
- gmm::copy(B, rhs);
- options.Fact = FACTORED;
- options.IterRefine = NOREFINE;
- switch (transp) {
+ StatInit(&stat);
+
+ Create_CompCol_Matrix(&SA, m, n, nz, (double *)(&(csc_A.pr[0])),
+ (int *)(&(csc_A.ir[0])), (int *)(&(csc_A.jc[0])));
+
+ Create_Dense_Matrix(&SB, m, 0, &rhs[0], m);
+ Create_Dense_Matrix(&SX, m, 0, &sol[0], m);
+ memset(&SL,0,sizeof SL);
+ memset(&SU,0,sizeof SU);
+ equed = 'B';
+ Rscale.resize(m); Cscale.resize(n); etree.resize(n);
+ ferr.resize(1); berr.resize(1);
+ R recip_pivot_gross, rcond;
+ perm_r.resize(m); perm_c.resize(n);
+ memory_used = SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, &equed /* output
*/,
+ &Rscale[0] /* row scale factors (output)
*/,
+ &Cscale[0] /* col scale factors (output)
*/,
+ &SL /* fact L (output)*/, &SU /* fact U
(output)*/,
+ NULL /* work
*/,
+ 0 /* lwork: superlu auto allocates (input)
*/,
+ &SB /* rhs */, &SX /* solution
*/,
+ &recip_pivot_gross /* reciprocal pivot growth
*/
+ /* factor max_j( norm(A_j)/norm(U_j) ).
*/,
+ &rcond /*estimate of the reciprocal condition
*/
+ /* number of the matrix A after equilibration
*/,
+ &ferr[0] /* estimated forward error
*/,
+ &berr[0] /* relative backward error
*/,
+ &stat, &info, T());
+
+ Destroy_SuperMatrix_Store(&SB);
+ Destroy_SuperMatrix_Store(&SX);
+ Create_Dense_Matrix(&SB, m, 1, &rhs[0], m);
+ Create_Dense_Matrix(&SX, m, 1, &sol[0], m);
+ StatFree(&stat);
+
+ GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
+ is_init = true;
+ }
+
+ template <class T> template <typename VECTX, typename VECTB>
+ void SuperLU_factor<T>::solve(const VECTX &X_, const VECTB &B,
+ int transp) const {
+ VECTX &X = const_cast<VECTX &>(X_);
+ gmm::copy(B, rhs);
+ options.Fact = FACTORED;
+ options.IterRefine = NOREFINE;
+ switch (transp) {
case LU_NOTRANSP: options.Trans = NOTRANS; break;
case LU_TRANSP: options.Trans = TRANS; break;
case LU_CONJUGATED: options.Trans = CONJ; break;
default: GMM_ASSERT1(false, "invalid value for transposition option");
- }
- StatInit(&stat);
- int info = 0;
- R recip_pivot_gross, rcond;
- SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
- &etree[0] /* output */, &equed /* output */,
- &Rscale[0] /* row scale factors (output) */,
- &Cscale[0] /* col scale factors (output) */,
- &SL /* fact L (output)*/, &SU /* fact U (output)*/,
- NULL /* work */,
- 0 /* lwork: superlu auto allocates (input) */,
- &SB /* rhs */, &SX /* solution */,
- &recip_pivot_gross /* reciprocal pivot growth */
- /* factor max_j( norm(A_j)/norm(U_j) ). */,
- &rcond /*estimate of the reciprocal condition */
- /* number of the matrix A after equilibration */,
- &ferr[0] /* estimated forward error */,
- &berr[0] /* relative backward error */,
- &stat, &info, T());
- StatFree(&stat);
- GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
- gmm::copy(sol, X);
}
+ StatInit(&stat);
+ int info = 0;
+ R recip_pivot_gross, rcond;
+ SuperLU_gssvx(&options, &SA, &perm_c[0], &perm_r[0],
+ &etree[0] /* output */, &equed /* output */,
+ &Rscale[0] /* row scale factors (output) */,
+ &Cscale[0] /* col scale factors (output) */,
+ &SL /* fact L (output)*/, &SU /* fact U (output)*/,
+ NULL /* work */,
+ 0 /* lwork: superlu auto allocates (input) */,
+ &SB /* rhs */, &SX /* solution */,
+ &recip_pivot_gross /* reciprocal pivot growth */
+ /* factor max_j( norm(A_j)/norm(U_j) ). */,
+ &rcond /*estimate of the reciprocal condition */
+ /* number of the matrix A after equilibration */,
+ &ferr[0] /* estimated forward error */,
+ &berr[0] /* relative backward error */,
+ &stat, &info, T());
+ StatFree(&stat);
+ GMM_ASSERT1(info == 0, "SuperLU solve failed: info=" << info);
+ gmm::copy(sol, X);
+ }
template <typename T, typename V1, typename V2> inline
void mult(const SuperLU_factor<T>& P, const V1 &v1, const V2 &v2) {
@@ -404,7 +404,7 @@ namespace gmm {
}
-
+
#endif // GMM_SUPERLU_INTERFACE_H
#endif // GMM_USES_SUPERLU
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Clean up,
Konstantinos Poulios <=