Word-around for a bug in solving Rosenbrock system using UMFPACK.

AMDiS/src/DOFMatrix.cc

@@ -241,7 +241,8 @@ namespace AMDiS {
       } else {
 	for (int j = 0; j < nCol; j++) {
 	  DegreeOfFreedom col = colIndices[j];
-	  ins[row][col] += elMat[i][j];
+	  if (fabs(elMat[i][j]) > 1e-10)
+	    ins[row][col] += elMat[i][j];
 	}
       }
     }

AMDiS/src/Debug.cc

@@ -437,6 +437,7 @@ namespace AMDiS {
 
       std::ofstream out;
       out.open(filename.c_str());
+      out.precision(20);
     
       for (cursor_type cursor = begin<major>(mat.getBaseMatrix()), 
 	     cend = end<major>(mat.getBaseMatrix()); cursor != cend; ++cursor)

AMDiS/src/ITL_OEMSolver.h

@@ -45,11 +45,11 @@ namespace AMDiS {
   template <typename ITLSolver>
   class ITL_OEMSolver : public MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > >  {
 
-    protected:
-    public:
+  protected:
+  public:
     /// The constructor reads needed parameters and sets solvers \ref name.
     ITL_OEMSolver(std::string name) : 
-	    MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > >(name)  {}
+      MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > >(name)  {}
 
   
     ~ITL_OEMSolver() {}
@@ -77,19 +77,19 @@ namespace AMDiS {
    */
   template <typename ITLSolver>
   class ITL_OEMSolver_para
-      : public MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > // ITL_OEMSolver<ITLSolver>
+    : public MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > // ITL_OEMSolver<ITLSolver>
   {
   public:
-      /// The constructor reads needed parameters and sets solvers \ref name.
-      ITL_OEMSolver_para(std::string name) 
-	  : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name) 
-      {
+    /// The constructor reads needed parameters and sets solvers \ref name.
+    ITL_OEMSolver_para(std::string name) 
+      : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name) 
+    {
 
-      }
+    }
       
-      ~ITL_OEMSolver_para() {}
+    ~ITL_OEMSolver_para() {}
 
-      /// Set parameter of iterative solver
+    /// Set parameter of iterative solver
 
     class Creator : public OEMSolverCreator
     {

AMDiS/src/ITL_OEMSolver.hh

@@ -9,51 +9,51 @@ namespace AMDiS {
     /// Pointer to the right preconditioner
     ITL_BasePreconditioner< Vector >* r;
     PreconPair():
-	    l(NULL), r(NULL)
-	  {}
+      l(NULL), r(NULL)
+    {}
 
   };
 
   struct BaseITL_runner { 
     OEMSolver& oem;
     BaseITL_runner(OEMSolver* oemPtr):
-	    oem(*oemPtr)
+      oem(*oemPtr)
     {
       TEST_EXIT(oemPtr != NULL)("need a real oemsolver\n");
     }
 
     template< typename Vector , typename Matrix >
-      void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
-      {
-	/// Creator for the left preconditioner
-	CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);
-
-	/// Creator for the right preconditioner
-	CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);
-
-	std::string preconType("no");
-	Parameters::get(oem.getName() + "->left precon", preconType);
-
-	leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
-	TEST_EXIT(leftCreator != NULL)
-		("there is no creator for the given left preconditioner");
-
-	preconType = "no";
-	Parameters::get(oem.getName() + "->right precon", preconType);
-	rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
-	TEST_EXIT(rightCreator != NULL)
-		("there is no creator for the given right preconditioner");
-
-	pair.l = leftCreator->create(matrix);
-	pair.r = rightCreator->create(matrix);
-      }
+    void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
+    {
+      /// Creator for the left preconditioner
+      CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);
+
+      /// Creator for the right preconditioner
+      CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);
+
+      std::string preconType("no");
+      Parameters::get(oem.getName() + "->left precon", preconType);
+
+      leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+      TEST_EXIT(leftCreator != NULL)
+	("there is no creator for the given left preconditioner");
+
+      preconType = "no";
+      Parameters::get(oem.getName() + "->right precon", preconType);
+      rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+      TEST_EXIT(rightCreator != NULL)
+	("there is no creator for the given right preconditioner");
+
+      pair.l = leftCreator->create(matrix);
+      pair.r = rightCreator->create(matrix);
+    }
   };
 
   template< typename ITLSolver >
   struct ITL_OEMSolver_para_runner : BaseITL_runner {
     ITL_OEMSolver_para_runner(OEMSolver* oem_):
-	      BaseITL_runner(oem_),
-	      ell(1)
+      BaseITL_runner(oem_),
+      ell(1)
     {
       ell=oem.getMaxIterations();
       Parameters::get(oem.getName() + "->ell", ell);
@@ -63,11 +63,11 @@ namespace AMDiS {
     int solve(const Matrix& A, Vector& x, Vector& b) 
     {
       itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(), 
-					     oem.getTolerance(), oem.getPrint_cycle());
+							      oem.getTolerance(), oem.getPrint_cycle());
 
       static PreconPair< Vector > preconPair;
       if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
-	     BaseITL_runner::setPrecon(preconPair, A); 
+	BaseITL_runner::setPrecon(preconPair, A); 
 
       int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter, ell);
       oem.setIterations(iter.iterations());
@@ -75,7 +75,7 @@ namespace AMDiS {
       oem.setErrorCode(error);
       return error;
     }
-    private:
+  private:
     int ell;
   };
 
