#222 ISTL solver causes segfault for special hierarchical matrix

Metadata

Property	Value
Reported by	Patrick Leidenberger (mailings.pl@gmail.com)
Reported at	Dec 13, 2006 13:27
Type	Bug Report
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Closed by	Markus Blatt (markus@dr-blatt.de)
Closed at	Dec 13, 2006 14:23
Closed in version	Unknown
Resolution	Fixed
Comment	Fixed in revision 707

Description

I create a three level nested matrix structure, the outer matrix is a 2x2 Matrix. If I fill only the [1][1] element with an BCRSMatrix the solver causes a segfault.

I modified the example.cc form istl/tutorial to demonstrate it. Here is the result:

[leidenberger@merlin00 tutorial]$ ./example ============================================. I [n=2,m=2,rowdim=6,coldim=6] row 0 1.00e+00 2.00e+00 . . . . row 1 2.00e+00 3.00e+00 . . . . row 2 . . 1.00e+00 2.00e+00 . . row 3 . . 2.00e+00 3.00e+00 . . row 4 . . . . 1.00e+00 2.00e+00 row 5 . . . . 2.00e+00 3.00e+00 =============================================. exact solution [blocks=2,dimension=6] entry 0 0.00e+00 0.00e+00 2.10e+01 0.00e+00 4.20e+01 0.00e+00 right hand side [blocks=2,dimension=6] entry 0 0.00e+00 0.00e+00 2.10e+01 4.20e+01 4.20e+01 8.40e+01 === LoopSolver Iter Defect Rate 0 1.0500E+02 Segmentation fault

And here is the code:

// start with including some headers #include "config.h"

#include // for input/output to shell #include // for input/output to files #include // STL vector class #include

#include<math.h> // Yes, we do some math here #include<stdio.h> // There is nothing better than sprintf #include<sys/times.h> // for timing measurements

#include<dune/istl/istlexception.hh> #include<dune/istl/basearray.hh> #include<dune/common/fvector.hh> #include<dune/common/fmatrix.hh> #include<dune/istl/bvector.hh> #include<dune/istl/vbvector.hh> #include<dune/istl/bcrsmatrix.hh> #include<dune/istl/io.hh> #include<dune/istl/gsetc.hh> #include<dune/istl/ilu.hh> #include<dune/istl/operators.hh> #include<dune/istl/solvers.hh> #include<dune/istl/preconditioners.hh> #include<dune/istl/scalarproducts.hh>

void test_Iter () {

const int GG=2, HH =3, II=2; typedef double ValueType; typedef Dune::BCRSMatrix<Dune::FieldMatrix<ValueType,1,1> > EmpMatrixTopoElemOrder; typedef Dune::BCRSMatrix EmpMatrixTopoElem; typedef Dune::BCRSMatrix EmpMatrixTopo;

typedef Dune::BlockVector<Dune::FieldVector<ValueType,1> > EmpVectorTopoElemOrder; typedef Dune::BlockVector EmpVectorTopoElem; typedef Dune::BlockVector EmpVectorTopo;

// Create matrix. Dune::FieldMatrix<ValueType,1,1> F; F = 1.0;

EmpMatrixTopoElemOrder G(GG,GG,GG*GG,EmpMatrixTopoElemOrder::row_wise); for (EmpMatrixTopoElemOrder::CreateIterator i=G.createbegin(); i!=G.createend(); ++i) for (int j=0; j<GG; ++j) i.insert(j);

for (EmpMatrixTopoElemOrder::RowIterator i=G.begin(); i!=G.end(); ++i) for (EmpMatrixTopoElemOrder::ColIterator j=(*i).begin(); j!=(*i).end(); ++j) *j = (F+i.index()+j.index());

// std::cout << "============================================." << std::endl; // printmatrix(std::cout,G,"G","row",10,2); // std::cout << "=============================================."<< std::endl;

EmpMatrixTopoElem H(HH,HH,HH,EmpMatrixTopoElem::row_wise); for (EmpMatrixTopoElem::CreateIterator i=H.createbegin(); i!=H.createend(); ++i) for (int j=0; j<HH; ++j) { if (i.index() == j) { i.insert(j); } }

for (EmpMatrixTopoElem::RowIterator i=H.begin(); i!=H.end(); ++i) for (EmpMatrixTopoElem::ColIterator j=(*i).begin(); j!=(*i).end(); ++j) *j = G;

// std::cout << "============================================." << std::endl; // printmatrix(std::cout,H,"H","row",10,2); // std::cout << "=============================================."<< std::endl;

EmpMatrixTopo I(II,II,1,EmpMatrixTopo::row_wise); for (EmpMatrixTopo::CreateIterator i=I.createbegin(); i!=I.createend(); ++i) for (int j=0; j<II; ++j) { if ((i.index() == 1) && (i.index() == j)) { i.insert(j); } }

for (EmpMatrixTopo::RowIterator i=I.begin(); i!=I.end(); ++i) for (EmpMatrixTopo::ColIterator j=(*i).begin(); j!=(*i).end(); ++j) *j = H;

std::cout << "============================================." << std::endl; printmatrix(std::cout,I,"I","row",10,2); std::cout << "=============================================."<< std::endl;

EmpVectorTopo x,b; x.resize(II);

EmpVectorTopo::iterator iter1 = x.begin(); for(iter1; iter1 != x.end(); ++iter1) { if (iter1.index() == 1) { iter1->resize(HH); EmpVectorTopoElem::iterator iter2 = (*iter1).begin(); for(iter2; iter2 != (*iter1).end(); ++iter2) { iter2->resize(GG); } } } x = 0; // x[0][0][0] = 1.0; x[1][1][0] = 21.0; x[1][2][0] = 42.0; b = x;

printvector(std::cout,x,"exact solution","entry",10,10,2); b=0; I.umv(x,b); // set right hand side printvector(std::cout,b,"right hand side","entry",10,10,2); x=0; // initial guess

const int recursion_level = 3; typedef Dune::SeqSSOR<EmpMatrixTopo,EmpVectorTopo, EmpVectorTopo,recursion_level> Preconditioner; typedef Dune::MatrixAdapter<EmpMatrixTopo,EmpVectorTopo,EmpVectorTopo> Operator; typedef Dune::LoopSolver Solver; typedef Dune::SeqScalarProduct ScalarProduct;

Operator op(I); ScalarProduct sp; Preconditioner pr(I,1,1); Solver solver(op, sp, pr, 10e-8, 120, 2); Dune::InverseOperatorResult r; solver.apply(x,b,r);

printvector(std::cout,x,"numerical solution x","entry",10,10,8);

}

int main (int argc , char ** argv) { try { test_Iter(); } catch (Dune::ISTLError& error) { std::cout << error << std::endl; } catch (Dune::Exception& error) { std::cout << error << std::endl; } catch (const std::bad_alloc& e) { std::cout << "memory exhausted" << std::endl; } catch (...) { std::cout << "unknown exception caught" << std::endl; }

return 0; }