1 files changed, 43 insertions, 11 deletions

diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 2c64171..d948f07 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -4,6 +4,9 @@
 //
 // Jacobian -- abstract base class to describe a Jacobian matrix
 //
+// print_Jacobian - print a Jacobian
+// print_Jacobians - print two Jacobians
+//
 // dense_Jacobian - Jacobian stored as a dense matrix
 //
 
@@ -47,10 +50,6 @@
 class	Jacobian
 	{
 public:
-	// basic meta-info
-	patch_system& my_patch_system();
-
-public:
 	// access a matrix element via row/column indices
 	virtual const fp& operator()(int II, int JJ) const = 0;	// rvalue
 	virtual       fp& operator()(int II, int JJ)       = 0;	// lvalue
@@ -72,22 +71,44 @@ public:
 		return operator()(II,JJ);
 		}
 
+	// is a given element explicitly stored, or implicitly 0 via sparsity
+	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_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;
+
 	// basic meta-info
-	const patch_system& my_patch_system() const { return ps_; }
+	patch_system& my_patch_system() const { return ps_; }
+	int NN() const { return NN_; }
 
 	// constructor, destructor
-	Jacobian(const patch_system& ps);
+	// ... for analyze-factor-solve sparse-matrix derived classes,
+	//     construction may include the "analyze" operation
+	Jacobian(patch_system& ps);
 	virtual ~Jacobian() { }
 
-private:
-	const patch_system& ps_;
+protected:
+	patch_system& ps_;
+	int NN_;	// ps_.N_grid_points()
 	};
 
 //******************************************************************************
+
+//
+// These functions print 1 or 2 Jacobians to a named output file;
+// in the latter case they also print the error in the SD Jacobian.
+//
+void print_Jacobian(const char file_name[], const Jacobian& Jac);
+void print_Jacobians(const char file_name[],
+		     const Jacobian& SD_Jac, const Jacobian& NP_Jac);
+
+//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -105,16 +126,27 @@ public:
 	fp& operator()(int II, int JJ)				// lvalue
 		{ return matrix_(JJ,II); }
 
+	bool is_explicitly_stored(int II, int JJ) const
+		{ return true; }
+
 	// zero the whole matrix or a single row
 	void zero();
 	void zero_row(int II);
 
+	// solve linear system J.x = rhs
+	// ... 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);
+
 	// constructor, destructor
 public:
-	dense_Jacobian(const patch_system& ps);
-	~dense_Jacobian() { }
+	dense_Jacobian(patch_system& ps);
+	~dense_Jacobian();
 
 private:
-	// subscripts are (JJ,II)
+	// subscripts are (JJ,II) (both 0-origin)
 	jtutil::array2d<fp> matrix_;
+
+	// pivot vector for LAPACK routines
+	integer *pivot_;
 	};