doc/html/hypre__UnitTest_8cpp_source.html

 /*@HEADER

 // ***********************************************************************

 //

 //       Ifpack: Object-Oriented Algebraic Preconditioner Package

 //                 Copyright (2002) Sandia Corporation

 //

 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive

 // license for use of this work by or on behalf of the U.S. Government.

 //

 // Redistribution and use in source and binary forms, with or without

 // modification, are permitted provided that the following conditions are

 // met:

 //

 // 1. Redistributions of source code must retain the above copyright

 // notice, this list of conditions and the following disclaimer.

 //

 // 2. Redistributions in binary form must reproduce the above copyright

 // notice, this list of conditions and the following disclaimer in the

 // documentation and/or other materials provided with the distribution.

 //

 // 3. Neither the name of the Corporation nor the names of the

 // contributors may be used to endorse or promote products derived from

 // this software without specific prior written permission.

 //

 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY

 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR

 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE

 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR

 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 //

 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)

 //

 // ***********************************************************************

 //@HEADER

 */


 #include "Ifpack.h"

 #include "Ifpack_Hypre.h"

 #include "AztecOO.h"

 #include "Galeri_Maps.h"

 #include "Galeri_CrsMatrices.h"

 #include "Galeri_Utils.h"

 #include "Epetra_MultiVector.h"


 #include "Teuchos_UnitTestHarness.hpp"

 #include "Teuchos_RCP.hpp"

 #include "Epetra_ConfigDefs.h"

 #include "Epetra_Vector.h"

 #include "Epetra_RowMatrix.h"

 #include "Epetra_MultiVector.h"

 #include "Epetra_CrsMatrix.h"

 #include "Epetra_Map.h"

 #ifdef HAVE_MPI

 #include "mpi.h"

 #include "Epetra_MpiComm.h"

 #else

 #include "Epetra_SerialComm.h"

 #endif


 #include "hypre_Helpers.hpp"


 #include "Teuchos_Array.hpp"

 #include <string>

 #include <stdio.h>

 #include <map>


 using Teuchos::RCP;

 using Teuchos::rcp;


 // Tests hypre interface's ability to initialize correctly

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, Construct ) {

   const int N = 10;


   Epetra_MpiComm comm(MPI_COMM_WORLD);

   Epetra_Map RowMap(N, 0, comm);

   Epetra_CrsMatrix Matrix(Copy, RowMap, 1);

   for(int i = 0; i < N; i++){

     int indices[1];

     double values[1];

     indices[0] = i;

     values[0] = 1.0;

     Matrix.InsertGlobalValues(i, 1, values, indices);

   }

   Matrix.FillComplete();

   Ifpack_Hypre preconditioner(&Matrix);

   TEST_EQUALITY(preconditioner.Initialize(),0);

 }


 // Tests hypre's ability to work when A has a funky row map, but the vectors have

 // a contiguous row map

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, ParameterList ){

   const double tol = 1e-7;


   Epetra_MpiComm Comm(MPI_COMM_WORLD);


   Epetra_Map*            Map;

   // pointer to the matrix to be created

   Epetra_CrsMatrix*      Matrix;

   // container for parameters

   Teuchos::ParameterList GaleriList;

   // here we specify the global dimension of the problem

   int nx = 10 * Comm.NumProc();

   int ny = 10 * Comm.NumProc();

   GaleriList.set("nx", nx);

   GaleriList.set("ny", ny);


   try

   {

     // Create the matrix

     Map = Galeri::CreateMap("Cartesian2D", Comm, GaleriList);

     Matrix   = Galeri::CreateCrsMatrix("Biharmonic2D", Map, GaleriList);


     // Create the parameter list

     Teuchos::ParameterList list("Preconditioner List");

     RCP<FunctionParameter> functs[11];

     functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 1000));               // max iterations

     functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, tol));                   // conv. tolerance

     functs[2] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTwoNorm, 1));                  // use the two norm as the stopping criteria

     functs[3] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetPrintLevel, 0));               // print solve info

     functs[4] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetLogging, 1));                  // needed to get run info later

     functs[5] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetPrintLevel, 1)); // print amg solution info

     functs[6] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetCoarsenType, 6));

     functs[7] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetRelaxType, 6));  //Sym G.S./Jacobi hybrid

     functs[8] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetNumSweeps, 1));

     functs[9] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetTol, 0.0));      // conv. tolerance zero

     functs[10] = rcp(new FunctionParameter(Preconditioner, &HYPRE_BoomerAMGSetMaxIter, 1));   //do only one iteration!


     list.set("Solver", PCG);

     list.set("Preconditioner", BoomerAMG);

     list.set("SolveOrPrecondition", Solver);

     list.set("SetPreconditioner", true);

     list.set("NumFunctions", 11);

     list.set<RCP<FunctionParameter>*>("Functions", functs);


     // Create the preconditioner

     Ifpack_Hypre preconditioner(Matrix);

     TEST_EQUALITY(preconditioner.SetParameters(list),0);

     TEST_EQUALITY(preconditioner.Compute(),0);


     // Create the RHS and solution vector

     int numVec = 2;

     Epetra_MultiVector X(preconditioner.OperatorDomainMap(), numVec);

     Epetra_MultiVector KnownX(preconditioner.OperatorDomainMap(), numVec);

     TEST_EQUALITY(KnownX.Random(),0);

     Epetra_MultiVector B(preconditioner.OperatorRangeMap(), numVec);

     TEST_EQUALITY(preconditioner.Apply(KnownX, B),0);


     TEST_EQUALITY(preconditioner.ApplyInverse(B,X),0);

     TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,Comm.NumProc())), true);


     //delete preconditioner;

     delete Map;

     delete Matrix;

   }

   catch (Galeri::Exception& rhs)

   {

     if (Comm.MyPID() == 0)

     {

       cerr << "Caught exception: ";

       rhs.Print();

     }

   }

 }


 // Tests the hypre interface's ability to work with both a preconditioner and linear

 // solver via ApplyInverse

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, Ifpack ){

   const double tol = 1E-7;


   //

   // Create Laplace2D with a contiguous row distribution

   //

   Epetra_CrsMatrix* Crs_Matrix = newCrsMatrix(4);


   //

   // Create the hypre preconditioner

   //

   Ifpack Factory;

   RCP<Ifpack_Preconditioner> preconditioner= rcp(Factory.Create("Hypre", Crs_Matrix));

   TEST_EQUALITY(preconditioner->Initialize(), 0);

   int NumProc = Crs_Matrix->Comm().NumProc();

   int MyPID = Crs_Matrix->Comm().MyPID();


   //

   // Create the solution vector and RHS

   //

   int numVec = 2;

   Epetra_MultiVector X(preconditioner->OperatorRangeMap(), numVec);

   Epetra_MultiVector KnownX(preconditioner->OperatorRangeMap(), numVec);

   Epetra_MultiVector B(preconditioner->OperatorDomainMap(), numVec);

   TEST_EQUALITY(KnownX.Random(), 0);

   TEST_EQUALITY(preconditioner->Apply(KnownX, B), 0);


   Teuchos::ParameterList list("New List");

   RCP<FunctionParameter> functs[5];

   functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 1000));

   functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, 1e-9));

   functs[2] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetLogging, 1));

   functs[3] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetPrintLevel, 0));

   functs[4] = rcp(new FunctionParameter(Preconditioner, &HYPRE_ParaSailsSetLogging, 0));

   list.set("NumFunctions", 5);

   list.set<RCP<FunctionParameter>*>("Functions", functs);

   list.set("SolveOrPrecondition", Solver);

   list.set("Solver", PCG);

   list.set("Preconditioner", ParaSails);

   list.set("SetPreconditioner", true);

   TEST_EQUALITY(preconditioner->SetParameters(list), 0);

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, PCG); // Use a PCG Solver

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Preconditioner, ParaSails); // Use a ParaSails Preconditioner

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGSetMaxIter, 1000); // Set maximum iterations

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGSetTol, 1e-9); // Set a tolerance

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver); // Solve the problem

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(true); // Use the preconditioner

   TEST_EQUALITY(preconditioner->Compute(),0);

   preconditioner->Condest(Ifpack_Cheap, 1000, 1e-9, Crs_Matrix);

   TEST_EQUALITY(preconditioner->ApplyInverse(B, X),0);

   TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,NumProc)), true);

   int numIters;

   double residual;

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGGetNumIterations, &numIters);

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGGetFinalRelativeResidualNorm, &residual);

   (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->CallFunctions();

   if(MyPID == 0) printf("It took %d iterations, and achieved %e residual.\n", numIters, residual);

   delete Crs_Matrix;

 }


 // This example uses contiguous maps, so hypre should not have problems

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, DiagonalMatrixInOrder ) {

   Epetra_MpiComm comm(MPI_COMM_WORLD);

   int numProcs = comm.NumProc();

   int myRank = comm.MyPID();


   //

   // Construct the contiguous map

   //

   Epetra_Map accMap(2*numProcs,0,comm);


   //

   // Construct the diagonal matrix

   //

   Epetra_CrsMatrix accMat(Copy,accMap,2,true);

   int col;

   double val;

   col = 2*myRank;

   val = 2*myRank+1;

   accMat.InsertGlobalValues(col,1,&val,&col);

   col = 2*myRank+1;

   val = 2*myRank+2;

   accMat.InsertGlobalValues(col,1,&val,&col);

   accMat.FillComplete();


   //

   // Create the initial guess and RHS

   //

   int numVec = 2;

   Epetra_MultiVector X(accMap, numVec);

   Epetra_MultiVector KnownX(accMap, numVec);

   KnownX.Random();

   Epetra_MultiVector B(accMap, numVec);

   accMat.Apply(KnownX, B);


   //

   // Create the parameter list

   //

   const double tol = 1e-7;

   Teuchos::ParameterList list("Preconditioner List");

   RCP<FunctionParameter> functs[5];

   functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 100));               // max iterations

   functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, tol));                   // conv. tolerance

   functs[2] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTwoNorm, 1));                  // use the two norm as the stopping criteria

   functs[3] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetPrintLevel, 2));               // print solve info

   functs[4] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetLogging, 1));

   list.set("Solver", PCG);

   list.set("SolveOrPrecondition", Solver);

   list.set("SetPreconditioner", false);

   list.set("NumFunctions", 5);

   list.set<RCP<FunctionParameter>*>("Functions", functs);


   //

   // Create the preconditioner (which is actually a PCG solver)

   //

   Ifpack_Hypre prec(&accMat);

   TEST_EQUALITY(prec.SetParameters(list),0);

   TEST_EQUALITY(prec.Compute(),0);


   //

   // Solve the linear system

   //

   TEST_EQUALITY(prec.ApplyInverse(B,X),0);


   TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,numProcs)), true);

 }


 // hypre does not like the distribution of the vectors in this example.

 // Our interface should detect that and remap as necessary.

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, DiagonalMatrixOutOfOrder ) {

   Epetra_MpiComm comm(MPI_COMM_WORLD);

   int numProcs = comm.NumProc();

   int myRank = comm.MyPID();


   //

   // Construct the map [0 0 1 1]

   // Note that the rows will be out of order

   //

   int elementList[2];

   elementList[0] = 2*myRank+1;

   elementList[1] = 2*myRank;

   Epetra_Map ncMap(2*numProcs,2,elementList,0,comm);


   //

   // Construct the diagonal matrix

   //

   Epetra_CrsMatrix accMat(Copy,ncMap,2,true);

   int col;

   double val;

   col = 2*myRank+1;

   val = 2*myRank+2;

   accMat.InsertGlobalValues(col,1,&val,&col);

   col = 2*myRank;

   val = 2*myRank+1;

   accMat.InsertGlobalValues(col,1,&val,&col);

   accMat.FillComplete();


   //

   // Create the parameter list

   //

   const double tol = 1e-7;

   Teuchos::ParameterList list("Preconditioner List");

   RCP<FunctionParameter> functs[5];

   functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 100));               // max iterations

   functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, tol));                   // conv. tolerance

   functs[2] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTwoNorm, 1));                  // use the two norm as the stopping criteria

   functs[3] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetPrintLevel, 2));               // print solve info

   functs[4] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetLogging, 1));

   list.set("Solver", PCG);

   list.set("SolveOrPrecondition", Solver);

   list.set("SetPreconditioner", false);

   list.set("NumFunctions", 5);

   list.set<RCP<FunctionParameter>*>("Functions", functs);


   //

   // Create the preconditioner (which is actually a PCG solver)

   //

   Ifpack_Hypre prec(&accMat);

   TEST_EQUALITY(prec.SetParameters(list),0);

   TEST_EQUALITY(prec.Compute(),0);


   //

   // Create the initial guess and RHS

   // Note that these use non-contiguous maps

   //

   int numVec = 2;

   Epetra_MultiVector X(ncMap, numVec);

   Epetra_MultiVector KnownX(ncMap, numVec);

   KnownX.Random();

   Epetra_MultiVector B(ncMap, numVec);

   accMat.Apply(KnownX, B);


   //

   // Solve the linear system

   //

   TEST_EQUALITY(prec.ApplyInverse(RHS,X),0);


   TEST_EQUALITY(EquivalentVectors(X, RHS, tol*10*pow(10.0,numProcs)), true);

 }


 // Creates the identity matrix with a non-contiguous row map

 // Even though the Epetra identity matrix has a map that hypre should not be happy with,

 // hypre should be able to redistribute it.  It should also be able to accept the

 // vectors we give it, since they're using the same distribution as the hypre matrix.

 // This tests hypre's ability to perform as a linear solver via ApplyInverse.

 TEUCHOS_UNIT_TEST( Ifpack_Hypre, NonContiguousRowMap ) {


   Epetra_MpiComm comm(MPI_COMM_WORLD);

   int numProcs = comm.NumProc();

   int myRank = comm.MyPID();


   //

   // Construct the map [0 1 0 1]

   //

   int elementList[2];

   elementList[0] = myRank;

   elementList[1] = myRank+numProcs;

   Epetra_Map badMap(2*numProcs,2,elementList,0,comm);

   Epetra_Map goodMap(2*numProcs,0,comm);


   //

   // Construct the identity matrix

   //

   Epetra_CrsMatrix badMat(Copy,badMap,2,true);

   int col;

   double val = 1.0;

   col = myRank;

   badMat.InsertGlobalValues(col,1,&val,&col);

   col = myRank+numProcs;

   badMat.InsertGlobalValues(col,1,&val,&col);

   badMat.FillComplete();


   //

   // Create the parameter list

   //

   const double tol = 1e-7;

   Teuchos::ParameterList list("Preconditioner List");

   RCP<FunctionParameter> functs[11];

   functs[0] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetMaxIter, 1));               // max iterations

   functs[1] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTol, tol));                   // conv. tolerance

   functs[2] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetTwoNorm, 1));                  // use the two norm as the stopping criteria

   functs[3] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetPrintLevel, 2));               // print solve info

   functs[4] = rcp(new FunctionParameter(Solver, &HYPRE_PCGSetLogging, 1));

   list.set("Solver", PCG);

   list.set("SolveOrPrecondition", Solver);

   list.set("SetPreconditioner", false);

   list.set("NumFunctions", 5);

   list.set<RCP<FunctionParameter>*>("Functions", functs);


   //

   // Create the preconditioner (which is actually a PCG solver)

   //

   Ifpack_Hypre prec(&badMat);

   TEST_EQUALITY(prec.SetParameters(list),0);

   TEST_EQUALITY(prec.Compute(),0);


   //

   // Create the initial guess and RHS

   //

   int numVec = 2;

   Epetra_MultiVector X(goodMap, numVec);

   Epetra_MultiVector RHS(goodMap, numVec);

   RHS.Random();


   //

   // Solve the linear system

   //

   TEST_EQUALITY(prec.ApplyInverse(RHS,X),0);


   TEST_EQUALITY(EquivalentVectors(X, RHS, tol*10*pow(10.0,numProcs)), true);

 }

