Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jakub.both/dune-istl
  • eduardo.bueno/dune-istl
  • tkoch/dune-istl
  • pipping/dune-istl
  • andreas.brostrom/dune-istl
  • stephan.hilb/dune-istl
  • max.kahnt/dune-istl
  • patrick.jaap/dune-istl
  • tobias.meyer.andersen/dune-istl
  • andreas.thune/dune-istl
  • dominic/dune-istl
  • ansgar/dune-istl
  • exadune/dune-istl
  • lars.lubkoll/dune-istl
  • govind.sahai/dune-istl
  • maikel.nadolski/dune-istl
  • markus.blatt/dune-istl
  • core/dune-istl
  • lisa_julia.nebel/dune-istl
  • michael.sghaier/dune-istl
  • liesel.schumacher/dune-istl
  • Xinyun.Li/dune-istl
  • lukas.renelt/dune-istl
  • simon.praetorius/dune-istl
  • rene.milk/dune-istl
  • lasse.hinrichsen/dune-istl
  • nils.dreier/dune-istl
  • claus-justus.heine/dune-istl
  • felix.mueller/dune-istl
  • kilian.weishaupt/dune-istl
  • lorenzo.cerrone/dune-istl
  • jakob.torben/dune-istl
  • liam.keegan/dune-istl
  • alexander.mueller/dune-istl
