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

1 files changed, 42 insertions, 35 deletions

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
@@ -2,6 +2,8 @@
 // $Id$
 
 //
+// new_Jacobian
+//
 // Jacobian::Jacobian
 //
 // dense_Jacobian::dense_Jacobian
@@ -10,8 +12,6 @@
 // dense_Jacobian::zero_row
 // dense_Jacobian::solve_linear_system
 //
-// new_Jacobian
-//
 
 #include <stdio.h>
 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$ */
@@ -111,6 +116,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.
 //
 Jacobian::Jacobian(patch_system& ps)
@@ -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*/
-}