Epetra_MultiVector

Ifpack_Preconditioner::Condest
virtual double Condest(const Ifpack_CondestType CT=Ifpack_Cheap, const int MaxIters=1550, const double Tol=1e-9, Epetra_RowMatrix *Matrix=0)=0
Computes the condition number estimate, returns its value.

Epetra_Map

Preconditioner

B

Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)

ParaSails

Epetra_CrsMatrix::InsertGlobalValues
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)

Epetra_CrsMatrix::Apply
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const

Epetra_MpiComm::NumProc
int NumProc() const

EquivalentVectors
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
Definition: hypre_Helpers.cpp:165

Teuchos::RCP::get
T * get() const

Epetra_Operator::OperatorDomainMap
virtual const Epetra_Map & OperatorDomainMap() const =0

Epetra_MpiComm

Ifpack_Preconditioner::Initialize
virtual int Initialize()=0
Computes all it is necessary to initialize the preconditioner.

Epetra_Comm::MyPID
virtual int MyPID() const =0

Epetra_CrsMatrix::FillComplete
int FillComplete(bool OptimizeDataStorage=true)

Solver
static bool Solver
Definition: performance.cpp:70

BoomerAMG

Ifpack_Hypre.h

Epetra_Operator::Apply
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0

Epetra_Operator::OperatorRangeMap
virtual const Epetra_Map & OperatorRangeMap() const =0