@@ -84,18 +84,18 @@ namespace AMDiS {
    
    
     ITL_OEMSolver_runner(OEMSolver* oem):
-	    BaseITL_runner(oem)
+      BaseITL_runner(oem)
     {}
 
     template< typename Matrix, typename Vector>
     int solve(const Matrix& A, Vector& x, Vector& b) 
     {
       itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(), 
-					     oem.getTolerance(), oem.getPrint_cycle());
+							      oem.getTolerance(), oem.getPrint_cycle());
 
       static PreconPair< Vector > preconPair;
       if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
-	     BaseITL_runner::setPrecon(preconPair, A); 
+	BaseITL_runner::setPrecon(preconPair, A); 
 
       int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter);
       oem.setErrorCode(error);

AMDiS/src/MTL4Solver.h

@@ -25,53 +25,49 @@
 #include <iostream>
 
 namespace AMDiS {
-  template< typename MTLMatrix, typename MTLVector, typename Worker >
-  class MTL4Solver : public OEMSolver {
-
-    protected:
 
-      Worker worker;
-
-      template< typename Matrix, typename Vector, typename Mapper >
-      int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper) 
-      {
-	 MTLMatrix matrix(mapper.nRow(), mapper.nCol());
-	 set_to_zero(matrix);
-	 MatMap< const Matrix, Mapper > matMap(A,mapper);
-	 matrix << matMap;
+  template< typename MTLMatrix, typename MTLVector, typename Worker >
+  class MTL4Solver : public OEMSolver {    
+  protected:    
+    Worker worker;
 
-         MTLVector mtl_x(mapper.nRow());
-	 set_to_zero(mtl_x);
-	 VecMap< Vector, Mapper > xVecMap(x, mapper);
-	 mtl_x << xVecMap;
+    template< typename Matrix, typename Vector, typename Mapper >
+    int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper) 
+    {
+      MTLMatrix matrix(mapper.nRow(), mapper.nCol());
+      set_to_zero(matrix);
+      MatMap< const Matrix, Mapper > matMap(A,mapper);
+      matrix << matMap;
 
-	 MTLVector mtl_b(mapper.nRow());
-	 set_to_zero(mtl_b);
-	 VecMap< Vector, Mapper> bVecMap(b, mapper);
-	 mtl_b << bVecMap;
+      MTLVector mtl_x(mapper.nRow());
+      set_to_zero(mtl_x);
+      VecMap< Vector, Mapper > xVecMap(x, mapper);
+      mtl_x << xVecMap;
 
-	 error = worker.solve(matrix, mtl_x, mtl_b);
+      MTLVector mtl_b(mapper.nRow());
+      set_to_zero(mtl_b);
+      VecMap< Vector, Mapper> bVecMap(b, mapper);
+      mtl_b << bVecMap;
 
-	 mtl_x >> xVecMap;
-	 return error;
-      }
+      error = worker.solve(matrix, mtl_x, mtl_b);
 
-    public:
-      MTL4Solver(std::string n):
-	OEMSolver(n),
-        worker(this)
-      {
-	FUNCNAME("MTL4Solver()");
+      mtl_x >> xVecMap;
+      return error;
+    }
 
-      }
+  public:
+    MTL4Solver(std::string n):
+      OEMSolver(n),
+      worker(this)
+    {}
 
-      virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+    virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
 			    SystemVector& x, 
 			    SystemVector& b,
 			    VectorialMapper& m)
-      {
-	return solve(A,x,b,m);
-      }
+    {
+      return solve(A,x,b,m);
+    }
 
       
   };

AMDiS/src/OEMSolver.h

