[Dune] Notification from Dune
Flyspray
dune at hal.iwr.uni-heidelberg.de
Wed Dec 13 14:27:02 CET 2006
THIS IS AN AUTOMATED MESSAGE, DO NOT REPLY.
A new Flyspray task has been opened. Details are below.
User who did this: - Patrick Leidenberger (leidenberger)
Attached to Project - Dune
Summary - ISTL solver causes segfault for special hierarchical matrix
Task Type - Bug Report
Category - ISTL
Status - Unconfirmed
Assigned To -
Operating System - All
Severity - Low
Priority - Normal
Reported Version - SVN
Due in Version -
Due Date - Undecided
Details - 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 at 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<iostream> // for input/output to shell
#include<fstream> // for input/output to files
#include<vector> // STL vector class
#include<complex>
#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<EmpMatrixTopoElemOrder>
EmpMatrixTopoElem;
typedef Dune::BCRSMatrix<EmpMatrixTopoElem>
EmpMatrixTopo;
typedef Dune::BlockVector<Dune::FieldVector<ValueType,1> >
EmpVectorTopoElemOrder;
typedef Dune::BlockVector<EmpVectorTopoElemOrder>
EmpVectorTopoElem;
typedef Dune::BlockVector<EmpVectorTopoElem> 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<EmpVectorTopo> Solver;
typedef Dune::SeqScalarProduct<EmpVectorTopo> 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;
}
More information can be found at the following URL:
http://hal.iwr.uni-heidelberg.de/flyspray/?do=details&id=222
You are receiving this message because you have requested it from the
Flyspray bugtracking system. You can be removed from future
notifications by visiting the URL shown above.
More information about the Dune
mailing list