Added pardiso test (Thanks to Bernd Flemisch)

[[Imported from SVN: r861]]

istl/Makefile.am

@@ -2,7 +2,11 @@
 
 SUBDIRS = . tutorial test paamg
 
-SOURCES = allocator.hh basearray.hh bcrsmatrix.hh bvector.hh gsetc.hh io.hh istlexception.hh vbcrsmatrix.hh vbvector.hh ilu.hh operators.hh preconditioners.hh solvers.hh indexset.hh communicator.hh remoteindices.hh mpitraits.hh interface.hh indicessyncer.hh matrixindexset.hh selection.hh plocalindex.hh localindex.hh scalarproducts.hh solvercategory.hh bdmatrix.hh matrixutils.hh schwarz.hh owneroverlapcopy.hh matrix.hh supermatrix.hh superlu.hh
+SOURCES = allocator.hh basearray.hh bcrsmatrix.hh bvector.hh gsetc.hh io.hh istlexception.hh vbcrsmatrix.hh vbvector.hh \
+	ilu.hh operators.hh preconditioners.hh solvers.hh indexset.hh communicator.hh remoteindices.hh mpitraits.hh \
+	interface.hh indicessyncer.hh matrixindexset.hh selection.hh plocalindex.hh localindex.hh scalarproducts.hh \
+	solvercategory.hh bdmatrix.hh matrixutils.hh schwarz.hh owneroverlapcopy.hh matrix.hh supermatrix.hh superlu.hh \
+	pardiso.hh
 
 istldir = $(includedir)/dune/istl
 istl_HEADERS = $(SOURCES)

