diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/elliptic/Jacobian.cc | 114 | ||||
-rw-r--r-- | src/elliptic/Jacobian.hh | 21 | ||||
-rw-r--r-- | src/elliptic/gecon_wrapper.F77 | 53 | ||||
-rw-r--r-- | src/elliptic/lapack.h | 32 | ||||
-rw-r--r-- | src/elliptic/make.code.defn | 2 |
5 files changed, 192 insertions, 30 deletions
diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 603abac..8762434 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -6,7 +6,7 @@ // // dense_Jacobian::dense_Jacobian // dense_Jacobian::~dense_Jacobian -// dense_Jacobian::zero +// dense_Jacobian::zero_matrix // dense_Jacobian::zero_row // dense_Jacobian::solve_linear_system // @@ -47,7 +47,7 @@ using jtutil::error_exit; #include "Jacobian.hh" //#include "lapack.h" //************************************** -/* lapack.h -- C/C++ prototypes for (some) LAPACK routines */ +/* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */ /* $Id$ */ /* @@ -60,6 +60,17 @@ using jtutil::error_exit; extern "C" { #endif +/* + * ***** BLAS ***** + */ +integer CCTK_FCALL + CCTK_FNAME(isamax)(const integer* N, const float SX[], const integer* incx); +integer CCTK_FCALL + CCTK_FNAME(idamax)(const integer* N, const double DX[], const integer* incx); + +/* + * ***** LAPACK ***** + */ void CCTK_FCALL CCTK_FNAME(sgesv)(const integer* N, const integer* NRHS, float A[], const int* LDA, @@ -71,6 +82,25 @@ void CCTK_FCALL integer IPIV[], double B[], const integer* LDB, integer* info); +/* + * ***** wrappers (for passing character-string args) ***** + */ +/* norm_int = 0 for infinity-norm, 1 for 1-norm */ +void CCTK_FCALL + CCTK_FNAME(sgecon_wrapper)(const integer* norm_int, + const integer* N, + float A[], const integer* LDA, + const float* anorm, float* rcond, + float WORK[], integer IWORK[], + integer* info); +void CCTK_FCALL + CCTK_FNAME(dgecon_wrapper)(const integer* norm_int, + const integer* N, + double A[], const integer* LDA, + const double* anorm, double* rcond, + double WORK[], integer IWORK[], + integer* info); + #ifdef __cplusplus }; /* extern "C" */ #endif @@ -98,7 +128,9 @@ Jacobian::Jacobian(patch_system& ps) dense_Jacobian::dense_Jacobian(patch_system& ps) : Jacobian(ps), matrix_(0,NN()-1, 0,NN()-1), - pivot_(new integer[NN()]) + pivot_(new integer[NN()]), + iwork_(new integer[NN()]), + rwork_(new fp[4*NN()]) // no comma { } //****************************************************************************** @@ -108,6 +140,8 @@ dense_Jacobian::dense_Jacobian(patch_system& ps) // dense_Jacobian::~dense_Jacobian() { +delete[] iwork_; +delete[] rwork_; delete[] pivot_; } @@ -116,7 +150,7 @@ delete[] pivot_; // // This function zeros a dense_Jacobian object. // -void dense_Jacobian::zero() +void dense_Jacobian::zero_matrix() { for (int JJ = 0 ; JJ < NN() ; ++JJ) { @@ -144,43 +178,79 @@ void dense_Jacobian::zero_row(int II) // // This function solves the linear system J.x = rhs, with rhs and x -// being nominal-grid gridfns. The computation is done using the LAPACK -// [sd]gesv() subroutines; which modify the Jacobian matrix J in-place -// for the LU decomposition. +// being nominal-grid gridfns, using LU decomposition. It returns the +// estimated infinity-norm condition number of the linear system, or +// 0.0 if the matrix is numerically singular // -void dense_Jacobian::solve_linear_system(int rhs_gfn, int x_gfn) +fp dense_Jacobian::solve_linear_system(int rhs_gfn, int x_gfn) { const fp *rhs = my_patch_system().gridfn_data(rhs_gfn); fp *x = my_patch_system().gridfn_data(x_gfn); +fp *A = matrix_.data_array(); +const integer zero = 0; +const integer one = 1; +const integer NRHS = 1; +const integer N = NN(); +const integer N2 = NN() * NN(); +integer info; + +// compute the infinity-norm of the matrix A +// ... posn = 1-origin index of A[] element with max absolute value +#if defined(FP_IS_FLOAT) + const integer posn = CCTK_FNAME(isamax)(&N2, A, &one); +#elif defined(FP_IS_DOUBLE) + const integer posn = CCTK_FNAME(idamax)(&N2, A, &one); +#else + #error "don't know fp datatype!" +#endif +const fp A_infnorm = A[posn-1]; + +// LU decompose and solve the linear system // -// [sd]gesv() use an "in out" design, where the same argument is used for -// both rhs and x ==> first copy rhs to x so we can pass that to [sd]gesv() +// ... [sd]gesv() use an "in out" design, where the same argument +// is used for both rhs and x ==> we must first copy rhs to x // for (int II = 0 ; II < NN() ; ++II) { x[II] = rhs[II]; } - -integer N = NN(); -integer NRHS = 1; -integer info; - -#ifdef FP_IS_FLOAT - CCTK_FNAME(sgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); -#endif -#ifdef FP_IS_DOUBLE - CCTK_FNAME(dgesv)(&N, &NRHS, matrix_.data_array(), &N, pivot_, x, &N, &info); +#if defined(FP_IS_FLOAT) + CCTK_FNAME(sgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info); +#elif defined(FP_IS_DOUBLE) + CCTK_FNAME(dgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info); +#else + #error "don't know fp datatype!" #endif -if (info != 0) +if (info < 0) then error_exit(ERROR_EXIT, "\n" "***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n" -" error return info=%d from [sd]gesv() LAPACK routine!\n" +" error return (bad argument) info=%d from [sd]gesv() LAPACK routine!\n" , rhs_gfn, x_gfn, int(info)); /*NOTREACHED*/ + +if (info > 0) + then return 0.0; // *** ERROR RETURN *** + // *** (singular matrix) + +// estimate infinity-norm condition number +fp rcond; +#if defined(FP_IS_FLOAT) + CCTK_FNAME(sgecon_wrapper)(&zero, &N, A, &N, &A_infnorm, &rcond, + rwork_, iwork_, &info); +#elif defined(FP_IS_DOUBLE) + CCTK_FNAME(dgecon_wrapper)(&zero, &N, A, &N, &A_infnorm, &rcond, + rwork_, iwork_, &info); +#else + #error "don't know fp datatype!" +#endif +if (rcond == 0.0) + then return 0.0; // *** ERROR RETURN *** + // *** (singular matrix) +return 1.0/rcond; } //****************************************************************************** diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh index 7bd37db..d5c2e1d 100644 --- a/src/elliptic/Jacobian.hh +++ b/src/elliptic/Jacobian.hh @@ -70,13 +70,16 @@ public: virtual bool is_explicitly_stored(int II, int JJ) const = 0; // zero the whole matrix or a single row - virtual void zero() = 0; + virtual void zero_matrix() = 0; virtual void zero_row(int II) = 0; // solve linear system J.x = rhs // ... rhs and x are nominal-grid gridfns // ... may modify Jacobian matrix (eg for LU decomposition) - virtual void solve_linear_system(int rhs_gfn, int x_gfn) = 0; + // ... returns 0.0 if matrix is numerically singular + // condition number (> 0.0) if known + // -1.0 if condition number is unknown + virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0; // basic meta-info patch_system& my_patch_system() const { return ps_; } @@ -113,13 +116,15 @@ public: { return true; } // zero the whole matrix or a single row - void zero(); + void zero_matrix(); void zero_row(int II); - // solve linear system J.x = rhs + // solve linear system J.x = rhs via LAPACK // ... rhs and x are nominal-grid gridfns - // ... may modify Jacobian matrix (eg for LU decomposition) - void solve_linear_system(int rhs_gfn, int x_gfn); + // ... overwrites Jacobian matrix with its LU decomposition + // ... returns 0.0 if matrix is numerically singular, or + // condition number (> 0.0) + fp solve_linear_system(int rhs_gfn, int x_gfn); // constructor, destructor public: @@ -132,6 +137,10 @@ private: // pivot vector for LAPACK routines integer *pivot_; + + // work vectors for LAPACK condition number computation + integer *iwork_; + fp *rwork_; }; //****************************************************************************** diff --git a/src/elliptic/gecon_wrapper.F77 b/src/elliptic/gecon_wrapper.F77 new file mode 100644 index 0000000..6feb1e8 --- /dev/null +++ b/src/elliptic/gecon_wrapper.F77 @@ -0,0 +1,53 @@ +c gecon_wrapper.f -- wrapper routines for LAPACK [sd]gecon() +c $Id$ + +c +c These subroutines are wrappers around the LAPACK [sd]gecon() subroutines. +c These subroutines take only integer/real/double precision arguments, +c avoiding problems with C/C++ --> Fortran passing of the character string +c arguments used by [sd]gecon(). +c +c Arguments: +c norm_int = (in) 0 ==> infinity-norm +c 1 ==> 1-norm +c + + subroutine sgecon_wrapper(norm_int, + $ N, A, LDA, anorm, rcond, + $ WORK, IWORK, info) + integer norm_int + integer N, LDA + real A(LDA,N) + real anorm, rcond + real WORK(*) + integer iwork(*) + integer info + if (norm_int .eq. 0) then + call sgecon('I', N,A,LDA, anorm,rcond, WORK,IWORK, info) + else if (norm_int .eq. 1) then + call sgecon('1', N,A,LDA, anorm,rcond, WORK,IWORK, info) + else + info = -1; + end if + return + end + + subroutine dgecon_wrapper(norm_int, + $ N, A, LDA, anorm, rcond, + $ WORK, IWORK, info) + integer norm_int + integer N, LDA + double precision A(LDA,N) + double precision anorm, rcond + double precision WORK(*) + integer iwork(*) + integer info + if (norm_int .eq. 0) then + call dgecon('I', N,A,LDA, anorm,rcond, WORK,IWORK, info) + else if (norm_int .eq. 1) then + call dgecon('1', N,A,LDA, anorm,rcond, WORK,IWORK, info) + else + info = -1; + end if + return + end diff --git a/src/elliptic/lapack.h b/src/elliptic/lapack.h index 205bc96..f0582c9 100644 --- a/src/elliptic/lapack.h +++ b/src/elliptic/lapack.h @@ -1,4 +1,4 @@ -/* lapack.h -- C/C++ prototypes for (some) LAPACK routines */ +/* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */ /* $Id$ */ /* @@ -11,6 +11,17 @@ extern "C" { #endif +/* + * ***** BLAS ***** + */ +integer CCTK_FCALL + CCTK_FNAME(isamax)(const integer* N, const float SX[], const integer* incx); +integer CCTK_FCALL + CCTK_FNAME(idamax)(const integer* N, const double DX[], const integer* incx); + +/* + * ***** LAPACK ***** + */ void CCTK_FCALL CCTK_FNAME(sgesv)(const integer* N, const integer* NRHS, float A[], const int* LDA, @@ -22,6 +33,25 @@ void CCTK_FCALL integer IPIV[], double B[], const integer* LDB, integer* info); +/* + * ***** wrappers (for passing character-string args) ***** + */ +/* norm_int = 0 for infinity-norm, 1 for 1-norm */ +void CCTK_FCALL + CCTK_FNAME(sgecon_wrapper)(const integer* norm_int, + const integer* N, + float A[], const integer* LDA, + const float* anorm, float* rcond, + float WORK[], integer IWORK[], + integer* info); +void CCTK_FCALL + CCTK_FNAME(dgecon_wrapper)(const integer* norm_int, + const integer* N, + double A[], const integer* LDA, + const double* anorm, double* rcond, + double WORK[], integer IWORK[], + integer* info); + #ifdef __cplusplus }; /* extern "C" */ #endif diff --git a/src/elliptic/make.code.defn b/src/elliptic/make.code.defn index 59e1bd7..3bc82f4 100644 --- a/src/elliptic/make.code.defn +++ b/src/elliptic/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = Jacobian.cc +SRCS = Jacobian.cc gecon_wrapper.F77 # Subdirectories containing source files SUBDIRS = |