Merge branch 'issue/spqr_block_matrix' into 'master'

spqr solver with block matrix

See merge request core/dune-istl!409

dune/istl/spqr.hh

@@ -152,17 +152,23 @@ namespace Dune {
     /** \copydoc InverseOperator::apply(X&, Y&, InverseOperatorResult&) */
     virtual void apply(domain_type& x, range_type& b, InverseOperatorResult& res)
     {
-      const std::size_t dimMat(spqrMatrix_.N());
+      const std::size_t numRows(spqrMatrix_.N());
       // fill B
-      for(std::size_t k = 0; k != dimMat; ++k)
-        (static_cast<T*>(B_->x))[k] = b[k];
+      for(std::size_t k = 0; k != numRows/n; ++k)
+        for (int l = 0; l < n; ++l)
+          (static_cast<T*>(B_->x))[n*k+l] = b[k][l];
+
       cholmod_dense* BTemp = B_;
       B_ = SuiteSparseQR_qmult<T>(0, spqrfactorization_, B_, cc_);
       cholmod_dense* X = SuiteSparseQR_solve<T>(1, spqrfactorization_, B_, cc_);
       cholmod_l_free_dense(&BTemp, cc_);
+
+      const std::size_t numCols(spqrMatrix_.M());
       // fill x
-      for(std::size_t k = 0; k != dimMat; ++k)
-        x [k] = (static_cast<T*>(X->x))[k];
+      for(std::size_t k = 0; k != numCols/m; ++k)
+        for (int l = 0; l < m; ++l)
+          x[k][l] = (static_cast<T*>(X->x))[m*k+l];
+
       cholmod_l_free_dense(&X, cc_);
       // this is a direct solver
       res.iterations = 1;

dune/istl/test/spqrtest.cc

@@ -12,64 +12,69 @@
 
 #include "laplacian.hh"
 
-int main(int argc, char** argv)
+template <class FIELD_TYPE, int BS>
+void run(std::size_t N)
 {
 #if HAVE_SUITESPARSE_SPQR
-  try
-  {
-    typedef double FIELD_TYPE;
-    //typedef std::complex<double> FIELD_TYPE;
+  std::cout << "testing for N=" << N << " BS=" << BS << std::endl;
 
-    const int BS=1;
-    std::size_t N=100;
+  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;
 
-    if (argc > 1)
-    {
-      N = atoi(argv[1]);
-    }
-    std::cout << "testing for N=" << N << " BS=" << BS << std::endl;
+  BCRSMat mat;
+  Operator fop(mat);
+  Vector b(N*N), x(N*N);
 
-    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;
+  setupLaplacian(mat,N);
+  b=1;
+  x=0;
 
-    BCRSMat mat;
-    Operator fop(mat);
-    Vector b(N*N), x(N*N);
+  Dune::Timer watch;
 
-    setupLaplacian(mat,N);
-    b=1;
-    x=0;
+  watch.reset();
 
-    Dune::Timer watch;
+  Dune::SPQR<BCRSMat> solver(mat,1);
 
-    watch.reset();
+  Dune::InverseOperatorResult res;
 
-    Dune::SPQR<BCRSMat> solver(mat,1);
+  Dune::SPQR<BCRSMat> solver1;
 
-    Dune::InverseOperatorResult res;
+  std::set<std::size_t> mrs;
+  for(std::size_t s=0; s < N/2; ++s)
+    mrs.insert(s);
 
-    Dune::SPQR<BCRSMat> solver1;
+  solver1.setSubMatrix(mat,mrs);
+  solver1.setVerbosity(true);
 
-    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);
+  Vector residuum(N*N);
+  residuum=0;
+  fop.apply(x,residuum);
+  residuum-=b;
+  std::cout<<"Residuum : "<<residuum.two_norm()<<std::endl;
 
-    solver.apply(x,b,res);
-    solver.free();
+  solver1.apply(x,b,res);
+#endif
+}
 
-    Vector residuum(N*N);
-    residuum=0;
-    fop.apply(x,residuum);
-    residuum-=b;
-    std::cout<<"Residuum : "<<residuum.two_norm()<<std::endl;
+int main(int argc, char** argv)
+{
+#if HAVE_SUITESPARSE_SPQR
+  try
+  {
+    std::size_t N=100;
+    if (argc > 1)
+      N = atoi(argv[1]);
 
-    solver1.apply(x,b,res);
+    run<double,1>(N);
+    run<double,2>(N);
+    // run<std::complex<double>,1>(N); // does not work
+    // run<std::complex<double>,2>(N);
 
     return 0;
   }