start adding support for sparse matrix Jacobians (doesn't work yet)

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

3 files changed, 195 insertions, 337 deletions

diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc
index 4ac18c3..048c6ac 100644
--- a/src/elliptic/Jacobian.cc
+++ b/src/elliptic/Jacobian.cc
@@ -1,21 +1,12 @@
-// Jacobian.cc -- data structures for the Jacobian matrix
+// Jacobian.cc -- generic routines for the Jacobian matrix
 // $Header$
 
 //
 // <<<literal contents of "lapack.hh">>>
 //
 // decode_Jacobian_type
-// new_Jacobian
-//
-// Jacobian::Jacobian
-//
-#ifdef HAVE_DENSE_JACOBIAN
-// dense_Jacobian::dense_Jacobian
-// dense_Jacobian::~dense_Jacobian
-// dense_Jacobian::zero_matrix
-// dense_Jacobian::zero_row
-// dense_Jacobian::solve_linear_system
-#endif
+// new_Jacobian_sparsity
+// new_Jacobian_matrix
 //
 
 #include <stdio.h>
@@ -45,76 +36,8 @@ using jtutil::error_exit;
 #include "../patch/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"
-//
-
-//***** begin "lapack.h" contents ******
-/* lapack.h -- C/C++ prototypes for (some) BLAS+LAPACK+wrapper routines */
-/* $Header$ */
-
-/*
- * prerequisites:
- *	"cctk.h"
- *	"config.hh"
- */
-
-#ifdef __cplusplus
-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,
-		    integer IPIV[],
-		    float B[], const integer* LDB, integer* info);
-void CCTK_FCALL
-  CCTK_FNAME(dgesv)(const integer* N, const integer* NRHS,
-		    double A[], const int* LDA,
-		    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
-//***** end "lapack.h" contents ********
+#include "dense_Jacobian.hh"
+#include "row_sparse_Jacobian.hh"
 
 //******************************************************************************
 //******************************************************************************
@@ -131,200 +54,101 @@ if	(STRING_EQUAL(Jacobian_type_string, "dense matrix"))
    then {
   #ifdef HAVE_DENSE_JACOBIAN
 	return Jacobian_type__dense_matrix;
-  #else
-	error_exit(ERROR_EXIT,
-"\n"
-"   decode_Jacobian_type():\n"
-"        Jacobian type \"dense matrix\" is not configured in this binary;\n"
-"        see \"src/include/config.hh\" for details on what Jacobian types\n"
-"        are configured and how to change this\n");		/*NOTREACHED*/
   #endif
 	}
+
+else if (STRING_EQUAL(Jacobian_type_string, "row-oriented sparse matrix"))
+   then {
+  #ifdef HAVE_ROW_SPARSE_JACOBIAN
+	return Jacobian_type__row_sparse_matrix;
+  #endif
+	}
+
 else	error_exit(ERROR_EXIT,
 "decode_Jacobian_type(): unknown Jacobian_type_string=\"%s\"!\n",
 		   Jacobian_type_string);			/*NOTREACHED*/
+
+// fall through to here ==> we recognize the matrix type,
+//                          but it's not configured
+error_exit(ERROR_EXIT,
+"\n"
+"   decode_Jacobian_type():\n"
+"        Jacobian type \"%s\" is not configured in this binary;\n"
+"        see \"src/include/config.hh\" for details on what Jacobian types\n"
+"        are configured and how to change this\n"
+	   ,
+	   Jacobian_type_string);				/*NOTREACHED*/
 }
 
 //******************************************************************************
 
 //
-// This function is an "object factory" for Jacobians: it constructs
-// and returns a new-allocated Jacobian object of a specified type.
+// This function is an "object factory" for Jacobian_sparsity objects:
+// it constructs and returns a pointer to a new-allocated Jacobian_sparsity
+// object of the specified derived 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, enum Jacobian_type Jac_type)
+Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type,
+					 patch_system& ps,
+					 bool print_msg_flag /* = false */)
 {
 switch	(Jac_type)
 	{
 #ifdef HAVE_DENSE_JACOBIAN
   case Jacobian_type__dense_matrix:
-	return *new dense_Jacobian(ps);
+	return new dense_Jacobian_sparsity(ps, print_msg_flag);
+#endif
+
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
+  case Jacobian_type__row_sparse_matrix:
+	return new row_sparse_Jacobian_sparsity(ps, print_msg_flag);
 #endif
+
   default:
 	error_exit(ERROR_EXIT,
-"new_Jacobian(): unknown Jacobian_type=(int)%d!\n",
+"new_Jacobian_sparsity(): unknown Jacobian_type=(int)%d!\n",
 		   Jac_type);					/*NOTREACHED*/
 	}
 }
 
 //******************************************************************************
-//******************************************************************************
-//******************************************************************************
 
 //
