* solve_linear_system() now returns *reciprocal* condition number

* add hooks for UMFPACK linear solver


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

2 files changed, 30 insertions, 4 deletions

diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc
index 08d8ecc..faaf767 100644
--- a/src/elliptic/Jacobian.cc
+++ b/src/elliptic/Jacobian.cc
@@ -66,6 +66,14 @@ else if (STRING_EQUAL(Jacobian_store_solve_method_string,
   #endif
 	}
 
+else if (STRING_EQUAL(Jacobian_store_solve_method_string,
+		      "row-oriented sparse matrix/UMFPACK"))
+   then {
+  #ifdef HAVE_ROW_SPARSE_JACOBIAN__UMFPACK
+	return Jacobian__row_sparse_matrix__UMFPACK;
+  #endif
+	}
+
 else	error_exit(ERROR_EXIT,
 "decode_Jacobian_store_solve_method():\n"
 "        unknown Jacobian_store_solve_method_string=\"%s\"!\n",
@@ -110,6 +118,11 @@ switch	(Jac_method)
 	return new row_sparse_Jacobian__ILUCG(ps, print_msg_flag);
 #endif
 
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__UMFPACK
+  case Jacobian__row_sparse_matrix__UMFPACK:
+	return new row_sparse_Jacobian__UMFPACK(ps, print_msg_flag);
+#endif
+
   default:
 	error_exit(ERROR_EXIT,
 		   "new_Jacobian(): unknown method=(int)%d!\n",
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 327cabf..803b1c7 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -30,6 +30,7 @@
 //		dense_Jacobian__LAPACK
 //	    row_sparse_Jacobian
 //		row_sparse_Jacobian__ILUCG
+//		row_sparse_Jacobian__UMFPACK
 // each derived class is inside a corresponding #ifdef (set or unset
 // as appropriate, in "../include/config.h").
 //
@@ -62,7 +63,9 @@
 // point number within the patch system.
 //
 // Since many of the derived classes use Fortran routines, we also use
-// 1-origin indices; these have a leading "f", eg fII/fJJ.
+// 1-origin indices; these have a leading "f", eg fII/fJJ.  Finally, we
+// use generic indices that might be either 0-origin or 1-origin; these
+// have a leading "g", eg gII/gJJ.
 //
 
 //
@@ -94,11 +97,13 @@ public:
 		{ return ps_.patch_irho_isigma_of_gpn(II, irho,isigma); }
 
 
+#ifdef NOT_USED
 	//
 	// convert C <--> Fortran indices
 	//
 	int csub(int f) const { return f-1; }
 	int fsub(int c) const { return c+1; }
+#endif
 
 
 	//
@@ -130,9 +135,11 @@ public:
 	// solve linear system J.x = rhs
 	// ... rhs and x are nominal-grid gridfns
 	// ... may modify Jacobian matrix (eg for LU decomposition)
-	// ... returns 0.0 if matrix is numerically singular
-	//             condition number (> 0.0) if known
+	// ... returns estimated reciprocal condition number if known,
 	//	       -1.0 if condition number is unknown
+	//	       ... rcond = 0.0         ==> exactly singular
+	//		   rcond < DBL_EPSILON ==> numerically singular
+	//		   rcond = 1.0         ==> orthogonal matrix
 	// ... once this has been called, the sparsity pattern should
 	//     not be changed, i.e. no new nonzeros should be introduced
 	//     into the matrix
@@ -182,6 +189,9 @@ struct	linear_solver_pars
 		fp  error_tolerance;
 		int max_CG_iterations;
 		} ILUCG_pars;
+	struct	UMFPACK_pars
+		{
+		} UMFPACK_pars;
 	};
 
 //******************************************************************************
@@ -192,7 +202,10 @@ enum	Jacobian_store_solve_method
 	Jacobian__dense_matrix__LAPACK,
 #endif
 #ifdef HAVE_ROW_SPARSE_JACOBIAN__ILUCG
-	Jacobian__row_sparse_matrix__ILUCG // no comma on last entry in enum
+	Jacobian__row_sparse_matrix__ILUCG,
+#endif
+#ifdef HAVE_ROW_SPARSE_JACOBIAN__UMFPACK
+	Jacobian__row_sparse_matrix__UMFPACK // no comma on last entry in enum
 #endif
 	};