corrected bug in block tridiag solver - mistake in left/right matrixmultiplication.

To test the implementation improved the test to not only include diagonal block matrices on the
diagonal, because with diagonal left/right multiplication does not make a differece.

dune/istl/btdmatrix.hh

@@ -124,8 +124,8 @@ namespace Dune {
         c[i] = (*this)[i][i+1];
 
       /* Modify the coefficients. */
-      block_type a_00_inv;
-      FMatrixHelp::invertMatrix((*this)[0][0], a_00_inv);
+      block_type a_00_inv = (*this)[0][0];
+      a_00_inv.invert();
 
       //c[0] /= (*this)[0][0];	/* Division by zero risk. */
       block_type c_0_tmp = c[0];
@@ -133,13 +133,13 @@ namespace Dune {
 
       // d = a^{-1} d        /* Division by zero would imply a singular matrix. */
       typename V::block_type d_0_tmp = d[0];
-      (*this)[0][0].solve(d[0], d_0_tmp);
+      a_00_inv.mv(d_0_tmp,d[0]);
 
       for (unsigned int i = 1; i < this->N(); i++) {
 
         // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
         block_type tmp;
-        FMatrixHelp::multMatrix(c[i-1], (*this)[i][i-1], tmp);
+        FMatrixHelp::multMatrix((*this)[i][i-1],c[i-1], tmp);
         block_type id = (*this)[i][i];
         id -= tmp;
         id.invert();         /* Division by zero risk. */
@@ -147,7 +147,7 @@ namespace Dune {
         if (i<c.size()) {
           // c[i] *= id
           tmp = c[i];
-          FMatrixHelp::multMatrix(tmp, id, c[i]);                      /* Last value calculated is redundant. */
+          FMatrixHelp::multMatrix(id,tmp, c[i]);                      /* Last value calculated is redundant. */
         }
 
         // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;

dune/istl/test/matrixtest.cc

@@ -287,6 +287,7 @@ int main()
   feenableexcept(FE_INVALID);
 #endif
 
+
   // ////////////////////////////////////////////////////////////
   //   Test the Matrix class -- a scalar dense dynamic matrix
   // ////////////////////////////////////////////////////////////
@@ -397,7 +398,42 @@ int main()
     btdMatrix[i][i] = ScaledIdentityMatrix<double,2>(1+i);
 
   for (size_type i=0; i<btdMatrix.N()-1; i++)
-    btdMatrix[i][i+1] = ScaledIdentityMatrix<double,2>(2+i);               // first off-diagonal
+    btdMatrix[i][i+1] = ScaledIdentityMatrix<double,2>(2+i);               // upper off-diagonal
+  for (size_type i=1; i<btdMatrix.N(); i++)
+    btdMatrix[i-1][i] = ScaledIdentityMatrix<double,2>(2+i);               // lower off-diagonal
+
+  // add some off diagonal stuff to the blocks in the matrix
+  // diagonals
+  btdMatrix[0][0][0][1] = 2;
+  btdMatrix[0][0][1][0] = -1;
+
+  btdMatrix[1][1][0][1] = 2;
+  btdMatrix[1][1][1][0] = 3;
+
+  btdMatrix[2][2][0][1] = 2;
+  btdMatrix[2][2][0][0] += sqrt(2.);
+  btdMatrix[2][2][1][0] = 3;
+
+  btdMatrix[3][3][0][1] = -1;
+  btdMatrix[3][3][0][0] -= 0.5;
+  btdMatrix[3][3][1][0] = 2;
+
+  // off diagonals
+  btdMatrix[0][1][0][1] = std::sqrt(2);
+  btdMatrix[1][0][0][1] = std::sqrt(2);
+
+  btdMatrix[1][0][1][0] = -13./17.;
+  btdMatrix[1][2][0][1] = -1./std::sqrt(2);
+  btdMatrix[1][2][1][0] = -13./17.;
+
+  btdMatrix[2][1][0][1] = -13./17.;
+  btdMatrix[2][1][1][0] = -13./17.;
+  btdMatrix[2][3][0][1] = -1./std::sqrt(2);
+  btdMatrix[2][3][1][0] = -17.;
+
+  btdMatrix[3][2][0][1] = 1.;
+  btdMatrix[3][2][1][0] = 1.;
+
 
   BTDMatrix<FieldMatrix<double,2,2> > btdMatrixThrowAway = btdMatrix;    // the test method overwrites the matrix
   testSuperMatrix(btdMatrixThrowAway);
@@ -460,5 +496,4 @@ int main()
   sIdMatrix = 3.1459;
 
   testMatrix(sIdMatrix, fvX, fvY);
-
 }