From e01a3c4dfd2cea40f95c7c49084ccc26d68f1a81 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 18 Jul 2002 17:40:09 +0000 Subject: storage for Jacobian matrix which knows ../util/ sparsity structure git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@631 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/elliptic/Jacobian.cc | 81 ++++++++++++++++++++++++++++++++ src/elliptic/Jacobian.hh | 120 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 201 insertions(+) create mode 100644 src/elliptic/Jacobian.cc create mode 100644 src/elliptic/Jacobian.hh (limited to 'src') diff --git a/src/elliptic/Jacobian.cc b/src/elliptic/Jacobian.cc new file mode 100644 index 0000000..1bd2b5b --- /dev/null +++ b/src/elliptic/Jacobian.cc @@ -0,0 +1,81 @@ +// Jacobian.cc -- data structures for the Jacobian matrix +// $Id$ + +// +// Jacobian::Jacobian +// +// dense_Jacobian::dense_Jacobian +// dense_Jacobian::zero +// dense_Jacobian::zero_row +// + +#include +#include +#include +#include + +#include "util_Table.h" +#include "cctk.h" +#include "cctk_Arguments.h" + +#include "stdc.h" +#include "config.hh" +#include "../jtutil/util.hh" +#include "../jtutil/array.hh" +#include "../jtutil/cpm_map.hh" +#include "../jtutil/linear_map.hh" +using jtutil::error_exit; + +#include "coords.hh" +#include "grid.hh" +#include "fd_grid.hh" +#include "patch.hh" +#include "patch_edge.hh" +#include "patch_interp.hh" +#include "ghost_zone.hh" +#include "patch_system.hh" +#include "Jacobian.hh" + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function constructs a Jacobian object. +// +Jacobian::Jacobian(const patch_system& ps) + : ps_(ps) +{ } + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function constructs a dense_Jacobian object. +// +dense_Jacobian::dense_Jacobian(const patch_system& ps) + : Jacobian(ps), + matrix_(0, ps.N_grid_points(), + 0, ps.N_grid_points()) +{ } + +//****************************************************************************** + +// +// This function zeros a dense_Jacobian object. +// +void dense_Jacobian::zero() +{ +// FIXME FIXME +} + +//****************************************************************************** + +// +// This function zeros a single row of a dense_Jacobian object. +// +void dense_Jacobian::zero_row(int II) +{ +// FIXME FIXME +} diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh new file mode 100644 index 0000000..2c64171 --- /dev/null +++ b/src/elliptic/Jacobian.hh @@ -0,0 +1,120 @@ +// Jacobian.hh -- data structures for the Jacobian matrix +// $Id$ + +// +// Jacobian -- abstract base class to describe a Jacobian matrix +// +// dense_Jacobian - Jacobian stored as a dense matrix +// + +// +// 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" +// + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// 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. +// +// A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid +// point number within the patch system. +// +// Note that the APIs here implicitly assume there is only a *single* gridfn +// in the Jacobian computation. (If we had N gridfns for this, then the +// Jacobian would really be a block matrix with N*N blocks.) This isn't +// very general, but matches our present use in this thorn. +// + +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 + + // 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) + 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); + } + + // zero the whole matrix or a single row + virtual void zero() = 0; + virtual void zero_row(int II) = 0; + + // basic meta-info + const patch_system& my_patch_system() const { return ps_; } + + // constructor, destructor + Jacobian(const patch_system& ps); + virtual ~Jacobian() { } + +private: + const patch_system& ps_; + }; + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This class stores the Jacobian as a dense matrix in Fortran (column) +// order. +// +class dense_Jacobian + : 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); } + + // zero the whole matrix or a single row + void zero(); + void zero_row(int II); + + // constructor, destructor +public: + dense_Jacobian(const patch_system& ps); + ~dense_Jacobian() { } + +private: + // subscripts are (JJ,II) + jtutil::array2d matrix_; + }; -- cgit v1.2.3