[Dune] Nested BCRSMatrix problems

[Patrick Leidenberger]

Hi Dune,

using the BCRSMatrix in a nested way shows me some problems:

I create a 4x4 BCRSMatrix of BCRSMatrices. Only in the second line and
column I insert an object, lets say a 5x5 BCRSMatrix. If I look now at
the construct, I get rowdim 10 and coldim 5: In the first row, second
column exists an empty 5x5 Matrix, which I have not created!
This problem can be avoid, if the empty diagonal elements in the 4x4
matrix are initialized with empty BCRSMatrices. If all elements in the
4x4 matrix are initialized with such empty matrices the total dimension
of the matrix is 0x0!
Perhaps this is a bug, if it isn't this behavior should be documented.
I modified the dune-istl/istl/tutorial/example.cc to show the problem.
The file is attached, this is the output:
============================================.
E [n=5,m=5,rowdim=5,coldim=5]
row    0   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    1   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    2   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    3   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    4   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
=============================================.
============================================.
D [n=4,m=4,rowdim=10,coldim=5]
row    0          .          .          .          .          .
row    1          .          .          .          .          .
row    2          .          .          .          .          .
row    3          .          .          .          .          .
row    4          .          .          .          .          .
row    5   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    6   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    7   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    8   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
row    9   4.20e+01   4.20e+01   4.20e+01   4.20e+01   4.20e+01
=============================================.


The printmatrix function for this nested BCRSMatrices works fine. But is
it possible to extend the writeMatrixToMatlab function so that it will
works as well?

Regards Patrick
-------------- next part --------------
// 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>

// a simple stop watch
class Timer
{
public:
  Timer ()
  {
	struct tms buf;
	cstart = times(&buf);
  }

  void start ()
  {
	struct tms buf;
	cstart = times(&buf);
  }

  double stop ()
  {
	struct tms buf;
	cend = times(&buf);
	return ((double)(cend-cstart))/100.0;
  }

  double gettime ()
  {
	return ((double)(cend-cstart))/100.0;
  }

private:
  clock_t cstart,cend;  
};


void test_BCRSMatrix ()
{
  const int N=5, M =4,M2 =5, L = 4;
  typedef double ValueType;
  typedef Dune::BCRSMatrix<Dune::FieldMatrix<ValueType,1,1> > EmpMatrixTopoElemOrder;
  typedef Dune::BCRSMatrix<EmpMatrixTopoElemOrder>            EmpMatrixTopoElem;


  Dune::FieldMatrix<ValueType,1,1> F;
  F = 42.0;


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

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

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


  EmpMatrixTopoElem D(M,M,2,EmpMatrixTopoElem::row_wise);
  for (EmpMatrixTopoElem::CreateIterator i=D.createbegin(); i!=D.createend(); ++i)
    for (int j=0; j<M; ++j)
    {
      if ((i.index() == 1) && (j == 1))
      {
        i.insert(j);
      }
    }
    
  for (EmpMatrixTopoElem::RowIterator i=D.begin(); i!=D.end(); ++i)
  for (EmpMatrixTopoElem::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
    *j = E;

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


}


int main (int argc , char ** argv)
{
  try {
 	test_BCRSMatrix();
  }
  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;
}