34 results
Show changes
Commits on Source (26)
Hide whitespace changes
Inline Side-by-side
Showing
with 879 additions and 514 deletions
LICENSE.md
View file @ 77390949
Copyright holders:
2014--2016 Marco Agnese
2003--2013 Peter Bastian
2004--2014 Markus Blatt
2014 Ansgar Burchardt
2004--2016 Markus Blatt
2014--2016 Ansgar Burchardt
2004--2005 Adrian Burri
2008--2014 Andreas Dedner
2003--2014 Christian Engwer
2005--2012 Jorrit Fahlke
2008--2013 Bernd Flemisch
2005--2013 Carsten Gräser
2012--2014 Christoph Grüninger
2008--2015 Andreas Dedner
2003--2016 Christian Engwer
2005--2016 Jorrit Fahlke
2008--2015 Bernd Flemisch
2005--2016 Carsten Gräser
2015--2016 Felix Gruber
2012--2016 Christoph Grüninger
2016 René Heß
2016 Stephan Hilb
2012--2013 Olaf Ippisch
2013--2014 Dominic Kempf
2004--2013 Robert Klöfkorn
2013--2014 Arne Morten Kvarving (SINTEF)
2013--2016 Dominic Kempf
2015 Emmanouil Kiagias
2004--2016 Robert Klöfkorn
2016 Timo Koch
2007 Sreejith Pulloor Kuttanikkad
2013--2016 Arne Morten Kvarving (SINTEF)
2009--2013 Andreas Lauser
2012--2014 Tobias Malkmus
2012--2016 Tobias Malkmus
2007--2009 Sven Marnach
2013 Rene Milk
2013--2014 Steffen Müthing
2013 René Milk
2013--2016 Steffen Müthing
2016 Maikel Nadolski
2003--2005 Thimo Neubauer
2010--2012 Rebecca Neumann
2012--2014 Andreas Nüßing
2008--2014 Martin Nolte
2014 Marian Piatkowski
2011--2013 Elias Pipping
2012--2015 Andreas Nüßing
2013--2014 Marian Piatkowski
2011--2016 Elias Pipping
2013 Jurgis Pods
2007 Sreejith Pulloor Kuttanikkad
2009 Atgeirr Rasmussen
2008--2013 Uli Sack
2004--2014 Oliver Sander
2004--2016 Oliver Sander
2013 Bård Skaflestad
2006--2010 Martin Weiser
2015 Sebastian Westerheide
2011--2012 Matthias Wohlmuth
2014 Jonathan Youett
2014--2016 Jonathan Youett
The DUNE library and headers are licensed under version 2 of the GNU
General Public License (see below), with a special exception for
......
cmake/modules/FindSuperLU.cmake
View file @ 77390949
......@@ -125,8 +125,20 @@ int main(void)
}"
SUPERLU_MIN_VERSION_5)
if(SUPERLU_MIN_VERSION_4_3)
include(CheckIncludeFiles)
set(HAVE_SLU_DDEFS_H 1)
check_include_files(slu_sdefs.h HAVE_SLU_SDEFS_H)
check_include_files(slu_cdefs.h HAVE_SLU_CDEFS_H)
check_include_files(slu_zdefs.h HAVE_SLU_ZDEFS_H)
endif(SUPERLU_MIN_VERSION_4_3)
cmake_pop_check_state()
set(SUPERLU_INT_TYPE "int" CACHE STRING
"The integer version that SuperLU was compiled for (Default is int.
Should be the same as int_t define in e.g. slu_sdefs.h")
if(NOT SUPERLU_MIN_VERSION_4)
set(SUPERLU_WITH_VERSION "SuperLU < 4.0" CACHE STRING
"Human readable string containing SuperLU version information.")
......
config.h.cmake
View file @ 77390949
......@@ -31,6 +31,22 @@
/* Define to ENABLE_SUPERLU if the SuperLU library is available */
#cmakedefine HAVE_SUPERLU ENABLE_SUPERLU
/* Define to the integer type that SuperLU was compiled for
See e.g. what int_t is defined to in slu_sdefs.h */
#cmakedefine SUPERLU_INT_TYPE @SUPERLU_INT_TYPE@
/* Define to 1 if header slu_sdefs.h is there. */
#cmakedefine01 HAVE_SLU_SDEFS_H
/* Define to 1 if header slu_ddefs.h is there. */
#cmakedefine01 HAVE_SLU_DDEFS_H
/* Define to 1 if header slu_cdefs.h is there. */
#cmakedefine01 HAVE_SLU_CDEFS_H
/* Define to 1 if header slu_zdefs.h is there. */
#cmakedefine01 HAVE_SLU_ZDEFS_H
/* Define to ENABLE_ARPACKPP if the ARPACK++ library is available */
#cmakedefine HAVE_ARPACKPP ENABLE_ARPACKPP
......
doc/doxygen/modules.txt
View file @ 77390949
/**
@defgroup ISTL Iterative Solvers Template Library (ISTL)
@brief Iterative Solvers supporting block recursive matrix and
vector classes at compile time.
The Iterative Solver Template Library applies generic programming
in C++ to the
The Iterative Solver Template Library applies generic programming
in C++ to the
domain of iterative solvers of linear systems stemming
from finite element discretizations. Those discretizations exhibit a
lot of structure, e.g:
......@@ -13,7 +13,7 @@
<ol>
<li>Certain discretizations for systems of PDEs or higher order
methods result in matrices where individual entries are replaced by
small blocks, say of size \f$2\times 2\f$ or \f$4\times 4\f$.
small blocks, say of size \f$2\times 2\f$ or \f$4\times 4\f$.
Dense blocks of different sizes
e.g. arise in \f$hp\f$ Discontinuous Galerkin discretization methods.
It is straightforward and
......@@ -30,13 +30,22 @@
Our matrix
and vector interface supports a block recursive structure. Each
sparse matrix entry can itself be either a sparse or a small
sparse matrix entry can itself be either a sparse or a small
dense matrix.
The solvers use this recursive block structure via template meta
The solvers use this recursive block structure via template meta
programming at compile time.
ISTL consists of the \ref ISTL_SPMV "matrix and vector API" and
ISTL consists of the \ref ISTL_SPMV "matrix and vector API" and
the \ref ISTL_Solvers "solvers" which use the \ref ISTL_Prec preconditioners.
*/
/** @defgroup ISTL_Solvers Iterative Solvers
@ingroup ISTL
*/
/** @defgroup ISTL_Eigenvalue Eigenvalue Solvers
@ingroup ISTL_Solvers
*/
/** @defgroup ISTL_Comm Communication Interface
@ingroup ISTL
*/
doc/istl.tex
View file @ 77390949
......@@ -83,8 +83,8 @@ there are many libraries available on the internet for doing sparse matrix/vecto
computations. A comprehensive overview is given in \cite{LALinks}.
The
widely availably Basic Linear Algebra Subprograms (BLAS) standard has
been extended to cover als sparse matrices \cite{BLASTForum}. BLAS
widely available Basic Linear Algebra Subprograms (BLAS) standard has
been extended to cover as sparse matrices \cite{BLASTForum}. BLAS
divides the available functions into level 1 (vector operations),
level 2 (vector/matrix operations) and level 3 (matrix/matrix
operations). BLAS for sparse matrices contains only level 1 and 2
......@@ -94,7 +94,7 @@ only a FORTRAN and C interface. As a consequence, the interface is
``coarse grained'', meaning that ``small'' functions such as access to
individual matrix elements is relatively slow.
Generic programming techniqes in C++ offer the possibility to combine
Generic programming techniques in C++ offer the possibility to combine
flexibility and reuse (``efficiency of the programmer'') with fast
execution (``efficieny of the program'') as has been demonstrated with
the Standard Template Library (STL), \cite{Stroustrup} or the Blitz++
......@@ -174,7 +174,7 @@ In mathemetics vectors are elements of a vector space. A vector space
$V(\K)$, defined over a field $\K$, is a set of elements with two
operations: (i) vector space addition $+ : V\times V \to V$ and (ii) scalar
multiplication $* : \K\times V \to V$. These operations obey certain formal
rules, see your favourite textbook on linear algebra,
rules, see your favorite textbook on linear algebra,
e.~g. \cite{LaBook}. In addition a
vector space may be normed, i.~e.~there is a function (obeying certain
rules) $\|.\| : V \to \R$ which measures distance in the vector
......@@ -187,7 +187,7 @@ field, such as $\K=\R$ or $\K=\C$ and take a tensor product:
V = \K^n = \underbrace{\K\times\K\times\ldots\times\K}_{\text{$n$ times}}.
\end{equation*}
$n\in\N$ is called the dimension of the vector space. There are also
infinite-dimensional vector spaceswhich are, however, not of interest
infinite-dimensional vector spaces which are, however, not of interest
in the context here. The idea of tensor products can be generalized.
If we have vector spaces $V_1(\K),\ldots,V_n(\K)$ we can construct a
new vector space by setting
......@@ -628,8 +628,8 @@ The base class
linear maps. The
template parameter \lstinline!X! is the type of the domain and
\lstinline!Y! is the type of the range of the operator. A linear
operator provides the methods \lstinline!apply(const X& x, Y& y)! and
apply \lstinline!applyscaledadd(const X& x, Y& y)! performing the
operator provides the methods \lstinline!apply(const X& x, Y& y)! and
\lstinline!applyscaledadd(const X& x, Y& y)! performing the
operations $y = A(x)$ and $y = y + \alpha A(x)$, respectively.
The subclass
\lstinline!template<class M, class X, class Y> AssembledLinearOperator!
......@@ -726,7 +726,7 @@ providing them with the vector implementation used.
\hline
\textbf{class}&\textbf{implements}\\\hline\hline
\lstinline!LoopSolver!& only apply precoditioner multiple time\\
\lstinline!GradientSolver!& preconditioned radient method\\
\lstinline!GradientSolver!& preconditioned gradient method\\
\lstinline!CGSolver!&preconditioned conjugate gradient method\\
\lstinline!BiCGStab!&preconditioned biconjugate gradient stabilized method\\\hline
\end{tabular}
......@@ -812,10 +812,10 @@ operation ($x = y + \alpha z$) in Tables
\hline
\end{tabular}}}
\end{table}
The code was comiled with the GNU C++
The code was compiled with the GNU C++
compiler version 4.0 with -O3 optimization. In the tables $N$ is the
number of
unknown blocks (equals the number of unknows for the scalar cases in
unknown blocks (equals the number of unknowns for the scalar cases in
Tables \ref{tab:perf_sp}, \ref{tab:perf_daxpy}, \ref{tab:perf_gs}).
The performance for the scalarproduct,
see Table \ref{tab:perf_sp},
......
dune.module
View file @ 77390949
Module: dune-istl
Version: 3.0-git
Version: 2.6-git
Maintainer: dune-devel@dune-project.org
Depends: dune-common (>= 3.0)
Depends: dune-common (>= 2.6)
Whitespace-Hook: Yes
dune/istl/CMakeLists.txt
View file @ 77390949
......@@ -42,6 +42,7 @@ install(FILES
solvertype.hh
spqr.hh
superlu.hh
superlufunctions.hh
supermatrix.hh
umfpack.hh
vbvector.hh
......
dune/istl/eigenvalue/arpackpp.hh
View file @ 77390949
......@@ -3,7 +3,7 @@
#ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
#if HAVE_ARPACKPP
#if HAVE_ARPACKPP || defined DOXYGEN
#include <cmath> // provides std::abs, std::pow, std::sqrt
......@@ -26,176 +26,182 @@
#endif
#include "arssym.h" // provides ARSymStdEig
namespace Dune
{
/**
* \brief Wrapper for a DUNE-ISTL BCRSMatrix which can be used
* together with those algorithms of the ARPACK++ library
* which solely perform the products A*v and/or A^T*A*v
* and/or A*A^T*v.
*
* \todo The current implementation is limited to DUNE-ISTL
* BCRSMatrix types with blocklevel 2. An extension to
* blocklevel >= 2 might be provided in a future version.
*
* \tparam BCRSMatrix Type of a DUNE-ISTL BCRSMatrix;
* is assumed to have blocklevel 2.
*
* \author Sebastian Westerheide.
*/
template <class BCRSMatrix>
class ArPackPlusPlus_BCRSMatrixWrapper
{
public:
//! Type of the underlying field of the matrix
typedef typename BCRSMatrix::field_type Real;
public:
//! Construct from BCRSMatrix A
ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
: A_(A),
m_(A_.M() * mBlock), n_(A_.N() * nBlock)
{
// assert that BCRSMatrix type has blocklevel 2
static_assert
(BCRSMatrix::blocklevel == 2,
"Only BCRSMatrices with blocklevel 2 are supported.");
// allocate memory for auxiliary block vector objects
// which are compatible to matrix rows / columns
domainBlockVector.resize(A_.N(),false);
rangeBlockVector.resize(A_.M(),false);
}
//! Perform matrix-vector product w = A*v
inline void multMv (Real* v, Real* w)
{
// get vector v as an object of appropriate type
arrayToDomainBlockVector(v,domainBlockVector);
/** @addtogroup ISTL_Eigenvalue
@{
*/
// perform matrix-vector product
A_.mv(domainBlockVector,rangeBlockVector);
// get vector w from object of appropriate type
rangeBlockVectorToArray(rangeBlockVector,w);
};
//! Perform matrix-vector product w = A^T*A*v
inline void multMtMv (Real* v, Real* w)
namespace Impl {
/**
* Wrapper for a DUNE-ISTL BCRSMatrix which can be used
* together with those algorithms of the ARPACK++ library
* which solely perform the products A*v and/or A^T*A*v
* and/or A*A^T*v.
*
* \todo The current implementation is limited to DUNE-ISTL
* BCRSMatrix types with blocklevel 2. An extension to
* blocklevel >= 2 might be provided in a future version.
*
* \tparam BCRSMatrix Type of a DUNE-ISTL BCRSMatrix;
* is assumed to have blocklevel 2.
*
* \author Sebastian Westerheide.
*/
template <class BCRSMatrix>
class ArPackPlusPlus_BCRSMatrixWrapper
{
// get vector v as an object of appropriate type
arrayToDomainBlockVector(v,domainBlockVector);
// perform matrix-vector product
A_.mv(domainBlockVector,rangeBlockVector);
A_.mtv(rangeBlockVector,domainBlockVector);
public:
//! Type of the underlying field of the matrix
typedef typename BCRSMatrix::field_type Real;
public:
//! Construct from BCRSMatrix A
ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
: A_(A),
m_(A_.M() * mBlock), n_(A_.N() * nBlock)
{
// assert that BCRSMatrix type has blocklevel 2
static_assert
(BCRSMatrix::blocklevel == 2,
"Only BCRSMatrices with blocklevel 2 are supported.");
// allocate memory for auxiliary block vector objects
// which are compatible to matrix rows / columns
domainBlockVector.resize(A_.N(),false);
rangeBlockVector.resize(A_.M(),false);
}
// get vector w from object of appropriate type
domainBlockVectorToArray(domainBlockVector,w);
};
//! Perform matrix-vector product w = A*v
inline void multMv (Real* v, Real* w)
{
// get vector v as an object of appropriate type
arrayToDomainBlockVector(v,domainBlockVector);
//! Perform matrix-vector product w = A*A^T*v
inline void multMMtv (Real* v, Real* w)
{
// get vector v as an object of appropriate type
arrayToRangeBlockVector(v,rangeBlockVector);
// perform matrix-vector product
A_.mv(domainBlockVector,rangeBlockVector);
// perform matrix-vector product
A_.mtv(rangeBlockVector,domainBlockVector);
A_.mv(domainBlockVector,rangeBlockVector);
// get vector w from object of appropriate type
rangeBlockVectorToArray(rangeBlockVector,w);
};
// get vector w from object of appropriate type
rangeBlockVectorToArray(rangeBlockVector,w);
};
//! Perform matrix-vector product w = A^T*A*v
inline void multMtMv (Real* v, Real* w)
{
// get vector v as an object of appropriate type
arrayToDomainBlockVector(v,domainBlockVector);
//! Return number of rows in the matrix
inline int nrows () const { return m_; }
// perform matrix-vector product
A_.mv(domainBlockVector,rangeBlockVector);
A_.mtv(rangeBlockVector,domainBlockVector);
//! Return number of columns in the matrix
inline int ncols () const { return n_; }
// get vector w from object of appropriate type
domainBlockVectorToArray(domainBlockVector,w);
};
protected:
// Number of rows and columns in each block of the matrix
constexpr static int mBlock = BCRSMatrix::block_type::rows;
constexpr static int nBlock = BCRSMatrix::block_type::cols;
// Type of vectors in the domain of the linear map associated with
// the matrix, i.e. block vectors compatible to matrix rows
constexpr static int dbvBlockSize = nBlock;
typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
// Type of vectors in the range of the linear map associated with
// the matrix, i.e. block vectors compatible to matrix columns
constexpr static int rbvBlockSize = mBlock;
typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
// Types for vector index access
typedef typename DomainBlockVector::size_type dbv_size_type;
typedef typename RangeBlockVector::size_type rbv_size_type;
typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
// Get vector v from a block vector object which is compatible to
// matrix rows
static inline void
//! Perform matrix-vector product w = A*A^T*v
inline void multMMtv (Real* v, Real* w)
{
// get vector v as an object of appropriate type
arrayToRangeBlockVector(v,rangeBlockVector);
// perform matrix-vector product
A_.mtv(rangeBlockVector,domainBlockVector);
A_.mv(domainBlockVector,rangeBlockVector);
// get vector w from object of appropriate type
rangeBlockVectorToArray(rangeBlockVector,w);
};
//! Return number of rows in the matrix
inline int nrows () const { return m_; }
//! Return number of columns in the matrix
inline int ncols () const { return n_; }
protected:
// Number of rows and columns in each block of the matrix
constexpr static int mBlock = BCRSMatrix::block_type::rows;
constexpr static int nBlock = BCRSMatrix::block_type::cols;
// Type of vectors in the domain of the linear map associated with
// the matrix, i.e. block vectors compatible to matrix rows
constexpr static int dbvBlockSize = nBlock;
typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
// Type of vectors in the range of the linear map associated with
// the matrix, i.e. block vectors compatible to matrix columns
constexpr static int rbvBlockSize = mBlock;
typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
// Types for vector index access
typedef typename DomainBlockVector::size_type dbv_size_type;
typedef typename RangeBlockVector::size_type rbv_size_type;
typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
// Get vector v from a block vector object which is compatible to
// matrix rows
static inline void
domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
{
for (dbv_size_type block = 0; block < dbv.N(); ++block)
for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
}
{
for (dbv_size_type block = 0; block < dbv.N(); ++block)
for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
}
// Get vector v from a block vector object which is compatible to
// matrix columns
static inline void
// Get vector v from a block vector object which is compatible to
// matrix columns
static inline void
rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
{
for (rbv_size_type block = 0; block < rbv.N(); ++block)
for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
}
public:
//! Get vector v as a block vector object which is compatible to
//! matrix rows
static inline void arrayToDomainBlockVector (const Real* v,
DomainBlockVector& dbv)
{
for (dbv_size_type block = 0; block < dbv.N(); ++block)
for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
}
{
for (rbv_size_type block = 0; block < rbv.N(); ++block)
for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
}
//! Get vector v as a block vector object which is compatible to
//! matrix columns
static inline void arrayToRangeBlockVector (const Real* v,
RangeBlockVector& rbv)
{
for (rbv_size_type block = 0; block < rbv.N(); ++block)
for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
}
public:
//! Get vector v as a block vector object which is compatible to
//! matrix rows
static inline void arrayToDomainBlockVector (const Real* v,
DomainBlockVector& dbv)
{
for (dbv_size_type block = 0; block < dbv.N(); ++block)
for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
}
protected:
// The DUNE-ISTL BCRSMatrix
const BCRSMatrix& A_;
//! Get vector v as a block vector object which is compatible to
//! matrix columns
static inline void arrayToRangeBlockVector (const Real* v,
RangeBlockVector& rbv)
{
for (rbv_size_type block = 0; block < rbv.N(); ++block)
for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
}
// Number of rows and columns in the matrix
const int m_, n_;
protected:
// The DUNE-ISTL BCRSMatrix
const BCRSMatrix& A_;
// Auxiliary block vector objects which are
// compatible to matrix rows / columns
mutable DomainBlockVector domainBlockVector;
mutable RangeBlockVector rangeBlockVector;
};
// Number of rows and columns in the matrix
const int m_, n_;
// Auxiliary block vector objects which are
// compatible to matrix rows / columns
mutable DomainBlockVector domainBlockVector;
mutable RangeBlockVector rangeBlockVector;
};
} // end namespace Impl
/**
* \brief A class template for performing some eigenvalue algorithms
* \brief Wrapper to use a range of ARPACK++ eigenvalue solvers
*
* A class template for performing some eigenvalue algorithms
* provided by the ARPACK++ library which is based on the implicitly
* restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
* Arnoldi/Lanczos process with the implicitily shifted QR technique.
......@@ -288,7 +294,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -390,7 +396,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -491,7 +497,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -608,7 +614,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A^T*A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -720,7 +726,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A^T*A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -828,7 +834,7 @@ namespace Dune
// use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
// and to perform the product A^T*A*v (LU decomposition is not used)
typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
WrappedMatrix A(m_);
// get number of rows and columns in A
......@@ -958,8 +964,9 @@ namespace Dune
const std::string blank_;
};
} // namespace Dune
/** @} */
} // namespace Dune
#endif // HAVE_ARPACKPP
......
dune/istl/eigenvalue/poweriteration.hh
View file @ 77390949
......@@ -24,154 +24,108 @@
#include <dune/istl/io.hh> // provides Dune::printvector(...)
#include <dune/istl/solvers.hh> // provides Dune::InverseOperatorResult
namespace Dune
{
/**
* \brief A linear operator scaling vectors by a scalar value.
* The scalar value can be changed as it is given in a
* form decomposed into an immutable and a mutable part.
*
* \author Sebastian Westerheide.
*/
template <class X, class Y = X>
class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
{
public:
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
enum {category = Dune::SolverCategory::sequential};
ScalingLinearOperator (field_type immutable_scaling,
const field_type& mutable_scaling)
: immutable_scaling_(immutable_scaling),
mutable_scaling_(mutable_scaling)
{}
/** @addtogroup ISTL_Eigenvalue
@{
*/
virtual void apply (const X& x, Y& y) const
{
y = x;
y *= immutable_scaling_*mutable_scaling_;
}
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
namespace Impl {
/**
* \brief A linear operator scaling vectors by a scalar value.
* The scalar value can be changed as it is given in a
* form decomposed into an immutable and a mutable part.
*
* \author Sebastian Westerheide.
*/
template <class X, class Y = X>
class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
{
X temp(x);
temp *= immutable_scaling_*mutable_scaling_;
y.axpy(alpha,temp);
}
public:
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
protected:
const field_type immutable_scaling_;
const field_type& mutable_scaling_;
};
enum {category = Dune::SolverCategory::sequential};
ScalingLinearOperator (field_type immutable_scaling,
const field_type& mutable_scaling)
: immutable_scaling_(immutable_scaling),
mutable_scaling_(mutable_scaling)
{}
/**
* \brief A linear operator representing the sum of two linear operators.
*
* \tparam OP1 Type of the first linear operator.
* \tparam OP2 Type of the second linear operator.
*
* \author Sebastian Westerheide.
*/
template <class OP1, class OP2>
class LinearOperatorSum
: public Dune::LinearOperator<typename OP1::domain_type,
typename OP1::range_type>
{
public:
typedef typename OP1::domain_type domain_type;
typedef typename OP1::range_type range_type;
typedef typename domain_type::field_type field_type;
virtual void apply (const X& x, Y& y) const
{
y = x;
y *= immutable_scaling_*mutable_scaling_;
}
enum {category = Dune::SolverCategory::sequential};
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
{
X temp(x);
temp *= immutable_scaling_*mutable_scaling_;
y.axpy(alpha,temp);
}
LinearOperatorSum (const OP1& op1, const OP2& op2)
: op1_(op1), op2_(op2)
{
static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
"Domain type of both operators doesn't match!");
static_assert(std::is_same<typename OP2::range_type,range_type>::value,
"Range type of both operators doesn't match!");
}
protected:
const field_type immutable_scaling_;
const field_type& mutable_scaling_;
};
virtual void apply (const domain_type& x, range_type& y) const
{
op1_.apply(x,y);
op2_.applyscaleadd(1.0,x,y);
}
virtual void applyscaleadd (field_type alpha,
const domain_type& x, range_type& y) const
/**
* \brief A linear operator representing the sum of two linear operators.
*
* \tparam OP1 Type of the first linear operator.
* \tparam OP2 Type of the second linear operator.
*
* \author Sebastian Westerheide.
*/
template <class OP1, class OP2>
class LinearOperatorSum
: public Dune::LinearOperator<typename OP1::domain_type,
typename OP1::range_type>
{
range_type temp(y);
op1_.apply(x,temp);
op2_.applyscaleadd(1.0,x,temp);
y.axpy(alpha,temp);
}
public:
typedef typename OP1::domain_type domain_type;
typedef typename OP1::range_type range_type;
typedef typename domain_type::field_type field_type;
protected:
const OP1& op1_;
const OP2& op2_;
};
enum {category = Dune::SolverCategory::sequential};
LinearOperatorSum (const OP1& op1, const OP2& op2)
: op1_(op1), op2_(op2)
{
static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
"Domain type of both operators doesn't match!");
static_assert(std::is_same<typename OP2::range_type,range_type>::value,
"Range type of both operators doesn't match!");
}
/**
* \brief Helper class for notifying a DUNE-ISTL linear solver about
* a change of the iteration matrix object in a unified way,
* i.e. independent from the solver's type (direct/iterative).
*
* \author Sebastian Westerheide.
*/
template <typename ISTLLinearSolver, typename BCRSMatrix>
class SolverHelper
{
public:
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
{
static const bool is_direct_solver
= Dune::IsDirectSolver<ISTLLinearSolver>::value;
SolverHelper<ISTLLinearSolver,BCRSMatrix>::
Implementation<is_direct_solver>::setMatrix(solver,matrix);
}
protected:
/**
* \brief Implementation that works together with iterative ISTL
* solvers, e.g. Dune::CGSolver or Dune::BiCGSTABSolver.
*/
template <bool is_direct_solver, typename Dummy = void>
struct Implementation
{
static void setMatrix (ISTLLinearSolver&,
const BCRSMatrix&)
{}
};
virtual void apply (const domain_type& x, range_type& y) const
{
op1_.apply(x,y);
op2_.applyscaleadd(1.0,x,y);
}
/**
* \brief Implementation that works together with direct ISTL
* solvers, e.g. Dune::SuperLU or Dune::UMFPack.
*/
template <typename Dummy>
struct Implementation<true,Dummy>
{
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
virtual void applyscaleadd (field_type alpha,
const domain_type& x, range_type& y) const
{
solver.setMatrix(matrix);
range_type temp(y);
op1_.apply(x,temp);
op2_.applyscaleadd(1.0,x,temp);
y.axpy(alpha,temp);
}
};
};
protected:
const OP1& op1_;
const OP2& op2_;
};
} // end namespace Impl
/**
* \brief A class template for performing some iterative eigenvalue algorithms
* based on power iteration.
* \brief Iterative eigenvalue algorithms based on power iteration.
*
* Given a square matrix whose eigenvalues shall be considered, this class
* template provides methods for performing the power iteration algorithm,
......@@ -213,8 +167,8 @@ namespace Dune
// Type definitions for type of iteration operator (m_ - mu_*I)
typedef typename Dune::MatrixAdapter<BCRSMatrix,BlockVector,BlockVector>
MatrixOperator;
typedef ScalingLinearOperator<BlockVector> ScalingOperator;
typedef LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
typedef Impl::ScalingLinearOperator<BlockVector> ScalingOperator;
typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
public:
//! Type of underlying field
......@@ -1078,7 +1032,8 @@ namespace Dune
const std::string blank_;
};
} // namespace Dune
/** @} */
} // namespace Dune
#endif // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
dune/istl/paamg/amg.hh
View file @ 77390949
......@@ -378,7 +378,9 @@ namespace Dune
enum SolverType { umfpack, superlu, none };
static constexpr SolverType solver =
#if HAVE_SUITESPARSE_UMFPACK
#if DISABLE_AMG_DIRECTSOLVER
none;
#elif HAVE_SUITESPARSE_UMFPACK
UMFPackMethodChooser< field_type > :: valid ? umfpack : none ;
#elif HAVE_SUPERLU
superlu ;
......
dune/istl/solver.hh
View file @ 77390949
......@@ -6,6 +6,7 @@
#include <iomanip>
#include <ostream>
#include "solvertype.hh"
namespace Dune
{
......@@ -155,6 +156,54 @@ namespace Dune
}
};
/**
* \brief Helper class for notifying a DUNE-ISTL linear solver about
* a change of the iteration matrix object in a unified way,
* i.e. independent from the solver's type (direct/iterative).
*
* \author Sebastian Westerheide.
*/
template <typename ISTLLinearSolver, typename BCRSMatrix>
class SolverHelper
{
public:
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
{
static const bool is_direct_solver
= Dune::IsDirectSolver<ISTLLinearSolver>::value;
SolverHelper<ISTLLinearSolver,BCRSMatrix>::
Implementation<is_direct_solver>::setMatrix(solver,matrix);
}
protected:
/**
* \brief Implementation that works together with iterative ISTL
* solvers, e.g. Dune::CGSolver or Dune::BiCGSTABSolver.
*/
template <bool is_direct_solver, typename Dummy = void>
struct Implementation
{
static void setMatrix (ISTLLinearSolver&,
const BCRSMatrix&)
{}
};
/**
* \brief Implementation that works together with direct ISTL
* solvers, e.g. Dune::SuperLU or Dune::UMFPack.
*/
template <typename Dummy>
struct Implementation<true,Dummy>
{
static void setMatrix (ISTLLinearSolver& solver,
const BCRSMatrix& matrix)
{
solver.setMatrix(matrix);
}
};
};
/**
* @}
*/
......
dune/istl/solvers.hh
View file @ 77390949
......@@ -12,6 +12,7 @@
#include <string>
#include <vector>
#include <array>
#include <type_traits>
#include "istlexception.hh"
#include "operators.hh"
......@@ -20,19 +21,16 @@
#include "preconditioner.hh"
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/conditional.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/timer.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/typetraits.hh>
namespace Dune {
/** @defgroup ISTL_Solvers Iterative Solvers
@ingroup ISTL
*/
/** @addtogroup ISTL_Solvers
@{
*/
/** \file
\brief Implementations of the inverse operator interface.
......@@ -40,9 +38,7 @@ namespace Dune {
This file provides various preconditioned Krylov methods.
*/
//=====================================================================
//=====================================================================
// Implementation of this interface
//=====================================================================
......@@ -173,7 +169,7 @@ namespace Dune {
this->printOutput(std::cout,i,defnew,def);
//std::cout << i << " " << defnew << " " << defnew/def << std::endl;
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
{
res.converged = true;
break;
......@@ -192,7 +188,7 @@ namespace Dune {
// fill statistics
res.iterations = i;
res.reduction = def/def0;
res.reduction = max_value(def/def0);
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
......@@ -315,7 +311,7 @@ namespace Dune {
this->printOutput(std::cout,i,defnew,def);
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
{
res.converged = true;
break;
......@@ -330,8 +326,8 @@ namespace Dune {
_prec.post(x); // postprocess preconditioner
res.iterations = i; // fill statistics
res.reduction = static_cast<double>(def/def0);
res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
res.reduction = static_cast<double>(max_value(def/def0));
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
if (_verbose>0) // final print
std::cout << "=== rate=" << res.conv_rate
......@@ -432,7 +428,7 @@ namespace Dune {
real_type def0 = _sp.norm(b); // compute norm
if (!isfinite(def0)) // check for inf or NaN
if (!all_true(isfinite(def0))) // check for inf or NaN
{
if (_verbose>0)
std::cout << "=== CGSolver: abort due to infinite or NaN initial defect"
......@@ -441,7 +437,7 @@ namespace Dune {
<< " is infinite or NaN");
}
if (def0<1E-30) // convergence check
if (max_value(def0)<1E-30) // convergence check
{
res.converged = true;
res.iterations = 0; // fill statistics
......@@ -460,7 +456,7 @@ namespace Dune {
std::cout << "=== CGSolver" << std::endl;
if (_verbose>1) {
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),def0);
this->printOutput(std::cout,0,def0);
}
}
......@@ -488,10 +484,10 @@ namespace Dune {
real_type defnew=_sp.norm(b); // comp defect norm
if (_verbose>1) // print
this->printOutput(std::cout,real_type(i),defnew,def);
this->printOutput(std::cout,i,defnew,def);
def = defnew; // update norm
if (!isfinite(def)) // check for inf or NaN
if (!all_true(isfinite(def))) // check for inf or NaN
{
if (_verbose>0)
std::cout << "=== CGSolver: abort due to infinite or NaN defect"
......@@ -500,7 +496,7 @@ namespace Dune {
"CGSolver: defect=" << def << " is infinite or NaN");
}
if (def<def0*_reduction || def<1E-30) // convergence check
if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
{
res.converged = true;
break;
......@@ -520,12 +516,12 @@ namespace Dune {
i=std::min(_maxit,i);
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,real_type(i),def);
this->printOutput(std::cout,i,def);
_prec.post(x); // postprocess preconditioner
res.iterations = i; // fill statistics
res.reduction = static_cast<double>(def/def0);
res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
res.iterations = i; // fill statistics
res.reduction = static_cast<double>(max_value(max_value(def/def0)));
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
if (_verbose>0) // final print
......@@ -624,6 +620,7 @@ namespace Dune {
{
using std::abs;
const real_type EPSILON=1e-80;
using std::abs;
double it;
field_type rho, rho_new, alpha, beta, h, omega;
real_type norm, norm_old, norm_0;
......@@ -671,7 +668,7 @@ namespace Dune {
}
}
if ( norm < (_reduction * norm_0) || norm<1E-30)
if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
{
res.converged = 1;
_prec.post(x); // postprocess preconditioner
......@@ -696,13 +693,13 @@ namespace Dune {
rho_new = _sp.dot(rt,r);
// look if breakdown occurred
if (abs(rho) <= EPSILON)
if (all_true(abs(rho) <= EPSILON))
DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - rho "
<< rho << " <= EPSILON " << EPSILON
<< rho << " <= EPSILON " << max_value(EPSILON)
<< " after " << it << " iterations");
if (abs(omega) <= EPSILON)
if (all_true(abs(omega) <= EPSILON))
DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - omega "
<< omega << " <= EPSILON " << EPSILON
<< omega << " <= EPSILON " << max_value(EPSILON)
<< " after " << it << " iterations");
......@@ -726,9 +723,9 @@ namespace Dune {
// alpha = rho_new / < rt, v >
h = _sp.dot(rt,v);
if (abs(h) < EPSILON)
if ( all_true(abs(h) < EPSILON) )
DUNE_THROW(SolverAbort,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
<< abs(h) << " < EPSILON " << EPSILON
<< abs(h) << " < EPSILON " << max_value(EPSILON)
<< " after " << it << " iterations");
alpha = rho_new / h;
......@@ -751,7 +748,7 @@ namespace Dune {
this->printOutput(std::cout,it,norm,norm_old);
}
if ( norm < (_reduction * norm_0) )
if ( all_true(norm < (_reduction * norm_0)) )
{
res.converged = 1;
break;
......@@ -790,7 +787,7 @@ namespace Dune {
this->printOutput(std::cout,it,norm,norm_old);
}
if ( norm < (_reduction * norm_0) || norm<1E-30)
if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
{
res.converged = 1;
break;
......@@ -807,8 +804,8 @@ namespace Dune {
_prec.post(x); // postprocess preconditioner
res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
res.reduction = static_cast<double>(norm/norm_0);
res.conv_rate = static_cast<double>(pow(res.reduction,1.0/it));
res.reduction = static_cast<double>(max_value(norm/norm_0));
res.conv_rate = pow(res.reduction,1.0/it);
res.elapsed = watch.elapsed();
if (_verbose>0) // final print
std::cout << "=== rate=" << res.conv_rate
......@@ -921,7 +918,7 @@ namespace Dune {
}
// check for convergence
if(def0 < 1e-30 ) {
if( max_value(def0) < 1e-30 ) {
res.converged = true;
res.iterations = 0;
res.reduction = 0;
......@@ -1048,7 +1045,8 @@ namespace Dune {
this->printOutput(std::cout,i,defnew,def);
def = defnew;
if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
if(all_true(def < def0*_reduction)
|| max_value(def) < 1e-30 || i == _maxit ) {
res.converged = true;
break;
}
......@@ -1061,8 +1059,8 @@ namespace Dune {
_prec.post(x);
// fill statistics
res.iterations = i;
res.reduction = static_cast<double>(def/def0);
res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
res.reduction = static_cast<double>(max_value(def/def0));
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
// final print
......@@ -1089,35 +1087,50 @@ namespace Dune {
private:
// helper function to extract the real value of a real or complex number
inline
real_type to_real(const real_type & v)
{
return v;
}
inline
real_type to_real(const std::complex<real_type> & v)
{
return v.real();
}
void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
{
using std::sqrt;
using std::abs;
using std::max;
using std::min;
const real_type eps = 1e-15;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
if(norm_dy < 1e-15) {
cs = 1.0;
sn = 0.0;
} else if(norm_dx < 1e-15) {
cs = 0.0;
sn = 1.0;
} else if(norm_dy > norm_dx) {
real_type temp = norm_dx/norm_dy;
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
cs *= temp;
sn *= dx/norm_dx;
// dy is real in exact arithmetic
// so we don't need to conjugate here
sn *= dy/norm_dy;
} else {
real_type temp = norm_dy/norm_dx;
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
sn *= dy/dx;
// dy and dx is real in exact arithmetic
// so we don't have to conjugate both of them
}
real_type norm_max = max(norm_dx, norm_dy);
real_type norm_min = min(norm_dx, norm_dy);
real_type temp = norm_min/norm_max;
// we rewrite the code in a vectorizable fashion
cs = cond(norm_dy < eps,
real_type(1.0),
cond(norm_dx < eps,
real_type(0.0),
cond(norm_dy > norm_dx,
real_type(1.0)/sqrt(1.0 + temp*temp)*temp,
// dy and dx are real in exact arithmetic
// thus dx*dy is real so we can explicitly enforce it
real_type(1.0)/sqrt(1.0 + temp*temp)*to_real(dx*dy)/norm_dx/norm_dy)));
sn = cond(norm_dy < eps,
field_type(0.0),
cond(norm_dx < eps,
field_type(1.0),
cond(norm_dy > norm_dx,
field_type(1.0)/sqrt(1.0 + temp*temp),
// dy and dx is real in exact arithmetic
// so we don't have to conjugate both of them
field_type(1.0)/sqrt(1.0 + temp*temp)*dy/dx)));
}
SeqScalarProduct<X> ssp;
......@@ -1239,7 +1252,7 @@ namespace Dune {
*/
virtual void apply (X& x, Y& b, InverseOperatorResult& res)
{
apply(x,b,_reduction,res);
apply(x,b,max_value(_reduction),res);
}
/*!
......@@ -1292,7 +1305,7 @@ namespace Dune {
}
}
if(norm_0 < EPSILON) {
if(all_true(norm_0 < EPSILON)) {
_W.post(x);
res.converged = true;
if(_verbose > 0) // final print
......@@ -1324,7 +1337,7 @@ namespace Dune {
w.axpy(-H[k][i],v[k]);
}
H[i+1][i] = _sp.norm(w);
if(abs(H[i+1][i]) < EPSILON)
if(all_true(abs(H[i+1][i]) < EPSILON))
DUNE_THROW(SolverAbort,
"breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
......@@ -1352,7 +1365,7 @@ namespace Dune {
norm_old = norm;
// check convergence
if(norm < reduction * norm_0)
if(all_true(norm < reduction * norm_0))
res.converged = true;
} // end for
......@@ -1388,8 +1401,8 @@ namespace Dune {
// save solver statistics
res.iterations = j-1; // it has to be j-1!!!
res.reduction = static_cast<double>(norm/norm_0);
res.conv_rate = static_cast<double>(pow(res.reduction,1.0/(j-1)));
res.reduction = static_cast<double>(max_value(norm/norm_0));
res.conv_rate = pow(res.reduction,1.0/(j-1));
res.elapsed = watch.elapsed();
if(_verbose>0)
......@@ -1435,35 +1448,51 @@ namespace Dune {
template<typename T>
typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
using std::conj;
return conj(t);
}
// helper function to extract the real value of a real or complex number
inline
real_type to_real(const real_type & v)
{
return v;
}
inline
real_type to_real(const std::complex<real_type> & v)
{
return v.real();
}
void
generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
{
using std::sqrt;
using std::abs;
using std::max;
using std::min;
const real_type eps = 1e-15;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
if(norm_dy < 1e-15) {
cs = 1.0;
sn = 0.0;
} else if(norm_dx < 1e-15) {
cs = 0.0;
sn = 1.0;
} else if(norm_dy > norm_dx) {
real_type temp = norm_dx/norm_dy;
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
cs *= temp;
sn *= dx/norm_dx;
sn *= conjugate(dy)/norm_dy;
} else {
real_type temp = norm_dy/norm_dx;
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
sn *= conjugate(dy/dx);
}
real_type norm_max = max(norm_dx, norm_dy);
real_type norm_min = min(norm_dx, norm_dy);
real_type temp = norm_min/norm_max;
// we rewrite the code in a vectorizable fashion
cs = cond(norm_dy < eps,
real_type(1.0),
cond(norm_dx < eps,
real_type(0.0),
cond(norm_dy > norm_dx,
real_type(1.0)/sqrt(1.0 + temp*temp)*temp,
real_type(1.0)/sqrt(1.0 + temp*temp)*to_real(dx*conjugate(dy))/norm_dx/norm_dy)));
sn = cond(norm_dy < eps,
field_type(0.0),
cond(norm_dx < eps,
field_type(1.0),
cond(norm_dy > norm_dx,
field_type(1.0)/sqrt(1.0 + temp*temp),
field_type(1.0)/sqrt(1.0 + temp*temp)*conjugate(dy/dx))));
}
......@@ -1569,7 +1598,7 @@ namespace Dune {
p[0].reset(new X(x));
real_type def0 = _sp.norm(b); // compute norm
if (def0<1E-30) // convergence check
if ( max_value(def0) < 1E-30 ) // convergence check
{
res.converged = true;
res.iterations = 0; // fill statistics
......@@ -1612,7 +1641,7 @@ namespace Dune {
if (_verbose>1) // print
this->printOutput(std::cout,++i,defnew,def);
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
{
res.converged = true;
if (_verbose>0) // final print
......@@ -1658,7 +1687,7 @@ namespace Dune {
this->printOutput(std::cout,++i,defNew,def);
def = defNew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
{
res.converged = true;
break;
......@@ -1677,7 +1706,7 @@ namespace Dune {
// fill statistics
res.iterations = i;
res.reduction = def/def0;
res.reduction = max_value(def/def0);
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
......
dune/istl/superlu.hh
View file @ 77390949
......@@ -4,32 +4,8 @@
#define DUNE_ISTL_SUPERLU_HH
#if HAVE_SUPERLU
#ifdef TRUE
#undef TRUE
#endif
#ifdef FALSE
#undef FALSE
#endif
#ifndef SUPERLU_NTYPE
#define SUPERLU_NTYPE 1
#endif
#if SUPERLU_NTYPE==0
#include "slu_sdefs.h"
#endif
#if SUPERLU_NTYPE==1
#include "slu_ddefs.h"
#endif
#if SUPERLU_NTYPE==2
#include "slu_cdefs.h"
#endif
#if SUPERLU_NTYPE>=3
#include "slu_zdefs.h"
#endif
#include "superlufunctions.hh"
#include "solvers.hh"
#include "supermatrix.hh"
#include <algorithm>
......@@ -81,7 +57,7 @@ namespace Dune
struct QuerySpaceChooser
{};
#if SUPERLU_NTYPE==0
#if HAVE_SLU_SDEFS_H
template<>
struct SuperLUDenseMatChooser<float>
{
......@@ -129,7 +105,7 @@ namespace Dune
#endif
#if SUPERLU_NTYPE==1
#if HAVE_SLU_DDEFS_H
template<>
struct SuperLUDenseMatChooser<double>
......@@ -176,7 +152,7 @@ namespace Dune
};
#endif
#if SUPERLU_NTYPE>=3
#if HAVE_SLU_ZDEFS_H
template<>
struct SuperLUDenseMatChooser<std::complex<double> >
{
......@@ -223,7 +199,7 @@ namespace Dune
};
#endif
#if SUPERLU_NTYPE==2
#if HAVE_SLU_CDEFS_H
template<>
struct SuperLUDenseMatChooser<std::complex<float> >
{
......
dune/istl/superlufunctions.hh 0 → 100644
View file @ 77390949
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_SUPERLUFUNCTIONS_HH
#define DUNE_ISTL_SUPERLUFUNCTIONS_HH
#if HAVE_SUPERLU
#define int_t SUPERLU_INT_TYPE
#include "supermatrix.h"
#include "slu_util.h"
#undef int_t
#if HAVE_SLU_SDEFS_H
extern "C" {
extern void
sgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
char *, float *, float *, SuperMatrix *, SuperMatrix *,
void *, int, SuperMatrix *, SuperMatrix *,
float *, float *, float *, float *,
#if SUPERLU_MIN_VERSION_5
GlobalLU_t*,
#endif
mem_usage_t *, SuperLUStat_t *, int *);
extern void
sCreate_Dense_Matrix(SuperMatrix *, int, int, float *, int,
Stype_t, Dtype_t, Mtype_t);
extern void
sCreate_CompCol_Matrix(SuperMatrix *, int, int, int, float *,
int *, int *, Stype_t, Dtype_t, Mtype_t);
extern int sQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
extern void sPrint_CompCol_Matrix(char *, SuperMatrix *);
}
#endif
#if HAVE_SLU_DDEFS_H
extern "C" {
extern void
dgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
char *, double *, double *, SuperMatrix *, SuperMatrix *,
void *, int, SuperMatrix *, SuperMatrix *,
double *, double *, double *, double *,
#if SUPERLU_MIN_VERSION_5
GlobalLU_t*,
#endif
mem_usage_t *, SuperLUStat_t *, int *);
extern void
dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,
int *, int *, Stype_t, Dtype_t, Mtype_t);
extern void
dCreate_Dense_Matrix(SuperMatrix *, int, int, double *, int,
Stype_t, Dtype_t, Mtype_t);
extern int dQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
extern void dPrint_CompCol_Matrix(char *, SuperMatrix *);
}
#endif
#if HAVE_SLU_CDEFS_H
#include "slu_scomplex.h"
extern "C" {
extern void
cgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
char *, float *, float *, SuperMatrix *, SuperMatrix *,
void *, int, SuperMatrix *, SuperMatrix *,
float *, float *, float *, float *,
#if SUPERLU_MIN_VERSION_5
GlobalLU_t*,
#endif
mem_usage_t *, SuperLUStat_t *, int *);
extern void
cCreate_Dense_Matrix(SuperMatrix *, int, int, ::complex *, int,
Stype_t, Dtype_t, Mtype_t);
extern void
cCreate_CompCol_Matrix(SuperMatrix *, int, int, int, ::complex *,
int *, int *, Stype_t, Dtype_t, Mtype_t);
extern int cQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
extern void cPrint_CompCol_Matrix(char *, SuperMatrix *);
}
#endif
#if HAVE_SLU_ZDEFS_H
#include "slu_dcomplex.h"
extern "C" {
extern void
zgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
char *, double *, double *, SuperMatrix *, SuperMatrix *,
void *, int, SuperMatrix *, SuperMatrix *,
double *, double *, double *, double *,
#if SUPERLU_MIN_VERSION_5
GlobalLU_t*,
#endif
mem_usage_t *, SuperLUStat_t *, int *);
extern void
zCreate_CompCol_Matrix(SuperMatrix *, int, int, int, doublecomplex *,
int *, int *, Stype_t, Dtype_t, Mtype_t);
extern void
zCreate_Dense_Matrix(SuperMatrix *, int, int, doublecomplex *, int,
Stype_t, Dtype_t, Mtype_t);
extern int zQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
extern void zPrint_CompCol_Matrix(char *, SuperMatrix *);
}
#endif
#endif
#endif
dune/istl/supermatrix.hh
View file @ 77390949
......@@ -5,26 +5,6 @@
#if HAVE_SUPERLU
#ifndef SUPERLU_NTYPE
#define SUPERLU_NTYPE 1
#endif
#if SUPERLU_NTYPE==0
#include "slu_sdefs.h"
#endif
#if SUPERLU_NTYPE==1
#include "slu_ddefs.h"
#endif
#if SUPERLU_NTYPE==2
#include "slu_cdefs.h"
#endif
#if SUPERLU_NTYPE>=3
#include "slu_zdefs.h"
#endif
#include "bcrsmatrix.hh"
#include "bvector.hh"
#include <dune/common/fmatrix.hh>
......@@ -34,6 +14,8 @@
#include"colcompmatrix.hh"
#include "superlufunctions.hh"
namespace Dune
{
......@@ -45,7 +27,7 @@ namespace Dune
struct SuperMatrixPrinter
{};
#if SUPERLU_NTYPE==0
#if HAVE_SLU_SDEFS_H
template<>
struct SuperMatrixCreateSparseChooser<float>
{
......@@ -68,7 +50,7 @@ namespace Dune
};
#endif
#if SUPERLU_NTYPE==1
#if HAVE_SLU_DDEFS_H
template<>
struct SuperMatrixCreateSparseChooser<double>
{
......@@ -91,7 +73,7 @@ namespace Dune
};
#endif
#if SUPERLU_NTYPE==2
#if HAVE_SLU_CDEFS_H
template<>
struct SuperMatrixCreateSparseChooser<std::complex<float> >
{
......@@ -114,7 +96,7 @@ namespace Dune
};
#endif
#if SUPERLU_NTYPE>=3
#if HAVE_SLU_ZDEFS_H
template<>
struct SuperMatrixCreateSparseChooser<std::complex<double> >
{
......
dune/istl/test/CMakeLists.txt
View file @ 77390949
......@@ -20,6 +20,11 @@ dune_add_test(SOURCES matrixutilstest.cc)
dune_add_test(SOURCES matrixtest.cc)
if(Vc_FOUND)
dune_add_test(SOURCES multirhstest.cc)
add_dune_vc_flags(multirhstest)
endif(Vc_FOUND)
dune_add_test(SOURCES bvectortest.cc)
dune_add_test(SOURCES vbvectortest.cc)
......@@ -50,35 +55,35 @@ dune_add_test(SOURCES solvertest.cc)
dune_add_test(SOURCES solveraborttest.cc)
# Pardiso tests
dune_add_test(SOURCES test_pardiso.cc)
dune_add_test(SOURCES test_pardiso.cc)
# SuperLU tests
dune_add_test(NAME superlustest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=0)
dune_add_test(NAME superlustest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=0)
dune_add_test(SOURCES superlutest.cc)
dune_add_test(SOURCES superlutest.cc)
dune_add_test(NAME superluctest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=2)
dune_add_test(NAME superluctest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=2)
dune_add_test(NAME superluztest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
dune_add_test(NAME superluztest
SOURCES superlutest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
dune_add_test(SOURCES complexrhstest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
dune_add_test(SOURCES complexrhstest.cc
COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
# SuiteSparse tests
dune_add_test(SOURCES umfpacktest.cc)
set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "umfpack_decomp")
dune_add_test(SOURCES umfpacktest.cc)
set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "umfpack_decomp")
dune_add_test(SOURCES ldltest.cc)
dune_add_test(SOURCES spqrtest.cc)
dune_add_test(SOURCES overlappingschwarztest.cc)
dune_add_test(SOURCES overlappingschwarztest.cc)
# MPI tests
dune_add_test(SOURCES matrixredisttest.cc
......@@ -87,6 +92,6 @@ dune_add_test(SOURCES matrixredisttest.cc
dune_add_test(SOURCES vectorcommtest.cc
CMAKE_GUARD MPI_FOUND)
dune_add_test(SOURCES matrixmarkettest.cc)
dune_add_test(SOURCES matrixmarkettest.cc)
exclude_from_headercheck(complexdata.hh)
dune/istl/test/multirhstest.cc 0 → 100644
View file @ 77390949
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// start with including some headers
#include "config.h"
#define DISABLE_AMG_DIRECTSOLVER 1
// #undef HAVE_VC
#include <iostream> // for input/output to shell
#include <fstream> // for input/output to files
#include <vector> // STL vector class
#include <complex>
#include <cmath> // Yes, we do some math here
#include <sys/times.h> // for timing measurements
#include <dune/common/classname.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/simd.hh>
#include <dune/common/timer.hh>
#include <dune/istl/istlexception.hh>
#include <dune/istl/basearray.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include "laplacian.hh"
template<typename T>
struct Random {
static T gen()
{
return drand48();
}
};
#if HAVE_VC
template<typename T, typename A>
struct Random<Vc::Vector<T,A>> {
static Vc::Vector<T,A> gen()
{
return Vc::Vector<T,A>::Random();
}
};
template<typename T, std::size_t N, typename V, std::size_t M>
struct Random<Vc::SimdArray<T,N,V,M>> {
static Vc::SimdArray<T,N,V,M> gen()
{
return Vc::SimdArray<T,N,V,M>::Random();
}
};
#endif
template <typename V>
V detectVectorType(Dune::LinearOperator<V,V> &);
template<typename Operator, typename Solver>
void run_test (std::string precName, std::string solverName, Operator & op, Solver & solver, unsigned int N, unsigned int Runs)
{
using Vector = decltype(detectVectorType(op));
using FT = typename Vector::field_type;
Dune::Timer t;
std::cout << "Trying " << solverName << "(" << precName << ")"
<< " with " << Dune::className<FT>() << std::endl;
for (unsigned int run = 0; run < Runs; run++) {
// set up system
Vector x(N),b(N);
for (unsigned int i=0; i<N; i++)
x[i] += Random<FT>::gen();
b=0; op.apply(x,b); // set right hand side accordingly
x=1; // initial guess
// call the solver
Dune::InverseOperatorResult r;
solver.apply(x,b,r);
}
std::cout << Runs << " run(s) took " << t.stop() << std::endl;
}
template<typename Operator, typename Prec>
void test_all_solvers(std::string precName, Operator & op, Prec & prec, unsigned int N, unsigned int Runs)
{
using Vector = decltype(detectVectorType(op));
double reduction = 1e-4;
int verb = 1;
Dune::LoopSolver<Vector> loop(op,prec,reduction,18000,verb);
Dune::CGSolver<Vector> cg(op,prec,reduction,8000,verb);
Dune::BiCGSTABSolver<Vector> bcgs(op,prec,reduction,8000,verb);
Dune::GradientSolver<Vector> grad(op,prec,reduction,18000,verb);
Dune::RestartedGMResSolver<Vector> gmres(op,prec,reduction,40,8000,verb);
Dune::MINRESSolver<Vector> minres(op,prec,reduction,8000,verb);
Dune::GeneralizedPCGSolver<Vector> gpcg(op,prec,reduction,8000,verb);
// run_test(precName, "Loop", op,loop,N,Runs);
run_test(precName, "CG", op,cg,N,Runs);
run_test(precName, "BiCGStab", op,bcgs,N,Runs);
run_test(precName, "Gradient", op,grad,N,Runs);
run_test(precName, "RestartedGMRes", op,gmres,N,Runs);
run_test(precName, "MINRes", op,minres,N,Runs);
run_test(precName, "GeneralizedPCG", op,gpcg,N,Runs);
}
template<typename FT>
void test_all(unsigned int Runs = 1)
{
// define Types
typedef Dune::FieldVector<FT,1> VB;
typedef Dune::FieldMatrix<double,1,1> MB;
typedef Dune::BlockVector<VB> Vector;
typedef Dune::BCRSMatrix<MB> Matrix;
// size
unsigned int size = 100;
unsigned int N = size*size;
// make a compressed row matrix with five point stencil
Matrix A;
setupLaplacian(A,size);
typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
Operator op(A); // make linear operator from A
// create all preconditioners
Dune::SeqJac<Matrix,Vector,Vector> jac(A,1,0.1); // Jacobi preconditioner
Dune::SeqGS<Matrix,Vector,Vector> gs(A,1,0.1); // GS preconditioner
Dune::SeqSOR<Matrix,Vector,Vector> sor(A,1,0.1); // SOR preconditioner
Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,0.1); // SSOR preconditioner
Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,0.1); // preconditioner object
Dune::SeqILUn<Matrix,Vector,Vector> ilu1(A,1,0.1); // preconditioner object
// AMG
typedef Dune::Amg::RowSum Norm;
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Matrix,Norm> >
Criterion;
typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = 1;
unsigned int coarsenTarget = 1000;
unsigned int maxLevel = 10;
Criterion criterion(15,coarsenTarget);
criterion.setDefaultValuesIsotropic(2);
criterion.setAlpha(.67);
criterion.setBeta(1.0e-4);
criterion.setMaxLevel(maxLevel);
criterion.setSkipIsolated(false);
criterion.setNoPreSmoothSteps(1);
criterion.setNoPostSmoothSteps(1);
Dune::SeqScalarProduct<Vector> sp;
typedef Dune::Amg::AMG<Operator,Vector,Smoother> AMG;
Smoother smoother(A,1,1);
AMG amg(op, criterion, smootherArgs);
// run the sub-tests
test_all_solvers("Jacobi", op,jac,N,Runs);
test_all_solvers("GaussSeidel", op,gs,N,Runs);
test_all_solvers("SOR", op,sor,N,Runs);
test_all_solvers("SSOR", op,ssor,N,Runs);
test_all_solvers("ILU0", op,ilu0,N,Runs);
test_all_solvers("ILU1", op,ilu1,N,Runs);
test_all_solvers("AMG", op,amg,N,Runs);
}
int main (int argc, char ** argv)
{
test_all<float>();
test_all<double>();
#if HAVE_VC
// test_all<Vc::float_v>();
test_all<Vc::double_v>();
test_all<Vc::Vector<double, Vc::VectorAbi::Scalar>>();
test_all<Vc::SimdArray<double,2>>();
test_all<Vc::SimdArray<double,2,Vc::Vector<double, Vc::VectorAbi::Scalar>,1>>();
test_all<Vc::SimdArray<double,8>>();
test_all<Vc::SimdArray<double,8,Vc::Vector<double, Vc::VectorAbi::Scalar>,1>>();
#endif
test_all<double>(8);
return 0;
}
dune/istl/test/superlutest.cc
View file @ 77390949
......@@ -14,36 +14,11 @@
#include "laplacian.hh"
#ifndef SUPERLU_NTYPE
#define SUPERLU_NTYPE 1
#endif
#if SUPERLU_NTYPE==1
typedef double FIELD_TYPE;
#endif
#if SUPERLU_NTYPE==0
typedef float FIELD_TYPE;
#endif
#if SUPERLU_NTYPE==2
typedef std::complex<float> FIELD_TYPE;
#endif
#if SUPERLU_NTYPE>=3
typedef std::complex<double> FIELD_TYPE;
#endif
int main(int argc, char** argv)
try
{
#if HAVE_SUPERLU
template<typename FIELD_TYPE>
void runSuperLU(std::size_t N)
{
const int BS=1;
std::size_t N=100;
if(argc>1)
N = atoi(argv[1]);
std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
typedef Dune::FieldMatrix<FIELD_TYPE,BS,BS> MatrixBlock;
......@@ -83,8 +58,34 @@ try
solver1.apply(x1,b1, res);
solver1.apply(reinterpret_cast<FIELD_TYPE*>(&x1[0]), reinterpret_cast<FIELD_TYPE*>(&b1[0]));
}
#endif
int main(int argc, char** argv)
try
{
std::size_t N=100;
if(argc>1)
N = atoi(argv[1]);
std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
#if HAVE_SUPERLU
#if HAVE_SLU_SDEFS_H
runSuperLU<float>(N);
#endif
#if HAVE_SLU_DDEFS_H
runSuperLU<double>(N);
#endif
#if HAVE_SLU_CDEFS_H
runSuperLU<std::complex<float> >(N);
#endif
#if HAVE_SLU_ZDEFS_H
runSuperLU<std::complex<double> >(N);
#endif
return 0;
#else // HAVE_SUPERLU
std::cerr << "You need SuperLU to run this test." << std::endl;
return 77;
......