Ifpack_Preconditioner::SetParameters
virtual int SetParameters(Teuchos::ParameterList &List)=0
Sets all parameters for the preconditioner.

Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)

tol
const double tol
Definition: LL/hypre_UnitTest.cpp:75

hypre_Helpers.hpp

Ifpack_Preconditioner::ApplyInverse
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Applies the preconditioner to vector X, returns the result in Y.

Epetra_CrsMatrix::Print
virtual void Print(std::ostream &os) const

TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Ifpack_Hypre, Construct)
Definition: hypre_UnitTest.cpp:78

N
const int N
Definition: LL/hypre_UnitTest.cpp:73

Teuchos::ParameterList

Ifpack.h

E

Ifpack_Cheap
cheap estimate
Definition: Ifpack_CondestType.h:49

Epetra_Comm::NumProc
virtual int NumProc() const =0

Ifpack
Ifpack: a function class to define Ifpack preconditioners.
Definition: Ifpack.h:138

Ifpack::Create
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition: Ifpack.cpp:253

Epetra_CrsMatrix

Copy

PCG

Ifpack_Preconditioner::Compute
virtual int Compute()=0
Computes all it is necessary to apply the preconditioner.

Teuchos::RCP

TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)

newCrsMatrix
Epetra_CrsMatrix * newCrsMatrix(int N)
Definition: hypre_Helpers.cpp:120

Epetra_CrsMatrix::Comm
const Epetra_Comm & Comm() const

Epetra_MpiComm::MyPID
int MyPID() const

RHS
#define RHS(a)
Definition: MatGenFD.c:60