-// This function constructs a  Jacobian  object.
+// This function is an "object factory" for Jacobian_matrix objects:
+// it constructs and returns a pointer to a new-allocated Jacobian_matrix
+// object of the specified derived type.
 //
-Jacobian::Jacobian(patch_system& ps)
-	: ps_(ps),
-	  NN_(ps.N_grid_points())
-{ }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// This function constructs a  dense_Jacobian  object.
-//
-dense_Jacobian::dense_Jacobian(patch_system& ps)
-	: Jacobian(ps),
-	  matrix_(0,NN()-1, 0,NN()-1),
-	  pivot_(new integer[NN()]),
-	  iwork_(new integer[NN()]),
-	  rwork_(new fp[4*NN()]) // no comma
-{ }
-#endif
-
-//******************************************************************************
-
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// THis function destroys a  dense_Jacobian  object.
-//
-dense_Jacobian::~dense_Jacobian()
-{
-delete[] iwork_;
-delete[] rwork_;
-delete[] pivot_;
-}
-#endif
-
-//******************************************************************************
-
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// This function zeros a  dense_Jacobian  object.
+// 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
 //
-void dense_Jacobian::zero_matrix()
+Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type,
+				     patch_system& ps,
+				     Jacobian_sparsity& JS,
+				     bool print_msg_flag /* = false */)
 {
-	for (int JJ = 0 ; JJ < NN() ; ++JJ)
+switch	(Jac_type)
 	{
-		for (int II = 0 ; II < NN() ; ++II)
-		{
-		operator()(II,JJ) = 0.0;
-		}
-	}
-}
-#endif
-
-//******************************************************************************
-
 #ifdef HAVE_DENSE_JACOBIAN
-//
-// This function zeros a single row of a  dense_Jacobian  object.
-//
-void dense_Jacobian::zero_row(int II)
-{
-	for (int JJ = 0 ; JJ < NN() ; ++JJ)
-	{
-	operator()(II,JJ) = 0.0;
-	}
-}
+  case Jacobian_type__dense_matrix:
+	return new dense_Jacobian_matrix
+			(ps, static_cast<dense_Jacobian_sparsity&>(JS),
+			 print_msg_flag);
 #endif
 
-//******************************************************************************
-
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// This function solves the linear system J.x = rhs, with rhs and x
-// being nominal-grid gridfns, using LAPACK LU-decomposition routines.
-// It returns the estimated infinity-norm condition number of the linear
-// system, or 0.0 if the matrix is numerically singular
-//
-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 infinity_norm_flag = 0;
-const integer stride = 1;
-const integer NRHS = 1;
-const integer N = NN();
-const integer N2 = NN() * NN();
-
-// compute the infinity-norm of the matrix A
-// ... max_posn = 1-origin index of A[] element with largest absolute value
-#if   defined(FP_IS_FLOAT)
-  const integer max_posn = CCTK_FNAME(isamax)(&N2, A, &stride);
-#elif defined(FP_IS_DOUBLE)
-  const integer max_posn = CCTK_FNAME(idamax)(&N2, A, &stride);
-#else
-  #error "don't know fp datatype!"
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
+  case Jacobian_type__row_sparse_matrix:
+	return new row_sparse_Jacobian_matrix
+			(ps, static_cast<row_sparse_Jacobian_sparsity&>(JS),
+			 print_msg_flag);
 #endif
