Fixed preconditioner bug in serial computations.

AMDiS/src/ITL_OEMSolver.h

@@ -44,19 +44,19 @@ namespace AMDiS {
 
   template <typename ITLSolver>
   class ITL_OEMSolver : public MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector> >  {
-
+    
   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, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >(name)  {}
-
-    int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+    
+      int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
 		    SystemVector& x, 
 		    SystemVector& b,
 		    VectorialMapper& m)
     {
-	    return MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
+      return MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
     }
 
   
@@ -74,7 +74,7 @@ namespace AMDiS {
 	return new ITL_OEMSolver<ITLSolver>(this->name); 
       }
     };
-  
+ 
   };
 
   /**
@@ -91,22 +91,20 @@ namespace AMDiS {
     /// The constructor reads needed parameters and sets solvers \ref name.
     ITL_OEMSolver_para(std::string name) 
       : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >(name) 
-    {
-
-    }
+    {}
 
     int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
 		    SystemVector& x, 
 		    SystemVector& b,
 		    VectorialMapper& m)
     {
-	    return MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
+      return MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
     }
-      
+    
     ~ITL_OEMSolver_para() {}
-
+    
     /// Set parameter of iterative solver
-
+    
     class Creator : public OEMSolverCreator
     {
     public:

AMDiS/src/ITL_OEMSolver.hh

@@ -58,13 +58,13 @@ namespace AMDiS {
       BaseITL_runner(oem_),
       ell(1)
     {
-      ell=oem.getMaxIterations();
+      ell = oem.getMaxIterations();
       Parameters::get(oem.getName() + "->ell", ell);
     }
 
     void init(const Matrix& A)
     {
-	BaseITL_runner::setPrecon(preconPair, A); 
+      BaseITL_runner::setPrecon(preconPair, A); 
     }
 
     int solve(const Matrix& A, Vector& x, Vector& b) 
@@ -84,6 +84,19 @@ namespace AMDiS {
 
       return error;
     }
+
+    void exit()
+    {
+      if (preconPair.l != NULL) {
+	delete preconPair.l;
+	preconPair.l = NULL;
+      }
+
+      if (preconPair.r != NULL) {
+	delete preconPair.r;
+	preconPair.r = NULL;
+      }
+    }
   private:
     PreconPair< Vector > preconPair;
     ITLSolver solver;
@@ -91,35 +104,47 @@ namespace AMDiS {
   };
 
   template< typename ITLSolver, typename Matrix, typename Vector >
-  struct ITL_OEMSolver_runner : BaseITL_runner {
-   
-   
+  struct ITL_OEMSolver_runner : BaseITL_runner {       
     ITL_OEMSolver_runner(OEMSolver* oem):
       BaseITL_runner(oem)
     {}
-
+    
     void init(const Matrix& A)
     {
-	BaseITL_runner::setPrecon(preconPair, A);
+      BaseITL_runner::setPrecon(preconPair, A);
     }
+    
     int solve(const Matrix& A , Vector& x, Vector& b) 
     {
       FUNCNAME("ITL_OEMSolver_runner::solve()");
       TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
       TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
-
+      
       Vector r(A * x - b); 
       itl::cyclic_iteration<typename Matrix::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(), 
 							      oem.getTolerance(), oem.getPrint_cycle());
-
+      
       int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter);
       oem.setErrorCode(error);
       oem.setIterations(iter.iterations());
       oem.setResidual(iter.resid());
-
+      
       return error;
     }
-    private:
+
+    void exit()
+    {
+      if (preconPair.l != NULL) {
+	delete preconPair.l;
+	preconPair.l = NULL;
+      }
+
+      if (preconPair.r != NULL) {
+	delete preconPair.r;
+	preconPair.r = NULL;
+      }
+    }
+  private:
     ITLSolver solver;
     PreconPair< Vector > preconPair;
   };

AMDiS/src/ITL_Preconditioner.h

@@ -34,10 +34,6 @@
 
 namespace AMDiS {
 
-  // ============================================================================
-  // ===== class ITL_Preconditioner ==============================================
-  // ============================================================================
-
  
   /**
    * \ingroup Solver
@@ -48,26 +44,25 @@ namespace AMDiS {
   class ITL_BasePreconditioner 
   {
   public:
-      virtual ~ITL_BasePreconditioner() {}
-
-      // Old style AMDiS
-      // virtual Vec member_solve(const Vec& vin) const =0 ;
-      // virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;
-
-     // New style without copying the result vector
-     virtual void solve(const Vec& x, Vec& y) const= 0;
-     virtual void adjoint_solve(const Vec& x, Vec& y) const= 0;
+    virtual ~ITL_BasePreconditioner() {}
+    
+    // Old style AMDiS
+    // virtual Vec member_solve(const Vec& vin) const =0 ;
+    // virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;
+    
+    // New style without copying the result vector
+    virtual void solve(const Vec& x, Vec& y) const= 0;
+    virtual void adjoint_solve(const Vec& x, Vec& y) const= 0;
   };
-
-
-
+  
+  
+  
   template< typename Vec >
   inline //Vec
   itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>
   solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
   {
     return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>(P, vin);
-    //return P.member_solve(vin);
   }
 
   template< typename Vec >
@@ -76,7 +71,6 @@ namespace AMDiS {
   adjoint_solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
   {
     return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, true>(P, vin);
-    //return P.member_adjoint_solve(vin);
   }
 
   /**
@@ -89,31 +83,16 @@ namespace AMDiS {
   {
   public:
     ITL_Preconditioner(const MTLMatrix& A) : precon(A) {}
-#if 0   
-    MTLVector
-    member_solve(const MTLVector& vin) const
-    {
-      return solve(precon, vin);
-    }
-
-    MTLVector
-    member_adjoint_solve(const MTLVector& vin) const
-    {
-      return adjoint_solve(precon, vin);
-    }
-#else
-      
     void solve(const MTLVector& vin, MTLVector& vout) const
     {
       precon.solve(vin, vout);
-   }
-
+    }
+    
     void adjoint_solve(const MTLVector& vin, MTLVector& vout) const
     {
       precon.adjoint_solve(vin, vout);
-   }
-#endif
-
+    }
+    
     /// Creator class used in the OEMSolverMap.
     class Creator : public CreatorInterface< ITL_BasePreconditioner< MTLVector > >
     {
@@ -127,16 +106,13 @@ namespace AMDiS {
 	return new ITL_Preconditioner<Preconditioner, MTLVector, MTLMatrix>(A);
       }
     };
-
+    
   private:
-      Preconditioner precon;
+    Preconditioner precon;
   };
+  
 
 
-  // ============================================================================
-  // ===== class DiagonalPreconditioner =========================================
-  // ============================================================================
-
   /**
    * \ingroup Solver
    * 
@@ -168,26 +144,22 @@ namespace AMDiS {
   };
 
 
-  // ============================================================================
-  // ===== class ILUPreconditioner ==============================================
-  // ============================================================================
-
   /**
    * \ingroup Solver
    * 
    * \brief
    * ILU (Incomplete LU factorization) preconditioner.
    *
-   * The preconditioner is used from ITL. It corresponds for instance to "Iterative Methods for
-   * Sparce Linear Systems", second edition, Yousef Saad. The preconditioner is
-   * described in chapter 10.3 (algorithm 10.4).
+   * The preconditioner is used from ITL. It corresponds for instance to
+   * "Iterative Methods for Sparce Linear Systems", second edition, Yousef Saad.
+   *  The preconditioner is described in chapter 10.3 (algorithm 10.4).
    */
   class ILUPreconditioner 
