From 8b0a74a327fbecb12f3d32ecff1fa28a6a6c2a15 Mon Sep 17 00:00:00 2001 From: jthorn Date: Wed, 22 Jan 2003 12:27:46 +0000 Subject: * add support for sparse-matrix Jacobians ==> works! * change default in param.ccl to use this * change default in src/include/config.h to default to no longer link in LAPACK routines * update documentation git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@931 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/dense_Jacobian.hh | 102 ++++++++++++++++++++--------------------- 1 file changed, 50 insertions(+), 52 deletions(-) (limited to 'src/elliptic/dense_Jacobian.hh') diff --git a/src/elliptic/dense_Jacobian.hh b/src/elliptic/dense_Jacobian.hh index 1e269d5..0c85f47 100644 --- a/src/elliptic/dense_Jacobian.hh +++ b/src/elliptic/dense_Jacobian.hh @@ -1,10 +1,12 @@ -// Jacobian.hh -- data structures for the Jacobian matrix +// dense_Jacobian.hh -- dense-matrix Jacobian // $Header$ // #ifdef HAVE_DENSE_JACOBIAN -// dense_Jacobian_sparsity -- (no-op) accumulate dense Jacobian sparsity info -// dense_Jacobian_matrix -- Jacobian stored as a dense matrix +// dense_Jacobian -- Jacobian stored as a dense matrix +#endif +#ifdef HAVE_DENSE_JACOBIAN__LAPACK +// dense_Jacobian__LAPACK -- dense_Jacobian with LAPACK linear solver #endif // @@ -18,44 +20,30 @@ #ifdef HAVE_DENSE_JACOBIAN // -// This class accumulates the sparsity pattern of a dense-matrix Jacobian. -// -// Since a dense-matrix Jacobian doesn't care about the sparsity pattern, -// the member functions of this class are just no-ops. +// This class stores the Jacobian as a dense matrix in Fortran (column) +// order. // -class dense_Jacobian_sparsity - : public Jacobian_sparsity +class dense_Jacobian + : public Jacobian { public: - // record the zeroing of the entire matrix - void zero_matrix() { } + // + // routines to access the matrix + // - // record the setting of a matrix element to a specified value - void set_element(int II, int JJ, fp value) { } + // get a matrix element + fp element(int II, int JJ) const + { return matrix_(JJ,II); } - // record the summing of a value into a matrix element - void sum_into_element(int II, int JJ, fp value) { } + // dense matrix ==> all elements are explicitly stored + bool is_explicitly_stored(int II, int JJ) const + { return true; } - // constructor, destructor - dense_Jacobian_sparsity(patch_system& ps, - bool print_msg_flag = false) - : Jacobian_sparsity(ps, print_msg_flag) - { } - ~dense_Jacobian_sparsity() { } - }; -#endif // HAVE_DENSE_JACOBIAN -//****************************************************************************** + // + // routines for setting values in the matrix + // -#ifdef HAVE_DENSE_JACOBIAN -// -// This class stores the Jacobian as a dense matrix in Fortran (column) -// order. -// -class dense_Jacobian_matrix - : public Jacobian_matrix - { -public: // zero the entire matrix void zero_matrix(); @@ -67,16 +55,30 @@ public: void sum_into_element(int II, int JJ, fp value) { matrix_(JJ,II) += value; } - // access (get) a matrix element - fp element(int II, int JJ) - const - { return matrix_(JJ,II); } + // + // constructor, destructor + // +protected: + dense_Jacobian(patch_system& ps, + bool print_msg_flag = false); + ~dense_Jacobian() { } + +protected: + // Fortran storage order ==> subscripts are (JJ,II) + jtutil::array2d matrix_; + }; +#endif // HAVE_DENSE_JACOBIAN - // this is a dense matrix ==> all values are stored explicitly - bool is_explicitly_stored(int II, int JJ) - const - { return true; } +//****************************************************************************** +#ifdef HAVE_DENSE_JACOBIAN__LAPACK +// +// This class defines the linear solver routine using LAPACK. +// +class dense_Jacobian__LAPACK + : public dense_Jacobian + { +public: // solve linear system J.x = rhs via LAPACK // ... rhs and x are nominal-grid gridfns // ... overwrites Jacobian matrix with its LU decomposition @@ -87,20 +89,16 @@ public: // constructor, destructor public: - dense_Jacobian_matrix(patch_system& ps, - dense_Jacobian_sparsity& JS, - bool print_msg_flag = false); - ~dense_Jacobian_matrix(); + dense_Jacobian__LAPACK(patch_system& ps, + bool print_msg_flag = false); + ~dense_Jacobian__LAPACK(); private: - // subscripts are (JJ,II) (both 0-origin) - jtutil::array2d matrix_; - // pivot vector for LAPACK routines - integer *pivot_; + integer *pivot_; // size N_rows_ // work vectors for LAPACK condition number computation - integer *iwork_; - fp *rwork_; + integer *iwork_; // size N_rows_ + fp *rwork_; // size 4*N_rows_ }; -#endif // HAVE_DENSE_JACOBIAN +#endif // HAVE_DENSE_JACOBIAN__LAPACK -- cgit v1.2.3