diff options
Diffstat (limited to 'src/elliptic/Jacobian.cc')
| -rw-r--r-- | src/elliptic/Jacobian.cc | 114 |
1 files changed, 92 insertions, 22 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; } //****************************************************************************** |