istl/pardiso.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_PARDISO_HH
+#define DUNE_PARDISO_HH
+
+#include <dune/istl/preconditioners.hh>
+
+/* Change this, if your Fortran compiler does not append underscores. */
+/* e.g. the AIX compiler:  #define F77_FUNC(func) func                */
+
+#ifdef AIX
+#define F77_FUNC(func)  func
+#else
+#define F77_FUNC(func)  func ## _
+#endif
+
+
+#ifdef HAVE_PARDISO
+/* PARDISO prototype. */
+extern "C" int F77_FUNC(pardisoinit)
+  (void *, int *, int *);
+
+extern "C" int F77_FUNC(pardiso)
+  (void *, int *, int *, int *, int *, int *,
+  double *, int *, int *, int *, int *, int *,
+  int *, double *, double *, int *);
+#endif
+
+namespace Dune {
+
+
+  /*! \brief The sequential Pardiso preconditioner.
+
+     Put the Pardiso direct solver into the preconditioner framework.
+   */
+  template<class M, class X, class Y>
+  class SeqPardiso : public Preconditioner<X,Y> {
+  public:
+    //! \brief The matrix type the preconditioner is for.
+    typedef M matrix_type;
+    //! \brief The domain type of the preconditioner.
+    typedef X domain_type;
+    //! \brief The range type of the preconditioner.
+    typedef Y range_type;
+    //! \brief The field type of the preconditioner.
+    typedef typename X::field_type field_type;
+
+    typedef typename M::RowIterator RowIterator;
+    typedef typename M::ColIterator ColIterator;
+
+    // define the category
+    enum {
+      //! \brief The category the preconditioner is part of
+      category=SolverCategory::sequential
+    };
+
+    /*! \brief Constructor.
+
+       Constructor gets all parameters to operate the prec.
+       \param A The matrix to operate on.
+       \param n The number of iterations to perform.
+       \param w The relaxation factor.
+     */
+    SeqPardiso (const M& A)
+      : A_(A)
+    {
+#ifdef HAVE_PARDISO
+
+      mtype_ = 11;
+      nrhs_ = 1;
+      num_procs_ = 1;
+      maxfct_ = 1;
+      mnum_   = 1;
+      msglvl_ = 0;
+      error_  = 0;
+
+      n_ = A_.rowdim();
+      int nnz = 0;
+      RowIterator endi = A_.end();
+      for (RowIterator i = A_.begin(); i != endi; ++i)
+      {
+        if (A_.rowdim(i.index()) != 1)
+          DUNE_THROW(NotImplemented, "SeqPardiso: row blocksize != 1.");
+        ColIterator endj = (*i).end();
+        for (ColIterator j = (*i).begin(); j != endj; ++j) {
+          if (A_.coldim(j.index()) != 1)
+            DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
+          nnz++;
+        }
+      }
+
+      std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
+
+      a_ = new double[nnz];
+      ia_ = new int[n_+1];
+      ja_ = new int[nnz];
+
+      int count = 0;
+      for (RowIterator i = A_.begin(); i != endi; ++i)
+      {
+        ia_[i.index()] = count+1;
+        ColIterator endj = (*i).end();
+        for (ColIterator j = (*i).begin(); j != endj; ++j) {
+          a_[count] = *j;
+          ja_[count] = j.index()+1;
+
+          count++;
+        }
+      }
+      ia_[n_] = count+1;
+
+      F77_FUNC(pardisoinit) (pt_,  &mtype_, iparm_);
+
+      int phase = 11;
+      int idum;
+      double ddum;
+      iparm_[2]  = num_procs_;
+
+      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
+                         &n_, a_, ia_, ja_, &idum, &nrhs_,
+                         iparm_, &msglvl_, &ddum, &ddum, &error_);
+
+      if (error_ != 0)
+        DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
+
+      std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
+
+#else
+      DUNE_THROW(NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
+#endif
+    }
+
+    /*!
+       \brief Prepare the preconditioner.
+
+       \copydoc Preconditioner::pre(X&,Y&)
+     */
+    virtual void pre (X& x, Y& b) {}
+
+    /*!
+       \brief Apply the preconditioner.
+
+       \copydoc Preconditioner::apply(X&,const Y&)
+     */
+    virtual void apply (X& v, const Y& d)
+    {
+#ifdef HAVE_PARDISO
+      int phase = 33;
+
+      iparm_[7] = 1;         /* Max numbers of iterative refinement steps. */
+      int idum;
+
+      double x[n_];
+      double b[n_];
+      for (int i = 0; i < n_; i++) {
+        x[i] = v[i];
+        b[i] = d[i];
+      }
+
+      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
+                         &n_, a_, ia_, ja_, &idum, &nrhs_,
+                         iparm_, &msglvl_, b, x, &error_);
+
+      if (error_ != 0)
+        DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
+
+      for (int i = 0; i < n_; i++)
+        v[i] = x[i];
+
+      std::cout << "SeqPardiso: Backsolve completed." << std::endl;
+#endif
+    }
+
+    /*!
+       \brief Clean up.
+
+       \copydoc Preconditioner::post(X&)
+     */
+    virtual void post (X& x) {}
+
+    ~SeqPardiso()
+    {
+#ifdef HAVE_PARDISO
+      int phase = -1;                   // Release internal memory.
+      int idum;
+      double ddum;
+
+      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
+                         &n_, &ddum, ia_, ja_, &idum, &nrhs_,
+                         iparm_, &msglvl_, &ddum, &ddum, &error_);
+#endif
+    }
+
+  private:
+    M A_; //!< The matrix we operate on.
+    int n_; //!< dimension of the system
+    double *a_; //!< matrix values
+    int *ia_; //!< indices to rows
+    int *ja_; //!< column indices
+    int mtype_; //!< matrix type, currently only 11 (real unsymmetric matrix) is supported
+    int nrhs_; //!< number of right hand sides
+    void *pt_[64]; //!< internal solver memory pointer
+    int iparm_[64]; //!< Pardiso control parameters.
+    int maxfct_;        //!< Maximum number of numerical factorizations.
+    int mnum_;  //!<        Which factorization to use.
+    int msglvl_;    //!< flag to print statistical information
+    int error_;      //!< error flag
+    int num_procs_; //!< number of processors.
+  };
+
+}
+
+
+
+
+
+
+#endif

istl/test/Makefile.am

@@ -9,15 +9,19 @@ if SUPERLU
   SUPERLUTESTS =  superlutest overlappingschwarztest
 endif
 
+if PARDISO
+  PARDISOTEST = test_pardiso
+endif
+
 # which tests where program to build and run are equal
 NORMALTESTS = matrixutilstest matrixtest bvectortest indexsettest \
 	bcrsbuildtest matrixiteratortest
 
 # list of tests to run (indicestest is special case)
