From 8ff42cf9a9432165fbc5b22ac2f7550ff6922a0e Mon Sep 17 00:00:00 2001
From: Marco Agnese <ma2413@imperial.ac.uk>
Date: Mon, 17 Nov 2014 16:35:28 +0000
Subject: [PATCH] [suitesparse] Add test for LDL wrapper
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Signed-off-by: Christoph Grüninger <gruenich@dune-project.org>
---
dune/istl/test/ldltest.cc | 79 +++++++++++++++++++++++++++++++++++++++
1 file changed, 79 insertions(+)
create mode 100644 dune/istl/test/ldltest.cc
diff --git a/dune/istl/test/ldltest.cc b/dune/istl/test/ldltest.cc
new file mode 100644
index 000000000..7187e5aea
--- /dev/null
+++ b/dune/istl/test/ldltest.cc
@@ -0,0 +1,79 @@
+#include "config.h"
+#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/ldl.hh>
+
+#include "laplacian.hh"
+
+int main(int argc, char** argv)
+{
+ try
+ {
+ typedef double FIELD_TYPE;
+
+ 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;
+ typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+ typedef Dune::FieldVector<FIELD_TYPE,BS> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> Vector;
+ typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
+
+ BCRSMat mat;
+ Operator fop(mat);
+ Vector b(N*N), x(N*N);
+
+ setupLaplacian(mat,N);
+ b=1;
+ x=0;
+
+ Dune::Timer watch;
+
+ watch.reset();
+
+ Dune::LDL<BCRSMat> solver(mat,1);
+
+ Dune::InverseOperatorResult res;
+
+ Dune::LDL<BCRSMat> solver1;
+
+ std::set<std::size_t> mrs;
+ for(std::size_t s=0; s < N/2; ++s)
+ mrs.insert(s);
+
+ solver1.setSubMatrix(mat,mrs);
+ solver1.setVerbosity(true);
+
+ solver.apply(x,b, res);
+ solver.free();
+
+ Vector residuum(N*N);
+ residuum=0;
+ fop.apply(x,residuum);
+ residuum-=b;
+ std::cout<<"Residuum : "<<residuum.two_norm()<<std::endl;
+
+ solver1.apply(x,b, res);
+ solver1.apply(reinterpret_cast<FIELD_TYPE*>(&x[0]), reinterpret_cast<FIELD_TYPE*>(&b[0]));
+
+ return 0;
+ }
+ catch(Dune::Exception &e)
+ {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...)
+ {}
+}
--
GitLab