@@ -73,9 +73,7 @@ namespace AMDiS {
 	print_cycle(100),
 	iterations(-1),
 	error(-1),
-	multipleRhs(false)/*,
-	leftPrecon(NULL),
-	rightPrecon(NULL)*/
+	multipleRhs(false)
     {}
 
     ///
@@ -99,47 +97,13 @@ namespace AMDiS {
 			    SystemVector& x, 
 			    SystemVector& b,
 			    VectorialMapper&) = 0;
-#if 0
-    {
-      int ns = x.getSize(),  // Number of systems.
-	size = x.getUsedSize(); // Size of all DOFVectors
-
-      // Copy vectors
-      mtl::dense_vector<value_type>   xx(size), bb(size);
-      xx << x;
-      bb << b
-/*      for (int i = 0, counter = 0; i < ns; i++) {
-	DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);	  
-	DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
-	
-	for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {	
-	  bb[counter] = *it_b;
-	  xx[counter] = *it_x;
-	  counter++;
-	}
-      }
-     */
-      // Solver on DOFVector for single system
-      DOFMatrix::base_matrix_type baseA;
-      baseA << A;
-      int r = solveSystem(baseA, xx, bb);
-      
-      for (int i = 0, counter = 0; i < ns; i++) {
-	DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);	  
-	for (it.reset(); !it.end(); ++it)	
-	  *it = xx[counter++];	
-      }
-      
-      return r;
-    }
-#endif
 
     inline int solveSystem(const SolverMatrix< Matrix< DOFMatrix* > >& A,
 		    SystemVector& x,
 		    SystemVector& b) 
     {
-	    VectorialMapper mapper(A);
-	    return solveSystem(A, x, b, mapper);
+      VectorialMapper mapper(A);
+      return solveSystem(A, x, b, mapper);
     }
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -306,13 +270,6 @@ namespace AMDiS {
      * as they can do the factorization of the matrix only once.
      */
     bool multipleRhs;
-#if 0
-    /// Left preconditioner (ignored by some iterative solvers and by UmfPack)
-    PreconCreator *leftPrecon;
-
-    /// Right preconditioner (ignored by some iterative solvers and by UmfPack)
-    PreconCreator *rightPrecon;
-#endif
   };
 
   /** 

AMDiS/src/ProblemStat.cc

@@ -634,6 +634,8 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemStat::buildAfterCoarsen()");
 
+    VtkWriter::writeFile(rhs->getDOFVector(0), "test.vtu");
+
 //     buildAfterCoarsen_sebastianMode(adaptInfo, flag);
 //     return;
 

AMDiS/src/UmfPackSolver.h

@@ -52,6 +52,8 @@ namespace AMDiS {
       store_symbolic(0),
       symmetric_strategy(0)
     {
+      FUNCNAME("Umfpack_runner::Umfpack_runner()");
+
       TEST_EXIT_DBG(oem_ != NULL)("Need real OEMSolver\n"); 
       Parameters::get(oem.getName() + "->store symbolic", store_symbolic);
       Parameters::get(oem.getName() + "->symmetric strategy", symmetric_strategy);
@@ -60,15 +62,14 @@ namespace AMDiS {
     template< typename Vector> 
     int solve(const Matrix& A, Vector& x, Vector& b) 
     { 
-     FUNCNAME("(UmfPackSolver.h)::solve()");
-     if (!solver) {
+      FUNCNAME("Umfpack_runner::solve()");
+      if (!solver) {
         if (!symmetric_strategy)
 	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
 	else
 	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
       } else {
 	if (!oem.getMultipleRhs()) {
-
  	  if (store_symbolic)
  	    solver->update_numeric();
  	  else
@@ -77,16 +78,14 @@ namespace AMDiS {
       }
 
       int code = (*solver)(x, b);
-      if( oem.getInfo() > 0) {
+      if (oem.getInfo() > 0) {
         mtl::dense_vector<typename matrix_type::value_type> r(b); 
         r -= A * x; 
         double residual = two_norm(r);
         oem.setResidual(residual);
         oem.setErrorCode(code);
         std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
-	if (residual > oem.getTolerance()) {
-	    ERROR_EXIT("Tolerance tol=%e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	  
-	}
+	TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	
       }
       return code;
 
@@ -94,11 +93,11 @@ namespace AMDiS {
 
     ~Umfpack_runner() 
     {
-      if(!solver)
+      if (!solver)
 	delete solver;
     }
 
-    private:
+  private:
     mtl::matrix::umfpack::solver<matrix_type> *solver;
 
     int store_symbolic;
@@ -108,7 +107,7 @@ namespace AMDiS {
   using namespace MTLTypes;
   class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
   {
-    protected:
+  protected:
 
   public:
     /// Creator class used in the OEMSolverMap.

AMDiS/src/time/RosenbrockStationary.cc

@@ -14,6 +14,7 @@
 #include "ProblemStat.h"
 #include "SystemVector.h"
 #include "OEMSolver.h"
+#include "Debug.h"
 
 namespace AMDiS {
 
@@ -79,7 +80,7 @@ namespace AMDiS {
       }
       
       ProblemStatSeq::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
-      solver->setMultipleRhs(i != 0);
+      //      solver->setMultipleRhs(i != 0);
       ProblemStatSeq::solve(adaptInfo);
       
       *(stageSolutions[i]) = *solution;