#62 Impossible to add sparser to denser Matrix

Metadata

Property	Value
Reported by	Oliver Sander (oliver.sander@tu-dresden.de)
Reported at	Dec 9, 2005 09:59
Type	Bug Report
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Closed by	Markus Blatt (markus@dr-blatt.de)
Closed at	Dec 9, 2005 11:41
Closed in version	Unknown
Resolution	Fixed
Comment	Yes this was a bug.

As yu told the prerequsite of the operation was wrong.

It should be all indices of the row to be added should be present in the row it is added to.

Fixed that.

Now operator += and operator-= work as excepcted.

in operator* (the euclidian scalar product) it is checked whether the two index sets match. |

Description

I would like to be able to add a BCRSMatrix foo to another one bar if I know that bar contains all the entries that foo has. First of all, there is no operator+= for BCRSMatrix, but I can loop over the matrix rows and than there is an operator+= for matrix lines. However, when writing

bar[i] += foo[i]

and exception is thrown if bar[i] contains an entry which foo[i] doesn't contain. Strangely enough, it works the other way around, even though in that case you lose data.

Here is a full test program:

#include "config.h"

#include <dune/common/fmatrix.hh> #include <dune/istl/bcrsmatrix.hh>

using namespace Dune;

int main() { typedef BCRSMatrix<FieldMatrix<double,1,1> > MatrixType;

// Make dense 2x2 matrix
MatrixType fullMatrix(2,2,MatrixType::random);

fullMatrix.setrowsize(0,2);
fullMatrix.setrowsize(1,2);
fullMatrix.endrowsizes();

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

fullMatrix.endindices();

// Make diagonal 2x2 matrix
MatrixType diagonalMatrix(2,2,MatrixType::random);

diagonalMatrix.setrowsize(0,1);
diagonalMatrix.setrowsize(1,1);
diagonalMatrix.endrowsizes();

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

diagonalMatrix.endindices();

// //////////////////////////////////////////////////////
//   What I really want is fullMatrix += diagonalMatrix;
// //////////////////////////////////////////////////////

// This works, but information is lost, since fullMatrix
// contains more entries than diagonalMatrix
for (int i=0; i<2; i++)
    diagonalMatrix[i] += fullMatrix[i];

// fullMatrix contains all entries of diagonalMatrix,
// yet an exception is thrown!
for (int i=0; i<2; i++)
    fullMatrix[i] += diagonalMatrix[i];

}