From 6c39b1711de0674b5dbe298391d9ba6c3cc6c632 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 20 Dec 2011 12:56:21 +0000
Subject: [PATCH] send Pardiso preconditioner to devel
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7166 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
CMakeLists.txt | 3 +-
configure.ac | 1 -
dumux/common/pardiso.hh | 456 ----------------------------
test/common/CMakeLists.txt | 2 +-
test/common/Makefile.am | 2 +-
test/common/pardiso/CMakeLists.txt | 8 -
test/common/pardiso/Makefile.am | 8 -
test/common/pardiso/test_pardiso.cc | 230 --------------
8 files changed, 3 insertions(+), 707 deletions(-)
delete mode 100644 dumux/common/pardiso.hh
delete mode 100644 test/common/pardiso/CMakeLists.txt
delete mode 100644 test/common/pardiso/Makefile.am
delete mode 100644 test/common/pardiso/test_pardiso.cc
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1882354dac..e80fac121f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,7 +2,7 @@
# general stuff
cmake_minimum_required(VERSION 2.6)
set(ProjectName "DuMuX")
-set(ProjectVersion "2.1svn")
+set(ProjectVersion "2.1-svn")
set(ProjectMaintainer "Bernd Flemisch")
set(ProjectMaintainerEmail "Bernd.Flemisch_at_iws dot uni-stuttgart dot de")
project(${ProjectName} CXX)
@@ -219,6 +219,5 @@ add_test(test_dec2p2c util/runTest.sh test/decoupled/2p2c/test_dec2p2c-reference
add_test(test_fluidsystems test/material/fluidsystems/test_fluidsystems)
add_test(test_fluidsystems test/material/ncpflash/test_ncpflash)
add_test(tutorial_coupled tutorial/tutorial_coupled 1 1)
-add_test(tutorial_coupled tutorial/tutorial_coupled 1 1)
add_test(tutorial_decoupled tutorial/tutorial_decoupled 1)
diff --git a/configure.ac b/configure.ac
index 426f9914f1..1fffbb86d6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -66,7 +66,6 @@ AC_CONFIG_FILES([dumux.pc
test/boxmodels/richards/Makefile
test/common/Makefile
test/common/generalproblem/Makefile
- test/common/pardiso/Makefile
test/common/propertysystem/Makefile
test/common/spline/Makefile
test/decoupled/Makefile
diff --git a/dumux/common/pardiso.hh b/dumux/common/pardiso.hh
deleted file mode 100644
index b4a4f4c0db..0000000000
--- a/dumux/common/pardiso.hh
+++ /dev/null
@@ -1,456 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2007-2009 by Bernd Flemisch *
- * Institute of Hydraulic Engineering *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Sequential preconditioner using the Pardiso direct solver.
- */
-#ifndef DUMUX_PARDISO_HH
-#define DUMUX_PARDISO_HH
-
-#include <dune/istl/preconditioners.hh>
-
-/* Change this, if your Fortran compiler does not append underscores. */
-/* e.g. the AIX compiler: #define F77_FUN(func) func */
-
-#ifdef HAVE_PARDISO
-
-#ifdef AIX
-#define F77_FUN(func) func
-#else
-#define F77_FUN(func) func ## _
-#endif
-
-/* PARDISO prototype. */
-
-extern "C" int F77_FUN(pardisoinit)
-(void *pt, int *mtype, int *solver, int *iparm, double *dparm, int *error);
-
-extern "C" int F77_FUN(pardiso)
-(void *pt, int *maxfct, int *mnum, int *mtype, int *phase, int *n,
- double *a, int *ia, int *ja, int *perm, int *nrhs, int *iparm,
- int *msglvl, double *b, double *x, int *error, double *dparm);
-
-#endif /* HAVE_PARDISO */
-
-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 Dune::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
- };
-
-
- SeqPardiso(bool verbose = false, int n = 0, double w = 0.0)
- : verbose_(verbose)
- {
-
-#ifdef HAVE_PARDISO
- mtype_ = 11;
- nrhs_ = 1;
- num_procs_ = 1;
- maxfct_ = 1;
- mnum_ = 1;
- msglvl_ = 0;
- error_ = 0;
-
- solver_ = 0; // solver_ = 0, choose sparse direct solver, = 1 multi-recursive iterative solver
-#else
- DUNE_THROW(Dune::NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
-#endif
- }
-
- void factorize (M& A)
- {
-#ifdef HAVE_PARDISO
-
- RowIterator i0 = A.begin();
- ColIterator j0 = (*i0).begin();
-
- systemsize_ = (*j0).N();
- n_ = A.N()*systemsize_;
- int nnz = 0;
- RowIterator endi = A.end();
- int rows = 0;
- for (RowIterator i = A.begin(); i != endi; ++i)
- {
- rows++;
- ColIterator endj = (*i).end();
- for (ColIterator j = (*i).begin(); j != endj; ++j) {
- nnz += systemsize_*systemsize_;
- }
- }
- //std::cout << "rows = " << rows;
- if (verbose_)
- std::cout << "SeqPardiso: 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)
- {
- for (int iComp = 0; iComp < systemsize_; iComp++) {
- ia_[i.index()*systemsize_ + iComp] = count+1;
- ColIterator endj = (*i).end();
- for (ColIterator j = (*i).begin(); j != endj; ++j) {
- for (int jComp = 0; jComp < systemsize_; jComp++) {
- a_[count] = (*j)[iComp][jComp];
- ja_[count] = j.index()*systemsize_ + jComp + 1;
-
- count++;
- }
- }
- }
- }
- ia_[n_] = count+1;
-
- F77_FUN(pardisoinit) (pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
-
- if (error_)
- {
- switch(error_)
- {
- case -10:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. No license file found. Error code " << error_);
- break;
- case -11:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. License has expired. Error code " << error_);
- break;
- case -12:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. Wrong username or hostname. Error code " << error_);
- break;
- default:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. Error code " << error_);
- break;
- }
- }
-
-
- phase_ = 11;
- int idum;
- double ddum;
- iparm_[2] = num_procs_;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, a_, ia_, ja_, &idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
-
- if (error_ != 0)
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: Reordering failed. Error code " << error_);
-
- if (verbose_)
- std::cout << " Reordering completed. Number of nonzeros in factors = " << iparm_[17] << std::endl;
-
- phase_ = 22;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, a_, ia_, ja_, idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
-
- if (error_ != 0)
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
-
- if (verbose_)
- std::cout << " Factorization completed." << std::endl;
-#endif
- }
-
- void initAndFactor(M& A)
- {
-#ifdef HAVE_PARDISO
-
- mtype_ = 11;
- nrhs_ = 1;
- num_procs_ = 1;
- maxfct_ = 1;
- mnum_ = 1;
- msglvl_ = 0;
- error_ = 0;
- solver_ = 0; // solver_ = 0, choose sparse direct solver, = 1 multi-recursive iterative solver
-
- RowIterator i0 = A.begin();
- ColIterator j0 = (*i0).begin();
-
- systemsize_ = (*j0).N();
- n_ = A.N()*systemsize_;
-
-
-
- int nnz = 0;
- RowIterator endi = A.end();
- int rows = 0;
- for (RowIterator i = A.begin(); i != endi; ++i)
- {
- rows++;
- ColIterator endj = (*i).end();
- for (ColIterator j = (*i).begin(); j != endj; ++j) {
- nnz += systemsize_*systemsize_;
- }
- }
- if (verbose_)
- std::cout << "SeqPardiso: 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)
- {
- for (int iComp = 0; iComp < systemsize_; iComp++) {
- ia_[i.index()*systemsize_ + iComp] = count+1;
- ColIterator endj = (*i).end();
- for (ColIterator j = (*i).begin(); j != endj; ++j) {
- for (int jComp = 0; jComp < systemsize_; jComp++) {
- a_[count] = (*j)[iComp][jComp];
- ja_[count] = j.index()*systemsize_ + jComp + 1;
-
- count++;
- }
- }
- }
- }
- ia_[n_] = count+1;
-
- F77_FUN(pardisoinit) (pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
-
- if (error_)
- {
- switch(error_)
- {
- case -10:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. No license file found. Error code " << error_);
- break;
- case -11:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. License has expired. Error code " << error_);
- break;
- case -12:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. Wrong username or hostname. Error code " << error_);
- break;
- default:
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: pardisoinit failed. Error code " << error_);
- break;
- }
- }
-
- phase_ = 11;
- int idum;
- double ddum;
- iparm_[2] = num_procs_;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, a_, ia_, ja_, &idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
-
- if (error_ != 0)
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: Reordering failed. Error code " << error_);
-
- if (verbose_)
- std::cout << " Reordering completed. Number of nonzeros in factors = " << iparm_[17] << std::endl;
-
- phase_ = 22;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, a_, ia_, ja_, &idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
-
- if (error_ != 0)
- DUNE_THROW(Dune::MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
-
- if (verbose_)
- std::cout << " Factorization completed." << std::endl;
-
-#else
- DUNE_THROW(Dune::NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
-#endif
- }
-
- /*! \brief Constructor.
-
- Constructor gets all parameters to operate the preconditioner.
- \param A The matrix to operate on
- \param verbose Write messages if true
- */
- SeqPardiso (M& A, bool verbose = false)
- : verbose_(verbose)
- {
- initAndFactor(A);
- }
-
- SeqPardiso (M& A, int iter, double relaxation)
- : verbose_(false)
- {
- initAndFactor(A);
- }
-
- /*! \brief Prepare the preconditioner.
-
- A solver solves a linear operator equation A(x)=b by applying
- one or several steps of the preconditioner. The method pre()
- is called before the first apply operation.
- b and x are right hand side and solution vector of the linear
- system respectively. It may. e.g., scale the system, allocate memory or
- compute a (I)LU decomposition.
- Note: The ILU decomposition could also be computed in the constructor
- or with a separate method of the derived method if several
- linear systems with the same matrix are to be solved.
-
- \param x The left hand side of the equation.
- \param b The right hand side of the equation.
- */
- virtual void pre (X& x, Y& b) {}
-
- /*! \brief Apply one step of the preconditioner to the system \f$ A(v)=d \f$.
-
- On entry \f$ v=0 \f$ and \f$ d=b-A(x) \f$ (although this might not be
- computed in that way. On exit v contains the update, i.e
- one step computes \f$ v = M^{-1} d \f$ where \f$ M \f$ is the
- approximate inverse of the operator \f$ A \f$ characterizing
- the preconditioner.
- \param[out] v The update to be computed
- \param d The current defect.
- */
- virtual void apply (X& v, const Y& d)
- {
-#ifdef HAVE_PARDISO
- phase_ = 33;
-
- iparm_[7] = 1; /* Max numbers of iterative refinement steps. */
- int idum;
- double x[2*n_];
- double b[2*n_];
- for (typename X::size_type i = 0; i < v.size(); i++) {
- for (int comp = 0; comp < systemsize_; comp++) {
- x[i*systemsize_ + comp] = v[i][comp];
- b[i*systemsize_ + comp] = d[i][comp];
- }
- }
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, a_, ia_, ja_, &idum, &nrhs_,
- //&n_, &ddum, &idum, &idum, &idum, &nrhs_,
- iparm_, &msglvl_, b, x, &error_, dparm_);
-
- if (error_ != 0)
- DUNE_THROW(Dune::MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
-
- for (typename X::size_type i = 0; i < v.size(); i++)
- for (int comp = 0; comp < systemsize_; comp++)
- v[i][comp] = x[i*systemsize_ + comp];
-#endif
- }
-
- /*! \brief Clean up.
-
- This method is called after the last apply call for the
- linear system to be solved. Memory may be deallocated safely
- here. x is the solution of the linear equation.
-
- \param x The right hand side of the equation.
- */
- virtual void post (X& x)
- {
-#ifdef HAVE_PARDISO
- phase_ = -1; // Release internal memory.
- int idum;
- double ddum;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, &ddum, ia_, ja_, &idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
- delete a_;
- delete ia_;
- delete ja_;
-#endif
- }
-
- ~SeqPardiso()
- {
-#ifdef HAVE_PARDISO
- if (phase_ != -1) {
- phase_ = -1; // Release internal memory.
- int idum;
- double ddum;
-
- F77_FUN(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase_,
- &n_, &ddum, ia_, ja_, &idum, &nrhs_,
- iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
- delete a_;
- delete ia_;
- delete ja_;
- }
-#endif
- }
-
-private:
- 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.
- int systemsize_;
- int phase_;
- bool verbose_;
- double dparm_[64];
- int solver_;
-};
-
-}
-
-
-
-
-
-
-#endif
-
diff --git a/test/common/CMakeLists.txt b/test/common/CMakeLists.txt
index 9459d3c589..733a9f761e 100644
--- a/test/common/CMakeLists.txt
+++ b/test/common/CMakeLists.txt
@@ -1,7 +1,7 @@
# directories which are commented out do not compile at the moment
# (they also don't compile if using the old build system)!
-add_subdirectory("pardiso")
+add_subdirectory("generalproblem")
add_subdirectory("propertysystem")
add_subdirectory("spline")
diff --git a/test/common/Makefile.am b/test/common/Makefile.am
index 2f543318c6..e90a6e2860 100644
--- a/test/common/Makefile.am
+++ b/test/common/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = generalproblem pardiso propertysystem spline .
+SUBDIRS = generalproblem propertysystem spline .
EXTRA_DIST = CMakeLists.txt
diff --git a/test/common/pardiso/CMakeLists.txt b/test/common/pardiso/CMakeLists.txt
deleted file mode 100644
index bc323e3c90..0000000000
--- a/test/common/pardiso/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-# add build targets
-ADD_EXECUTABLE("test_pardiso" test_pardiso.cc)
-TARGET_LINK_LIBRARIES("test_pardiso" ${DumuxLinkLibraries})
-
-# add required libraries and includes to the build flags
-LINK_DIRECTORIES(${DumuxLinkDirectories})
-INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
-
diff --git a/test/common/pardiso/Makefile.am b/test/common/pardiso/Makefile.am
deleted file mode 100644
index 2f3eee0c04..0000000000
--- a/test/common/pardiso/Makefile.am
+++ /dev/null
@@ -1,8 +0,0 @@
-# tests where program to build and program to run are equal
-check_PROGRAMS = test_pardiso
-
-EXTRA_DIST = CMakeLists.txt
-
-test_pardiso_SOURCES = test_pardiso.cc
-
-include $(top_srcdir)/am/global-rules
diff --git a/test/common/pardiso/test_pardiso.cc b/test/common/pardiso/test_pardiso.cc
deleted file mode 100644
index 1fef48a84c..0000000000
--- a/test/common/pardiso/test_pardiso.cc
+++ /dev/null
@@ -1,230 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2008 by Bernd Flemisch *
- * Institute of Hydraulic Engineering *
- * University of Stuttgart, Germany *
- * email: andreas.lauser _at_ iws.uni-stuttgart.de *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief test for the Pardiso preconditioner
- */
-#include "config.h"
-#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 "dumux/common/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];
- double dparm[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. */
- int solver = 0;
-
- /* -------------------------------------------------------------------- */
- /* .. Setup Pardiso control parameters. */
- /* -------------------------------------------------------------------- */
-
- F77_FUN(pardisoinit) (pt, &mtype, &solver, iparm, dparm, &error);
- if (error != 0) {
- printf("\nERROR during initialization: %d", error);
- exit(3);
- }
-
- 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_FUN(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
- &n, a, ia, ja, &idum, &nrhs,
- iparm, &msglvl, b, x, &error, dparm);
-
- 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_FUN(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
- &n, &ddum, ia, ja, &idum, &nrhs,
- iparm, &msglvl, &ddum, &ddum, &error, dparm);
-
- 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;
- }
-}
--
GitLab