Richard.
2004Mar23 Richard Guenther <address@hidden>
* src/Transform/PETSc.h: new.
examples/Solvers/PETSc/PETSc.cpp: new.
examples/Solvers/PETSc/include.mk: new.
examples/Solvers/PETSc/makefile: new.
diff -Nru a/r2/examples/Solvers/PETSc/PETSc.cpp
b/r2/examples/Solvers/PETSc/PETSc.cpp
--- /dev/null Wed Dec 31 16:00:00 1969
+++ b/r2/examples/Solvers/PETSc/PETSc.cpp Tue Mar 23 10:25:34 2004
@@ -0,0 +1,100 @@
+// Example on how to use PETSc linear solvers from inside a Pooma program.
+// This solves the Poisson equation for a density distribution in a 3D
+// Array to create the corresponding gravitational potential.
+
+#include <iostream>
+#include "Transform/PETSc.h"
+
+#include "petscda.h"
+#include "petscksp.h"
+
+void doit()
+{
+ typedef Array<2, double, MultiPatch<GridTag, Remote<Brick> > > Array_t;
+ typedef Array_t::Layout_t::const_iterator layout_iter;
+
+ Interval<2> domain(8, 8);
+ Loc<2> blocks = makeRBlocks(domain, Pooma::contexts());
+ GridLayout<2> layout(domain, blocks, GuardLayers<2>(1), DistributedTag());
+
+ // Create DA
+ Pooma::PoomaDA da(layout, DA_NONPERIODIC, DA_STENCIL_STAR, 1);
+
+ // Assemble matrix
+ Mat jac;
+ DAGetMatrix(da, MATMPIAIJ, &jac);
+ {
+ int mx, my;
+ PetscScalar Hx, Hy, HxdHy, HydHx;
+ DAGetInfo(da,0,&mx,&my,PETSC_NULL,0,0,0,0,0,0,0);
+ Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1);
+ HxdHy = Hx/Hy; HydHx = Hy/Hx;
+ int xs, ys, xm, ym;
+ MatStencil row,col[5];
+ DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
+ PetscScalar v[5];
+ for (int j=ys; j<ys+ym; j++){
+ for(int i=xs; i<xs+xm; i++){
+ row.i = i; row.j = j;
+ if (i==0 || j==0 || i==mx-1 || j==my-1) {
+ v[0] = 2.0*(HxdHy + HydHx);
+ MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
+ } else {
+ v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
+ v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
+ v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
+ v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
+ v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
+ MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
+ }
+ }
+ }
+ }
+ MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
+
+ // problem
+ Array_t rh(layout), ph(layout);
+ rh(rh.totalDomain()) = 1.0;
+ ph(ph.totalDomain()) = 0.0;
+ Pooma::pinfo << "Density distribution:\n" << rh << std::endl;
+
+ // assemble rhs and solution
+ Vec pph, prh;
+ DACreateGlobalVector(da, &prh);
+ VecDuplicate(prh, &pph);
+
+ da.assign(prh, rh);
+
+ // create solver
+ KSP ksp;
+ KSPCreate(PETSC_COMM_WORLD, &ksp);
+ KSPSetTolerances(ksp, 1e-10, 1e-5, 1e5, 10);
+ KSPSetOperators(ksp, jac, jac, DIFFERENT_NONZERO_PATTERN);
+ KSPSetFromOptions(ksp);
+
+ // solve and copy solution
+ KSPSetRhs(ksp, prh);
+ KSPSetSolution(ksp, pph);
+ KSPSolve(ksp);
+ da.assign(ph, pph);
+
+ // show soultion
+ Pooma::pinfo << "Solution:\n" << ph << std::endl;
+}
+
+int main(int argc, char **argv)
+{
+ Pooma::initialize(argc, argv);
+
+ PetscSetCommWorld(MPI_COMM_WORLD);
+ PetscInitialize(&argc, &argv, NULL, NULL);
+
+ doit();
+
+ // cleanup
+ PetscFinalize();
+ Pooma::finalize();
+ return 0;
+}
+
diff -Nru a/r2/examples/Solvers/PETSc/include.mk
b/r2/examples/Solvers/PETSc/include.mk
--- /dev/null Wed Dec 31 16:00:00 1969
+++ b/r2/examples/Solvers/PETSc/include.mk Tue Mar 23 10:25:34 2004
@@ -0,0 +1,59 @@
+# Generated by mm.pl: Mon Mar 9 13:58:39 MST 1998
+# ACL:license
+# ----------------------------------------------------------------------
+# This software and ancillary information (herein called "SOFTWARE")
+# called POOMA (Parallel Object-Oriented Methods and Applications) is
+# made available under the terms described here. The SOFTWARE has been
+# approved for release with associated LA-CC Number LA-CC-98-65.
+#
+# Unless otherwise indicated, this SOFTWARE has been authored by an
+# employee or employees of the University of California, operator of the
+# Los Alamos National Laboratory under Contract No. W-7405-ENG-36 with
+# the U.S. Department of Energy. The U.S. Government has rights to use,
+# reproduce, and distribute this SOFTWARE. The public may copy, distribute,
+# prepare derivative works and publicly display this SOFTWARE without
+# charge, provided that this Notice and any statement of authorship are
+# reproduced on all copies. Neither the Government nor the University
+# makes any warranty, express or implied, or assumes any liability or
+# responsibility for the use of this SOFTWARE.
+#
+# If SOFTWARE is modified to produce derivative works, such modified
+# SOFTWARE should be clearly marked, so as not to confuse it with the
+# version available from LANL.
+#
+# For more information about POOMA, send e-mail to address@hidden,
+# or visit the POOMA web page at http://www.acl.lanl.gov/pooma/.
+# ----------------------------------------------------------------------
+# ACL:license
+
+
+# Wrap make components from SHARED_ROOT and the current directory in the
+# proper order so that variables like ODIR have the correct directory-specific
+# value at the right moment. See the included files for details of what they
+# are doing. This file should NOT be manually edited.
+
+# Set NEXTDIR, THISDIR and DIR_LIST
+include $(SHARED_ROOT)/include1.mk
+
+# Include list of subdirectories to process
+-include $(THISDIR)/subdir.mk
+
+# Set ODIR, PROJECT_INCLUDES, UNIQUE
+include $(SHARED_ROOT)/include2.mk
+
+# Set list of object files, relative to ODIR
+-include $(THISDIR)/objfile.mk
+
+# Set rules for the ODIR directory
+include $(SHARED_ROOT)/compilerules.mk
+
+# Remove current dir from DIR_LIST
+DIR_LIST :=$(filter-out $(firstword $(DIR_LIST)), $(DIR_LIST))
+
+
+# ACL:rcsinfo
+# ----------------------------------------------------------------------
+# $RCSfile: include.mk,v $ $Author: swhaney $
+# $Revision: 1.3 $ $Date: 2000/03/07 13:15:37 $
+# ----------------------------------------------------------------------
+# ACL:rcsinfo
diff -Nru a/r2/examples/Solvers/PETSc/makefile
b/r2/examples/Solvers/PETSc/makefile
--- /dev/null Wed Dec 31 16:00:00 1969
+++ b/r2/examples/Solvers/PETSc/makefile Tue Mar 23 10:25:34 2004
@@ -0,0 +1,65 @@
+# Generated by mm.pl: Mon Mar 9 13:58:39 MST 1998
+# ACL:license
+# ----------------------------------------------------------------------
+# This software and ancillary information (herein called "SOFTWARE")
+# called POOMA (Parallel Object-Oriented Methods and Applications) is
+# made available under the terms described here. The SOFTWARE has been
+# approved for release with associated LA-CC Number LA-CC-98-65.
+#
+# Unless otherwise indicated, this SOFTWARE has been authored by an
+# employee or employees of the University of California, operator of the
+# Los Alamos National Laboratory under Contract No. W-7405-ENG-36 with
+# the U.S. Department of Energy. The U.S. Government has rights to use,
+# reproduce, and distribute this SOFTWARE. The public may copy, distribute,
+# prepare derivative works and publicly display this SOFTWARE without
+# charge, provided that this Notice and any statement of authorship are
+# reproduced on all copies. Neither the Government nor the University
+# makes any warranty, express or implied, or assumes any liability or
+# responsibility for the use of this SOFTWARE.
+#
+# If SOFTWARE is modified to produce derivative works, such modified
+# SOFTWARE should be clearly marked, so as not to confuse it with the
+# version available from LANL.
+#
+# For more information about POOMA, send e-mail to address@hidden,
+# or visit the POOMA web page at http://www.acl.lanl.gov/pooma/.
+# ----------------------------------------------------------------------
+# ACL:license
+
+# This file is user-editable
+
+PROJECT_ROOT = $(shell cd ../../..; pwd)
+include $(PROJECT_ROOT)/config/head.mk
+
+PASS=APP
+
+default:: PETSc
+
+.PHONY: PETSc
+
+PETSc:: $(ODIR)/PETSc
+
+include $(SHARED_ROOT)/tail.mk
+
+$(ODIR)/PETSc.o: PETSc.cpp
+ $(CXX) -o $@ -c $< $(RULE_CXXOPTS) $(RULE_INCLUDES) $(SUITE_DEFINES) \
+ -DPETSC_USE_EXTERN_CXX \
+ -I$(PETSC_DIR)/bmake/linux-gnu \
+ -I$(PETSC_DIR)/include
+
+$(ODIR)/PETSc: $(ODIR)/PETSc.o
+ $(CXX) -o $@ $^ $(RULE_LD_OPTS) \
+ -L$(PETSC_DIR)/lib/libg/linux-gnu \
+ -lpetscmat -lpetscksp -lpetscdm -lpetscvec -lpetscmat -lpetscdm
-lpetsc -lpetscts -lpetscvec -lpetscsnes -lblas -llapack -lblas -lg2c -ldl
-static
+
+# ACL:rcsinfo
+# ----------------------------------------------------------------------
+# $RCSfile: makefile,v $ $Author: julianc $
+# $Revision: 1.6 $ $Date: 2000/07/21 20:43:14 $
+# ----------------------------------------------------------------------
+# ACL:rcsinfo
diff -Nru a/r2/src/Transform/PETSc.h b/r2/src/Transform/PETSc.h
--- /dev/null Wed Dec 31 16:00:00 1969
+++ b/r2/src/Transform/PETSc.h Tue Mar 23 10:25:34 2004
@@ -0,0 +1,396 @@
+#ifndef POOMA_TRANSFORM_PETSC_H
+#define POOMA_TRANSFORM_PETSC_H
+
+// PETSc interfacing with POOMA using the PETSc DA interface.
+//
+// Copyright (c) 2004 by Richard Guenther <address@hidden>
+//
+// This file is in the public domain.
+
+/** @file
+ * @ingroup Utilities
+ * @brief
+ * Interfacing with the PETSc library of (non-)linear solvers.
+ *
+ * Interfacing supports the PETSc DA (distributed arrays) notion
+ * for creating (non-)linear solvers for implicit finite difference
+ * methods. Using this wrappers you can fill your right-hand-side
+ * vector from a POOMA engine and transfer the result-vector to
+ * a POOMA engine.
+ *
+ * You are going to use the PetscDA class and its methods.
+ * See examples/Solver/PETSc for how to use this.
+ */
+
+#include "Pooma/Arrays.h"
+#include "petscda.h"
+
+
+template <class MeshTag, class T, class EngineTag>
+class Field;
+
+
+namespace Pooma {
+
+
+/**
+ * Helper to convert DALocalInfo domain info to appropriate
+ * Pooma Interval
+ */
+
+template <int Dim>
+struct PoomaDAGetDomain;
+
+template <>
+struct PoomaDAGetDomain<1> {
+ static inline
+ Interval<1> innerDomain(DALocalInfo& i)
+ {
+ return Interval<1>(i.xs, i.xs+i.xm-1);
+ }
+ static inline
+ Interval<1> totalDomain(DALocalInfo& i)
+ {
+ return Interval<1>(i.gxs, i.gxs+i.gxm-1);
+ }
+};
+
+template <>
+struct PoomaDAGetDomain<2> {
+ static inline
+ Interval<2> innerDomain(DALocalInfo& i)
+ {
+ return Interval<2>(Interval<1>(i.xs, i.xs+i.xm-1),
+ Interval<1>(i.ys, i.ys+i.ym-1));
+ }
+ static inline
+ Interval<2> totalDomain(DALocalInfo& i)
+ {
+ return Interval<2>(Interval<1>(i.gxs, i.gxs+i.gxm-1),
+ Interval<1>(i.gys, i.gys+i.gym-1));
+ }
+};
+
+template <>
+struct PoomaDAGetDomain<3> {
+ static inline
+ Interval<3> innerDomain(DALocalInfo& i)
+ {
+ return Interval<3>(Interval<1>(i.xs, i.xs+i.xm-1),
+ Interval<1>(i.ys, i.ys+i.ym-1),
+ Interval<1>(i.zs, i.zs+i.zm-1));
+ }
+ static inline
+ Interval<3> totalDomain(DALocalInfo& i)
+ {
+ return Interval<3>(Interval<1>(i.gxs, i.gxs+i.gxm-1),
+ Interval<1>(i.gys, i.gys+i.gym-1),
+ Interval<1>(i.gzs, i.gzs+i.gzm-1));
+ }
+};
+
+
+
+/**
+ * Helper to ease brick-engine -> vector copy
+ */
+
+template <int Dim>
+struct PoomaDACopy;
+
+template <>
+struct PoomaDACopy<1> {
+ template <class T>
+ static
+ void copy(Vec v, const Engine<1, T, Brick>& e)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<1> d(e.domain());
+ for (int I=d.first(); I<=d.last(); ++I)
+ pa[idx++] = e(I);
+ VecRestoreArray(v, &pa);
+ }
+ template <class T>
+ static
+ void copy(const Engine<1, T, Brick>& e, Vec v)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<1> d(e.domain());
+ for (int I=d.first(); I<=d.last(); ++I)
+ e(I) = pa[idx++];
+ VecRestoreArray(v, &pa);
+ }
+};
+
+template <>
+struct PoomaDACopy<2> {
+ template <class T>
+ static
+ void copy(Vec v, const Engine<2, T, Brick>& e)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<2> d(e.domain());
+ for (int J=d[1].first(); J<=d[1].last(); ++J)
+ for (int I=d[0].first(); I<=d[0].last(); ++I)
+ pa[idx++] = e(I, J);
+ VecRestoreArray(v, &pa);
+ }
+ template <class T>
+ static
+ void copy(const Engine<2, T, Brick>& e, Vec v)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<2> d(e.domain());
+ for (int J=d[1].first(); J<=d[1].last(); ++J)
+ for (int I=d[0].first(); I<=d[0].last(); ++I)
+ e(I, J) = pa[idx++];
+ VecRestoreArray(v, &pa);
+ }
+};
+
+template <>
+struct PoomaDACopy<3> {
+ template <class T>
+ static
+ void copy(Vec v, const Engine<3, T, Brick>& e)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<3> d(e.domain());
+ for (int K=d[2].first(); K<=d[2].last(); ++K)
+ for (int J=d[1].first(); J<=d[1].last(); ++J)
+ for (int I=d[0].first(); I<=d[0].last(); ++I)
+ pa[idx++] = e(I, J, K);
+ VecRestoreArray(v, &pa);
+ }
+ template <class T>
+ static
+ void copy(const Engine<3, T, Brick>& e, Vec v)
+ {
+ PetscScalar *pa;
+ VecGetArray(v, &pa);
+ int idx=0;
+ Interval<3> d(e.domain());
+ for (int K=d[2].first(); K<=d[2].last(); ++K)
+ for (int J=d[1].first(); J<=d[1].last(); ++J)
+ for (int I=d[0].first(); I<=d[0].last(); ++I)
+ e(I, J, K) = pa[idx++];
+ VecRestoreArray(v, &pa);
+ }
+};
+
+
+
+/**
+ * Struct to wrap extra global information about a DA.
+ */
+
+struct PoomaDA {
+
+ /// Creates a PETSc DA from the specified layout.
+ /// Extra arguments are like DACreateNd, namely the periodicity
+ /// and stencil type and the stencil width.
+
+ template <class Layout>
+ PoomaDA(const Layout &l, DAPeriodicType pt, DAStencilType st, int sw);
+
+ ~PoomaDA()
+ {
+ delete[] info;
+ DADestroy(da);
+ }
+
+
+ /// Can use this as PETSc DA type.
+
+ operator DA() const { return da; }
+
+
+ /// Assign from POOMA engine to PETSc vector.
+
+ template <int Dim, class T, class EngineTag>
+ void assign(Vec v, const Engine<Dim, T, EngineTag> &e);
+
+ /// Assign from POOMA array to PETSc vector.
+
+ template <int Dim, class T, class EngineTag>
+ void assign(Vec v, const Array<Dim, T, EngineTag> &a)
+ {
+ this->assign(v, a.engine());
+ }
+
+ /// Assign from POOMA field to PETSc vector.
+
+ template <class MeshTag, class T, class EngineTag>
+ void assign(Vec v, const Field<MeshTag, T, EngineTag> &f)
+ {
+ this->assign(v, f.fieldEngine().engine());
+ }
+
+
+ /// Assign from PETSc vector to POOMA engine.
+
+ template <int Dim, class T, class EngineTag>
+ void assign(const Engine<Dim, T, EngineTag> &e, Vec v);
+
+ /// Assign from PETSc vector to POOMA array.
+
+ template <int Dim, class T, class EngineTag>
+ void assign(const Array<Dim, T, EngineTag> &a, Vec v)
+ {
+ this->assign(a.engine(), v);
+ }
+
+ /// Assign from PETSc vector to POOMA field.
+
+ template <class MeshTag, class T, class EngineTag>
+ void assign(const Field<MeshTag, T, EngineTag> &f, Vec v)
+ {
+ this->assign(f.fieldEngine().engine(), v);
+ }
+
+
+private:
+ DA da;
+ int n;
+ DALocalInfo *info;
+
+};
+
+
+template <class Layout>
+PoomaDA::PoomaDA(const Layout &l, DAPeriodicType pt, DAStencilType st, int sw)
+{
+ enum { Dim = Layout::dimensions };
+
+ // create DA
+ if (Dim == 1) {
+ DACreate1d(PETSC_COMM_WORLD, /* MPI communicator */
+ pt, /* grid periodicity */
+ l.innerDomain()[0].size(), /* global domain size */
+ 1, /* degrees of freedom */
+ sw, /* stencil width */
+ PETSC_NULL, /* local domain sizes per dimension */
+ &this->da);
+ } else if (Dim == 2) {
+ DACreate2d(PETSC_COMM_WORLD, /* MPI communicator */
+ pt, /* grid periodicity */
+ st, /* stencil type */
+ l.innerDomain()[0].size(),
+ l.innerDomain()[1].size(), /* global domain size */
+ PETSC_DECIDE, PETSC_DECIDE,/* # processors */
+ 1, /* degrees of freedom */
+ sw, /* stencil width */
+ PETSC_NULL, PETSC_NULL, /* local domain sizes per dimension */
+ &this->da);
+ } else if (Dim == 3) {
+ DACreate3d(PETSC_COMM_WORLD, /* MPI communicator */
+ pt, /* grid periodicity */
+ st, /* stencil type */
+ l.innerDomain()[0].size(), l.innerDomain()[1].size(),
+ l.innerDomain()[2].size(), /* global domain size */
+ PETSC_DECIDE, PETSC_DECIDE,
+ PETSC_DECIDE, /* # processors */
+ 1, /* degrees of freedom */
+ sw, /* stencil width */
+ PETSC_NULL, PETSC_NULL,
+ PETSC_NULL, /* local domain sizes per dimension */
+ &this->da);
+ }
+
+ // collect local information
+ int m, n, p;
+ DAGetInfo(this->da, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL,
+ &m, &n, &p,
+ PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);
+ this->n = m*n*p;
+ PInsist(Pooma::contexts() == this->n, "nr patches");
+ this->info = new DALocalInfo[this->n];
+ DAGetLocalInfo(this->da, &this->info[Pooma::context()]);
+
+ // distribute local information
+ // fixme - MPI_Allgather wrapper missing
+#if POOMA_MPI
+ MPI_Allgather(&this->info[Pooma::context()], sizeof(DALocalInfo), MPI_CHAR,
+ this->info, sizeof(DALocalInfo), MPI_CHAR,
+ MPI_COMM_WORLD);
+#endif
+}
+
+template <int Dim, class T, class EngineTag>
+void PoomaDA::assign(Vec v, const Engine<Dim, T, EngineTag> &e)
+{
+ typedef Engine<Dim, T, EngineTag> Engine_t;
+ typedef typename NewEngine<Engine_t, Interval<Dim> >::Type_t ViewEngine_t;
+
+ // our local brick engine
+ Engine<Dim, T, Brick> local_e;
+ Interval<Dim> local_I;
+
+ // loop over all DA patches
+ for (int i=0; i<this->n; ++i) {
+ // create DA patch context local pooma array
+ Interval<Dim>
lPatch(PoomaDAGetDomain<Dim>::innerDomain(this->info[i]));
+ Array<Dim, T, Remote<Brick> > a;
+ a.engine() = Engine<Dim, T, Remote<Brick> >(i, lPatch);
+ Array<Dim, T, typename ViewEngine_t::Tag_t> e_array(ViewEngine_t(e,
lPatch));
+ a = e_array;
+
+ // remember local engine
+ if (i == Pooma::context()) {
+ local_e = a.engine().localEngine();
+ local_I = lPatch;
+ }
+ }
+ Pooma::blockAndEvaluate();
+
+ // do the copy locally
+ PoomaDACopy<Dim>::copy(v, local_e);
+}
+
+template <int Dim, class T, class EngineTag>
+void PoomaDA::assign(const Engine<Dim, T, EngineTag> &e, Vec v)
+{
+ typedef Engine<Dim, T, EngineTag> Engine_t;
+ typedef typename NewEngine<Engine_t, Interval<Dim> >::Type_t ViewEngine_t;
+
+ // our local brick engine
+ Interval<Dim>
local_I(PoomaDAGetDomain<Dim>::innerDomain(this->info[Pooma::context()]));
+ Engine<Dim, T, Brick> local_e(local_I);
+
+ // copy into the local brick
+ // if it were not the different array index ordering we could construct
+ // the brick engine with external data and avoid the double copying
+ PoomaDACopy<Dim>::copy(local_e, v);
+
+ // loop over all DA patches
+ for (int i=0; i<this->n; ++i) {
+ // create DA patch context local pooma array
+ Interval<Dim> lPatch(PoomaDAGetDomain<Dim>::innerDomain(this->info[i]));
+ Array<Dim, T, Remote<Brick> > a;
+ a.engine() = Engine<Dim, T, Remote<Brick> >(i, lPatch);
+
+ // override local engine with our one
+ if (Pooma::context() == i)
+ a.engine().localEngine() = local_e;
+
+ // distribute the copy
+ Array<Dim, T, typename ViewEngine_t::Tag_t> e_array;
+ e_array.engine() = ViewEngine_t(e, lPatch);
+ e_array = a;
+ }
+}
+
+
+} // namespace Pooma
+
+#endif