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 <logar...@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