#81 CGSolver crashes on doubly nested BCRSMatrix

Metadata

Property	Value
Reported by	Oliver Sander (oliver.sander@tu-dresden.de)
Reported at	Jan 19, 2006 09:53
Type	Bug Report
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Last edited by	Markus Blatt (markus@dr-blatt.de)
Last edited at	Jan 23, 2006 15:33
Closed by	Markus Blatt (markus@dr-blatt.de)
Closed at	Aug 11, 2006 13:00
Closed in version	1.0
Resolution	Fixed
Comment	This is fixed. Concerning the testing of the ISTL kernels there bis another task available

Description

The ISTL CGSolver segfaults when operating on a doubly nested BCRSMatrix. The reason seems to be an incomplete resize of an internal variable.

Here's a sample program which tries to solve a system with the identity matrix:

#include <config.h>

#include <dune/common/matrix.hh> #include <dune/istl/solvers.hh>

using namespace Dune;

int main (int argc, char *argv[]) { typedef BlockVector<BlockVector<FieldVector<double,1> > > VectorType; typedef BCRSMatrix<BCRSMatrix<FieldMatrix<double,1,1> > > MatrixType; typedef BCRSMatrix<FieldMatrix<double,1,1> > MatrixEntryType;

VectorType rhs(2);
VectorType x(2);

rhs[0].resize(2);
rhs[1].resize(2);
x[0].resize(2);
x[1].resize(2);

// Outer matrix
MatrixType stiffnessMatrix(2,2, MatrixType::random);
stiffnessMatrix.setrowsize(0,1);
stiffnessMatrix.setrowsize(1,1);
stiffnessMatrix.endrowsizes();

stiffnessMatrix.addindex(0,0);
stiffnessMatrix.addindex(1,1);
stiffnessMatrix.endindices();

stiffnessMatrix = 0;

// Make 2x2 identity matrix for the diagonal blocks
MatrixEntryType identity(2,2, MatrixEntryType::random);
identity.setrowsize(0,1);
identity.setrowsize(1,1);
identity.endrowsizes();

identity.addindex(0,0);
identity.addindex(1,1);

identity.endindices();

identity[0][0] = 1;
identity[1][1] = 1;

// ///////////////////////////////////////////////
//   Set diagonal of BlockMatrix to Identity
// ///////////////////////////////////////////////

stiffnessMatrix[0][0] = identity;
stiffnessMatrix[1][1] = identity;
    
x = 0;
rhs = 1;

// /////////////////////////
//   Compute solution
// /////////////////////////

MatrixAdapter<MatrixType,VectorType,VectorType> op(stiffnessMatrix);
SeqSSOR<MatrixType,VectorType,VectorType,2> ssor(stiffnessMatrix,1,1.0);
CGSolver<VectorType> cg(op,ssor,1E-4,50,2); 
InverseOperatorResult statistics;

// Solve!
cg.apply(x, rhs, statistics);

}