-const fp A_infnorm = jtutil::abs(A[max_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 ==> we must first copy rhs to x
-//
-	for (int II = 0 ; II < NN() ; ++II)
-	{
-	x[II] = rhs[II];
+  default:
+	error_exit(ERROR_EXIT,
+"new_Jacobian_matrix(): unknown Jacobian_type=(int)%d!\n",
+		   Jac_type);					/*NOTREACHED*/
 	}
-integer 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)
-   then error_exit(ERROR_EXIT,
-"\n"
-"***** dense_Jacobian::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\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)(&infinity_norm_flag,
-			     &N, A, &N, &A_infnorm, &rcond,
-			     rwork_, iwork_, &info);
-#elif defined(FP_IS_DOUBLE)
-  CCTK_FNAME(dgecon_wrapper)(&infinity_norm_flag,
-			     &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;
 }
-#endif
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 3faf2d1..d4b90ab 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -1,11 +1,10 @@
-// Jacobian.hh -- data structures for the Jacobian matrix
+// Jacobian.hh -- generic data structures for the Jacobian matrix
 // $Header$
 
 //
-// Jacobian -- abstract base class to describe a Jacobian matrix
-#ifdef HAVE_DENSE_JACOBIAN
-// dense_Jacobian -- Jacobian stored as a dense matrix
-#endif
+// Jacobian -- ABC to describe Jacobian matrix
+// Jacobian_sparsity -- ABC to accumulate Jacobian sparsity info
+// Jacobian_matrix -- ABC to store/use Jacobian matrix
 //
 // decode_Jacobian_type - decode string into internal enum
 // new_Jacobian - factory method
@@ -13,29 +12,53 @@
 
 //
 // prerequisites:
-//	"stdc.h"
-//	"config.hh"
-//	"../jtutil/util.hh"
-//	"../jtutil/array.hh"
-//	"../jtutil/cpm_map.hh"
-//	"../jtutil/linear_map.hh"
-//	"coords.hh"
-//	"grid.hh"
-//	"fd_grid.hh"
-//	"patch.hh"
-//	"patch_edge.hh"
-//	"patch_interp.hh"
-//	"ghost_zone.hh"
-//	"patch_system.hh"
+//	"../patch/patch_system.hh"
 //
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
-// A Jacobian object stores the (a) Jacobian matrix for a patch system.
-// Jacobian is an abstract base class, from which we derive a separate
-// concrete class for each type of sparsity pattern (more precisely for
-// each type of storage scheme) of the Jacobian matrix.
+// These classes are used to store and manipulate Jacobian matrices
+// for a patch system.
+//
+// The inheritance diagram for the Jacobian classes looks like this:
+//	Jacobian
+//	    Jacobian_sparsity
+//		dense_Jacobian_sparsity		#ifdef HAVE_DENSE_JACOBIAN
+//		row_sparse_Jacobian_sparsity	#ifdef HAVE_ROW_SPARSE_JACOBIAN
+//	    Jacobian_matrix
+//		dense_Jacobian_matrix		#ifdef HAVE_DENSE_JACOBIAN
+//		row_sparse_Jacobian_matrix	#ifdef HAVE_ROW_SPARSE_JACOBIAN
+// where the most-derived classes are all conditional on the noted #ifdef
+// symbols.
+//
+// The  Jacobian_sparsity  objects are used to accumulate information
+// on the Jacobian's sparsity pattern; this is then used in constructing
+// the corresponding  Jacobian_matrix  object (which actually stores
+// the Jacobian, and contains member functions to solve elliptic systems
+// using the Jacobian as the LHS matrix).
+//
+// We assume a "traversal" interface for computing Jacobians, defined
+// by the  Jacobian  API.  Because this API is shared by (all) the
+// Jacobian_sparsity and Jacobian_matrix classes, the same traversal
+// routines (prototyped to update  Jacobian  objects) can be used both
+// to record sparsity patterns (ignoring the actual floating-point values)
+// and to store the Jacobian matrix.
+//
+// Thus the typical code to construct and use a Jacobian matrix looks
+// like this:
+//	// prototype
+//	void trverse_Jacobian(const patch_system& ps, Jacobian& Jac);
+//
+//	Jacobian_sparsity& JS = new_Jacobian_sparsity(Jac_type, ps);
+//	traverse_Jacobian(ps, JS);
+//	Jacobian_matrix& Jac = new_Jacobian_matrix(Jac_type, ps, JS);
+//	traverse_Jacobian(ps, Jac);
+//	Jac.solve_linear_system(...);
+//
+
 //
 // A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid
 // point number within the patch system.
@@ -46,55 +69,49 @@
 // very general, but matches our present use in this thorn.
 //
 
+//******************************************************************************
+
+// ABC to define Jacobian-traversal interface
 class	Jacobian
 	{
 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
+	// basic meta-info
+	patch_system& my_patch_system() const { return ps_; }
+	int NN() const { return NN_; }
 
-	// access a matrix element via row index and column (patch,irho,isigma)
-	const fp& operator()					// rvalue
-		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
+	// convert (patch,irho,isigma) <--> row/column index
+	int II_of_patch_irho_isigma(const patch& p, int irho, int isigma)
 		const
-		{
-		const int JJ
-		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
-		return operator()(II,JJ);
-		}
-	fp& operator()						// lvalue
-		(int II,   const patch& JJ_patch, int JJ_irho, int JJ_isigma)
-		{
-		const int JJ
-		   = ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
-		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;
+		{ return ps_.gpn_of_patch_irho_isigma(p, irho,isigma); }
+	const patch& patch_irho_isigma_of_II(int II, int& irho, int& isigma)
+		const
+		{ return ps_.patch_irho_isigma_of_gpn(II, irho,isigma); }
 
-	// zero the whole matrix or a single row
+	// zero the entire matrix
 	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)
-	// ... 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;
+	// set a matrix element to a specified value
+	virtual void set_element(int II, int JJ, fp value) = 0;
 
-	// basic meta-info
-	patch_system& my_patch_system() const { return ps_; }
-	int NN() const { return NN_; }
+	// sum a value into a matrix element
+	virtual void sum_into_element(int II, int JJ, fp value) = 0;
 
+protected:
 	// constructor, destructor
-	// ... for analyze-factor-solve sparse-matrix derived classes,
-	//     construction may include the "analyze" operation
-	Jacobian(patch_system& ps);
+	Jacobian(patch_system& ps)
+		: ps_(ps),
+		  NN_(ps.N_grid_points())
+		{ }
+public:
 	virtual ~Jacobian() { }
 
+private:
+	// we forbid copying and passing by value
+	// by declaring the copy constructor and assignment operator
+	// private, but never defining them
+	Jacobian(const Jacobian &rhs);
+	Jacobian& operator=(const Jacobian &rhs);
+
 protected:
 	patch_system& ps_;
 	int NN_;	// ps_.N_grid_points()
@@ -102,65 +119,80 @@ protected:
 
 //******************************************************************************
 
-#ifdef HAVE_DENSE_JACOBIAN
-//
-// This class stores the Jacobian as a dense matrix in Fortran (column)
-// order.
-//
-class	dense_Jacobian
+// ABC to accumulate Jacobian sparsity info
+class	Jacobian_sparsity
 	: public Jacobian
 	{
 public:
-	// access a matrix element via row/column indices
-	const fp& operator()(int II, int JJ) const		// rvalue
-		{ return matrix_(JJ,II); }
-	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_matrix();
-	void zero_row(int II);
+	Jacobian_sparsity(patch_system& ps,
+			  bool print_msg_flag = false)
+		: Jacobian(ps)
+		{ }
+	virtual ~Jacobian_sparsity() { }
+	};
 
-	// solve linear system J.x = rhs via LAPACK
-	// ... rhs and x are nominal-grid gridfns
-	// ... 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
+// Jacobian_matrix -- ABC to store/use Jacobian matrix
+class	Jacobian_matrix
+	: public Jacobian
+	{
 public:
-	dense_Jacobian(patch_system& ps);
-	~dense_Jacobian();
+	// access (get) a matrix element
+	virtual fp element(int II, int JJ) const = 0;
 
-private:
-	// subscripts are (JJ,II) (both 0-origin)
-	jtutil::array2d<fp> matrix_;
+	// is a given element explicitly stored, or implicitly 0 via sparsity
+	virtual bool is_explicitly_stored(int II, int JJ) const = 0;
 
-	// pivot vector for LAPACK routines
-	integer *pivot_;
+	// 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
+	//	       -1.0 if condition number is unknown
+	virtual fp solve_linear_system(int rhs_gfn, int x_gfn) = 0;
 
-	// work vectors for LAPACK condition number computation
-	integer *iwork_;
-	fp *rwork_;
+	// constructor, destructor
+	// ... note Jacobian_sparsity argument of constructor is non-const
+	//     to allow us to take over ownership of data structures
+	Jacobian_matrix(patch_system& ps,
+			Jacobian_sparsity& JS,
+			bool print_msg_flag = false)
+		: Jacobian(ps)
+		{ }
+	virtual ~Jacobian_matrix() { }
 	};
-#endif	// HAVE_DENSE_JACOBIAN
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 enum	Jacobian_type
 	{
 #ifdef HAVE_DENSE_JACOBIAN
-	Jacobian_type__dense_matrix	// no comma on last entry in enum
+	Jacobian_type__dense_matrix,
+#endif
+#ifdef HAVE_ROW_SPARSE_JACOBIAN
+	Jacobian_type__row_sparse_matrix // no comma on last entry in enum
 #endif
 	};
 
+//
+// prototypes of various Jacobian-related functions
+//
+
 // decode string into internal enum
 enum Jacobian_type
   decode_Jacobian_type(const char Jacobian_type_string[]);
 
-// construct and return new-allocated Jacobian object of specified type
-Jacobian& new_Jacobian(patch_system& ps, enum Jacobian_type Jac_type);
+// "object factory" routines to construct and return
+// pointers to new-allocated objects of specified derived types
+class Jacobian_sparsity;
+class Jacobian_matrix;
+Jacobian_sparsity* new_Jacobian_sparsity(enum Jacobian_type Jac_type,
+					 patch_system& ps,
+					 bool print_msg_flag = false);
+Jacobian_matrix* new_Jacobian_matrix(enum Jacobian_type Jac_type,
+				     patch_system& ps,
+				     Jacobian_sparsity& JS,
+				     bool print_msg_flag = false);
diff --git a/src/elliptic/make.code.defn b/src/elliptic/make.code.defn
index 3bc82f4..28308a0 100644
--- a/src/elliptic/make.code.defn
+++ b/src/elliptic/make.code.defn
@@ -2,7 +2,9 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Jacobian.cc gecon_wrapper.F77
+SRCS = Jacobian.cc \
+       dense_Jacobian.cc gecon_wrapper.F77 \
+       row_sparse_Jacobian.cc
 
 # Subdirectories containing source files
 SUBDIRS =