-      : public ITL_Preconditioner< itl::pc::ilu_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
+    : public ITL_Preconditioner< itl::pc::ilu_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
   {};
-
+  
   class ICPreconditioner 
-      : public ITL_Preconditioner< itl::pc::ic_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
+    : public ITL_Preconditioner< itl::pc::ic_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
   {};
 
 

AMDiS/src/MTL4Solver.h

@@ -28,7 +28,7 @@
 #include <boost/numeric/mtl/utility/is_distributed.hpp>
 
 #if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
-# include <boost/numeric/mtl/par/distribution.hpp>
+#include <boost/numeric/mtl/par/distribution.hpp>
 #endif
 
 namespace AMDiS {
@@ -71,7 +71,7 @@ namespace AMDiS {
       FUNCNAME("MTL4Solver::solve()");
 
       Timer t;
-      if(num_rows(matrix) == 0 || !getMultipleRhs() ) {
+      if (num_rows(matrix) == 0 || !getMultipleRhs()) {
 	init(mapper, mtl::traits::is_distributed<MTLMatrix>());
 
         set_to_zero(matrix);
@@ -81,13 +81,13 @@ namespace AMDiS {
       }
 
       MTLVector mtl_x;
-      init(mtl_x,mtl::traits::is_distributed<MTLVector>());
+      init(mtl_x, mtl::traits::is_distributed<MTLVector>());
       set_to_zero(mtl_x);
       VecMap< Vector, Mapper > xVecMap(x, mapper);
       mtl_x << xVecMap;
 
       MTLVector mtl_b;
-      init(mtl_b,mtl::traits::is_distributed<MTLVector>());
+      init(mtl_b, mtl::traits::is_distributed<MTLVector>());
       set_to_zero(mtl_b);
       VecMap< Vector, Mapper> bVecMap(b, mapper);
       mtl_b << bVecMap;
@@ -104,6 +104,9 @@ namespace AMDiS {
       MSG("MTL4Solver: ||b-Ax||= %e\n", residual);
 
       mtl_x >> xVecMap;
+
+      worker.exit();
+
       return error;
     }
 

AMDiS/src/OEMSolver.h

@@ -45,7 +45,6 @@
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
 #include "parallel/ParallelMapper.h"
 #endif
-//#include "ITL_Preconditioner.h"
 
 namespace AMDiS {
 
@@ -99,7 +98,8 @@ namespace AMDiS {
 			    VectorialMapper&) 
     {
       FUNCNAME("OEMSolver::solveSystem()");
-      TEST_EXIT(false)("This linear solver is not suitable for sequentiell problems\n");
+      TEST_EXIT(false)
+	("This linear solver is not suitable for sequentiell problems\n");
       return -1;
     }
 
@@ -118,7 +118,8 @@ namespace AMDiS {
 			    ParallelMapper&) 
     { 
       FUNCNAME("OEMSolver::solveSystem()");
-      TEST_EXIT(false)("This linear solver is not suitable for domaindecomposition problems\n");
+      TEST_EXIT(false)
+	("This linear solver is not suitable for domaindecomposition problems\n");
       return -1;
     }
 #endif
@@ -269,11 +270,9 @@ namespace AMDiS {
     /// Error code  in last solver (not set by UmfPack)
     int error;
 
-    /** \brief
-     * If true, the solver has to solve multiple right hand sides with the same
-     * matrix. Some solvers, e.g. direct once, may take advanteges from this knowledge,
-     * as they can do the factorization of the matrix only once.
-     */
+    /// If true, the solver has to solve multiple right hand sides with the same
+    /// matrix. Some solvers, e.g. direct once, may take advanteges from this knowledge,
+    /// as they can do the factorization of the matrix only once.
     bool multipleRhs;
   };
 
@@ -295,10 +294,8 @@ namespace AMDiS {
     }
 
   protected:
-    /** \brief
-     * Pointer to the problem the solver will belong to. Needed as parameter
-     * when constructing an OEMSolver
-     */
+    /// Pointer to the problem the solver will belong to. Needed as parameter
+    /// when constructing an OEMSolver
     std::string name;
   };
 }

AMDiS/src/UmfPackSolver.h

@@ -94,9 +94,11 @@ namespace AMDiS {
 	TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	
       }
       return code;
-
     }
 
+    void exit()
+    {}
+
     ~Umfpack_runner() 
     {
       if (solver != NULL)
@@ -135,17 +137,15 @@ namespace AMDiS {
     /// Constructor
     UmfPackSolver(std::string name) 
       : MTL4Solver< MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >(name)
-    {
-    }
-
+    {}
+    
     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);
+    }    
     
   private:
   };