-TESTS = $(NORMALTESTS) $(MPITESTS)
+TESTS = $(NORMALTESTS) $(MPITESTS) $(SUPERLUTESTS) $(PARDISOTEST)
 
 # programs just to build when "make check" is used
-check_PROGRAMS = $(NORMALTESTS) $(MPITESTS) $(SUPERLUTESTS)
+check_PROGRAMS = $(NORMALTESTS) $(MPITESTS) $(SUPERLUTESTS) $(PARDISOTEST)
 
 # define the programs
 
@@ -32,6 +36,12 @@ if SUPERLU
   overlappingschwarztest_LDFLAGS= $(AM_LDFLAGS) $(SUPERLU_LDFLAGS)
   overlappingschwarztest_CPPFLAGS=$(AM_CPPFLAGS) $(SUPERLU_CPPFLAGS)
 endif
+
+if PARDISO
+  test_pardiso_SOURCES = test_pardiso.cc
+  test_pardiso_LDADD = $(PARDISO_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+endif
+
 bcrsbuildtest_SOURCES = bcrsbuild.cc
 
 bvectortest_SOURCES = bvectortest.cc

istl/test/test_pardiso.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "config.h"
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/pardiso.hh>
+
+int main(int argc, char** argv)
+{
+  try
+  {
+#ifdef HAVE_PARDISO
+    /* Matrix data. */
+    int n = 8;
+    int ia[ 9] = { 0, 4, 7, 9, 11, 12, 15, 17, 20 };
+    int ja[20] = { 0,    2,       5, 6,
+                   1, 2,    4,
+                   2,             7,
+                   3,       6,
+                   4,
+                   2,       5,    7,
+                   1,             6,
+                   2,          6, 7 };
+    double a[20] = { 7.0,      1.0,           2.0, 7.0,
+                     -4.0, 8.0,      2.0,
+                     1.0,                     5.0,
+                     7.0,           9.0,
+                     -4.0,
+                     7.0,           3.0,      8.0,
+                     1.0,                    11.0,
+                     -3.0,                2.0, 5.0 };
+
+    int nnz = ia[n];
+    int mtype = 11;                     /* Real unsymmetric matrix */
+
+    /* RHS and solution vectors. */
+    double b[8], x[8];
+    int nrhs = 1;                       /* Number of right hand sides. */
+
+    /* Internal solver memory pointer pt,                  */
+    /* 32-bit: int pt[64]; 64-bit: long int pt[64]         */
+    /* or void *pt[64] should be OK on both architectures  */
+    void    *pt[64];
+
+    /* Pardiso control parameters. */
+    int iparm[64];
+    int maxfct, mnum, phase, error, msglvl;
+
+    /* Number of processors. */
+    int num_procs;
+
+    /* Auxiliary variables. */
+    char    *var;
+    int i;
+
+    double ddum;                        /* Double dummy */
+    int idum;                           /* Integer dummy. */
+
+    /* -------------------------------------------------------------------- */
+    /* ..  Setup Pardiso control parameters.                                */
+    /* -------------------------------------------------------------------- */
+
+    F77_FUNC(pardisoinit) (pt,  &mtype, iparm);
+
+    iparm[2]  = 1;
+
+    maxfct = 1;         /* Maximum number of numerical factorizations.  */
+    mnum   = 1;         /* Which factorization to use. */
+
+    msglvl = 0;         /* Print statistical information  */
+    error  = 0;         /* Initialize error flag */
+
+    /* -------------------------------------------------------------------- */
+    /* ..  Convert matrix from 0-based C-notation to Fortran 1-based        */
+    /*     notation.                                                        */
+    /* -------------------------------------------------------------------- */
+    for (i = 0; i < n+1; i++) {
+      ia[i] += 1;
+    }
+    for (i = 0; i < nnz; i++) {
+      ja[i] += 1;
+    }
+
+    phase = 13;
+
+    iparm[7] = 1;       /* Max numbers of iterative refinement steps. */
+
+    /* Set right hand side to one. */
+    for (i = 0; i < n; i++) {
+      b[i] = 1;
+    }
+
+    F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
+                       &n, a, ia, ja, &idum, &nrhs,
+                       iparm, &msglvl, b, x, &error);
+
+    if (error != 0) {
+      printf("\nERROR during solution: %d", error);
+      exit(3);
+    }
+
+    printf("\nSolve completed ... ");
+    printf("\nThe solution of the system is: ");
+    for (i = 0; i < n; i++) {
+      printf("\n x [%d] = % f", i, x[i] );
+    }
+    printf ("\n\n");
+
+    /* -------------------------------------------------------------------- */
+    /* ..  Termination and release of memory.                               */
+    /* -------------------------------------------------------------------- */
+    phase = -1;                 /* Release internal memory. */
+
+    F77_FUNC(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
+                       &n, &ddum, ia, ja, &idum, &nrhs,
+                       iparm, &msglvl, &ddum, &ddum, &error);
+
+    typedef Dune::FieldMatrix<double,1,1> M;
+    Dune::BCRSMatrix<M> B(8,8,Dune::BCRSMatrix<M>::random);
+
+    // initially set row size for each row
+    B.setrowsize(0,4);
+    B.setrowsize(1,3);
+    B.setrowsize(2,2);
+    B.setrowsize(3,2);
+    B.setrowsize(4,1);
+    B.setrowsize(5,3);
+    B.setrowsize(6,2);
+    B.setrowsize(7,3);
+
+    // finalize row setup phase
+    B.endrowsizes();
+
+    // add column entries to rows
+    B.addindex(0,0); B.addindex(0,2); B.addindex(0,5); B.addindex(0,6);
+    B.addindex(1,1); B.addindex(1,2); B.addindex(1,4);
+    B.addindex(2,2); B.addindex(2,7);
+    B.addindex(3,3); B.addindex(3,6);
+    B.addindex(4,4);
+    B.addindex(5,2); B.addindex(5,5); B.addindex(5,7);
+    B.addindex(6,1); B.addindex(6,6);
+    B.addindex(7,2); B.addindex(7,6); B.addindex(7,7);
+
+    // finalize column setup phase
+    B.endindices();
+
+    // set entries using the random access operator
+    B[0][0] = 7; B[0][2] = 1; B[0][5] = 2; B[0][6] = 7;
+    B[1][1] = -4; B[1][2] = 8; B[1][4] = 2;
+    B[2][2] = 1; B[2][7] = 5;
+    B[3][3] = 7; B[3][6] = 9;
+    B[4][4] = -4;
+    B[5][2] = 7; B[5][5] = 3; B[5][7] = 8;
+    B[6][1] = 1; B[6][6] = 11;
+    B[7][2] = -3; B[7][6] = 2; B[7][7] = 5;
+
+    //printmatrix(std::cout, B, "matrix B", "row", 9, 1);
+
+    typedef Dune::FieldVector<double, 1> VB;
+    typedef Dune::BlockVector<VB> Vector;
+    typedef Dune::BCRSMatrix<M> Matrix;
+    Dune::MatrixAdapter<Matrix,Vector,Vector> op(B);        // make linear operator from A
+    Dune::SeqPardiso<Matrix,Vector,Vector> pardiso(B);        // preconditioner object
+    Dune::LoopSolver<Vector> loop(op, pardiso, 1E-14, 2, 1);         // an inverse operator
+
+    Vector f(n);
+    f = 1;
+    Vector y(n);
+    y = 0;
+    Dune::InverseOperatorResult r;
+    loop.apply(y, f, r);
+
+    std::cout << "\nSolve completed ... ";
+    std::cout << "\nThe solution of the system is: ";
+    for (i = 0; i < n; i++) {
+      std::cout << "\n x [" << i << "] = " << y[i];
+    }
+    std::cout << "\n";
+
+#else
+    DUNE_THROW(Dune::NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
+#endif
+
+
+    return 0;
+  }
+  catch (Dune::Exception &e) {
+    std::cerr << "Dune reported error: " << e << std::endl;
+  }
+  catch (...) {
+    std::cerr << "Unknown exception thrown!" << std::endl;
+  }
+}