From 9bf8e92fc7b790a23f3f760d8aac4400f763c21c Mon Sep 17 00:00:00 2001 From: jthorn Date: Wed, 31 Jul 2002 13:12:20 +0000 Subject: fix a bug where we took A(i,j) as the infinity-norm instead of |A(i,j)| git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@676 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 77 ++++++++++++++++++++++++++---------------------- 1 file changed, 42 insertions(+), 35 deletions(-) (limited to 'src/elliptic/Jacobian.cc') diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc index 8762434..651170d 100644 --- a/src/elliptic/Jacobian.cc +++ b/src/elliptic/Jacobian.cc @@ -1,6 +1,8 @@ // Jacobian.cc -- data structures for the Jacobian matrix // $Id$ +// +// new_Jacobian // // Jacobian::Jacobian // @@ -10,8 +12,6 @@ // dense_Jacobian::zero_row // dense_Jacobian::solve_linear_system // -// new_Jacobian -// #include using std::fopen; @@ -45,7 +45,12 @@ using jtutil::error_exit; #include "../util/patch_system.hh" #include "Jacobian.hh" +// +// FIXME: Cactus's CCTK_FCALL() isn't expanded in .h files (this is a bug), +// so we include the contents of "lapack.h" inline here. +// //#include "lapack.h" +// //************************************** /* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */ /* $Id$ */ @@ -110,6 +115,28 @@ void CCTK_FCALL //****************************************************************************** //****************************************************************************** +// +// This function is an "object factory" for Jacobians: it constructs +// and returns a new-allocated Jacobian object of a specified type. +// +// FIXME: the patch system shouldn't really have to be non-const, but +// the Jacobian constructors all require this to allow the +// linear solvers to directly update gridfns +// +Jacobian& new_Jacobian(patch_system& ps, + const char Jacobian_type[]) +{ +if (STRING_EQUAL(Jacobian_type, "dense matrix")) + then return *new dense_Jacobian(ps); +else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "unknown Jacobian_type=\"%s\"!", + Jacobian_type); /*NOTREACHED*/ +} + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + // // This function constructs a Jacobian object. // @@ -156,7 +183,7 @@ void dense_Jacobian::zero_matrix() { for (int II = 0 ; II < NN() ; ++II) { - matrix_(JJ,II) = 0.0; + operator()(II,JJ) = 0.0; } } } @@ -170,7 +197,7 @@ void dense_Jacobian::zero_row(int II) { for (int JJ = 0 ; JJ < NN() ; ++JJ) { - matrix_(JJ,II) = 0.0; + operator()(II,JJ) = 0.0; } } @@ -188,23 +215,22 @@ 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 infinity_norm_flag = 0; +const integer stride = 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 +// ... max_posn = 1-origin index of A[] element with largest absolute value #if defined(FP_IS_FLOAT) - const integer posn = CCTK_FNAME(isamax)(&N2, A, &one); + const integer max_posn = CCTK_FNAME(isamax)(&N2, A, &stride); #elif defined(FP_IS_DOUBLE) - const integer posn = CCTK_FNAME(idamax)(&N2, A, &one); + const integer max_posn = CCTK_FNAME(idamax)(&N2, A, &stride); #else #error "don't know fp datatype!" #endif -const fp A_infnorm = A[posn-1]; +const fp A_infnorm = jtutil::abs(A[max_posn-1]); // LU decompose and solve the linear system // @@ -215,6 +241,7 @@ const fp A_infnorm = A[posn-1]; { x[II] = rhs[II]; } +integer info; #if defined(FP_IS_FLOAT) CCTK_FNAME(sgesv)(&N, &NRHS, A, &N, pivot_, x, &N, &info); #elif defined(FP_IS_DOUBLE) @@ -239,10 +266,12 @@ if (info > 0) // estimate infinity-norm condition number fp rcond; #if defined(FP_IS_FLOAT) - CCTK_FNAME(sgecon_wrapper)(&zero, &N, A, &N, &A_infnorm, &rcond, + CCTK_FNAME(sgecon_wrapper)(&infinity_norm_flag, + &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, + CCTK_FNAME(dgecon_wrapper)(&infinity_norm_flag, + &N, A, &N, &A_infnorm, &rcond, rwork_, iwork_, &info); #else #error "don't know fp datatype!" @@ -252,25 +281,3 @@ if (rcond == 0.0) // *** (singular matrix) return 1.0/rcond; } - -//****************************************************************************** -//****************************************************************************** -//****************************************************************************** - -// -// This function is an "object factory" for Jacobians: it constructs -// and returns a new-allocated Jacobian object of a specified type. -// -// FIXME: the patch system shouldn't really have to be non-const, but -// the Jacobian constructors all require this to allow the -// linear solvers to directly update gridfns -// -Jacobian& new_Jacobian(patch_system& ps, - const char Jacobian_type[]) -{ -if (STRING_EQUAL(Jacobian_type, "dense matrix")) - then return *new dense_Jacobian(ps); -else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, - "unknown Jacobian_type=\"%s\"!", - Jacobian_type); /*NOTREACHED*/ -} -- cgit v1.2.3