add support for estimating condition number when we solve a linear system

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@668 f88db872-0e4f-0410-b76b-b9085cfa78c5

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 =