Refactor UMFPack test

To prepare for testing with more matrix types.

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;
+}