From 72f16adce25d95006d847943cfe658bb6c38ddad Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 22 Jan 2019 21:44:48 +0100
Subject: [PATCH] Refactor UMFPack test
To prepare for testing with more matrix types.
---
dune/istl/test/umfpacktest.cc | 101 ++++++++++++++++------------------
1 file changed, 47 insertions(+), 54 deletions(-)
diff --git a/dune/istl/test/umfpacktest.cc b/dune/istl/test/umfpacktest.cc
index 158dc4d05..60e824bc1 100644
--- a/dune/istl/test/umfpacktest.cc
+++ b/dune/istl/test/umfpacktest.cc
@@ -2,91 +2,84 @@
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
-#include <complex>
#include <iostream>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/timer.hh>
#include <dune/istl/bvector.hh>
-#include <dune/istl/colcompmatrix.hh>
-#include <dune/istl/io.hh>
-#include <dune/istl/operators.hh>
#include <dune/istl/umfpack.hh>
#include "laplacian.hh"
-int main(int argc, char** argv)
+template<typename Matrix, typename Vector>
+void runUMFPack(std::size_t N)
{
-#if HAVE_SUITESPARSE_UMFPACK
- try
- {
- typedef double FIELD_TYPE;
- //typedef std::complex<double> FIELD_TYPE;
+ typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
- const int BS=1;
- std::size_t N=100;
+ Matrix mat;
+ Operator fop(mat);
+ Vector b(N*N), x(N*N), b1(N/2), x1(N/2);
- if(argc>1)
- N = atoi(argv[1]);
- std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+ setupLaplacian(mat,N);
+ b=1;
+ b1=1;
+ x=0;
+ x1=0;
- typedef Dune::FieldMatrix<FIELD_TYPE,BS,BS> MatrixBlock;
- typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
- typedef Dune::FieldVector<FIELD_TYPE,BS> VectorBlock;
- typedef Dune::BlockVector<VectorBlock> Vector;
- typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
+ Dune::Timer watch;
- BCRSMat mat;
- Operator fop(mat);
- Vector b(N*N), x(N*N), b1(N/2), x1(N/2);
+ watch.reset();
- setupLaplacian(mat,N);
- b=1;
- b1=1;
- x=0;
- x1=0;
+ Dune::UMFPack<Matrix> solver(mat,1);
- Dune::Timer watch;
+ Dune::InverseOperatorResult res;
- watch.reset();
+ solver.apply(x, b, res);
+ solver.free();
- Dune::UMFPack<BCRSMat> solver(mat,1);
+ Dune::UMFPack<Matrix> solver1;
- Dune::InverseOperatorResult res;
+ std::set<std::size_t> mrs;
+ for(std::size_t s=0; s < N/2; ++s)
+ mrs.insert(s);
- solver.apply(x, b, res);
- solver.free();
+ solver1.setSubMatrix(mat,mrs);
+ solver1.setVerbosity(true);
- Dune::UMFPack<BCRSMat> solver1;
+ solver1.apply(x1,b1, res);
+ solver1.apply(reinterpret_cast<typename Matrix::field_type*>(&x1[0]),
+ reinterpret_cast<typename Matrix::field_type*>(&b1[0]));
- std::set<std::size_t> mrs;
- for(std::size_t s=0; s < N/2; ++s)
- mrs.insert(s);
+ Dune::UMFPack<Matrix> save_solver(mat,"umfpack_decomp",0);
+ Dune::UMFPack<Matrix> load_solver(mat,"umfpack_decomp",0);
+}
- solver1.setSubMatrix(mat,mrs);
- solver1.setVerbosity(true);
+int main(int argc, char** argv) try
+{
+#if HAVE_SUITESPARSE_UMFPACK
+ std::size_t N=100;
- solver1.apply(x1,b1, res);
- solver1.apply(reinterpret_cast<FIELD_TYPE*>(&x1[0]), reinterpret_cast<FIELD_TYPE*>(&b1[0]));
+ if(argc>1)
+ N = atoi(argv[1]);
- Dune::UMFPack<BCRSMat> save_solver(mat,"umfpack_decomp",0);
- Dune::UMFPack<BCRSMat> load_solver(mat,"umfpack_decomp",0);
- return 0;
- }
- catch (std::exception &e)
+ // ------------------------------------------------------------------------------
+ std::cout<<"testing for N="<<N<<" BCRSMatrix<FieldMatrix<double,1,1> >"<<std::endl;
{
- std::cout << "ERROR: " << e.what() << std::endl;
- return 1;
- }
- catch (...)
- {
- std::cerr << "Dune reported an unknown error." << std::endl;
- exit(1);
+ using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double,1> >;
+ runUMFPack<Matrix,Vector>(N);
}
+
+ return 0;
#else // HAVE_SUITESPARSE_UMFPACK
std::cerr << "You need SuiteSparse's UMFPack to run this test." << std::endl;
return 77;
#endif // HAVE_SUITESPARSE_UMFPACK
}
+catch (std::exception &e)
+{
+ std::cout << "ERROR: " << e.what() << std::endl;
+ return 1;
+}
--
GitLab