]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Refactor AMG preconditioner and add AMG solver.
authorturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Jun 2014 20:49:24 +0000 (20:49 +0000)
committerturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 26 Jun 2014 20:49:24 +0000 (20:49 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_paralution@33091 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/block_vector.h
deal.II/include/deal.II/lac/paralution_precondition.h
deal.II/include/deal.II/lac/paralution_solver.h
deal.II/include/deal.II/lac/paralution_vector.h
deal.II/source/lac/paralution_precondition.cc
deal.II/source/lac/paralution_solver.cc
deal.II/source/lac/paralution_vector.cc

index 77108b55258da4980eb1554d83cffa98b8ea8eba..64a19b44cc629346c9c3058a88f2575da598747c 100644 (file)
@@ -112,7 +112,7 @@ public:
   /**
    * Copy-Constructor. Dimension set to
    * that of V, all components are copied
-   * from V
+   * from V.
    */
   BlockVector (const BlockVector<Number> &V);
 
@@ -195,7 +195,7 @@ public:
                const InputIterator              end);
 
   /**
-   * Destructor. Clears memory
+   * Destructor. Clears memory.
    */
   ~BlockVector ();
 
index f6ffe9e21c149e73012c772684e4219290bfb346..4e94a075cc069047ea1fb740776431a8a3739e07 100644 (file)
@@ -35,7 +35,7 @@ DEAL_II_NAMESPACE_OPEN
 namespace ParalutionWrappers
 {
   /**
-   * Cycle for multigrid preconditioner.
+   * Cycle for multigrid.
    */
   enum mg_cycle {V_cycle, W_cycle, K_cycle, F_cycle};
 
@@ -45,14 +45,13 @@ namespace ParalutionWrappers
   enum mg_interpolation {aggregation, smoothed_aggregation};
 
   /**
-   * Precondition for the smoother.
+   * Preconditioner for the smoother of the multigrid.
    */
   enum mg_preconditioner {jacobi, sgs, multicolored_sgs, multicolored_sor, ilu,
-                          ilut, multicolored_ilu
-                         };
+      ilut, multicolored_ilu};
 
   /**
-   * Solver used as smoother and to solve the coarse system.
+   * Solver used as smoother and to solve the coarse system of the multigrid.
    */
   enum mg_solver {richardson, cg, cr, bicgstab, gmres, fgmres};
 
@@ -375,79 +374,69 @@ namespace ParalutionWrappers
      */
     struct AdditionalData
     {
-      /**
-       * Constructor. This constructor lets Paralution choose the parameters
-       * of the algebraic multigrid.
-       */
-      AdditionalData (const unsigned int verbose = 0, const unsigned int max_iter = 2);
-
-      /**
-       * Constructor. The relaxation parameter is used only when smoothed
-       * aggregation is selected. Otherwise, this parameter is disregarded.
-       */
-      AdditionalData (const unsigned int verbose,
-                      const unsigned int max_iter,
-                      const mg_cycle cycle,
-                      const mg_interpolation interpolation,
-                      const mg_preconditioner preconditioner,
-                      const mg_solver smoother,
-                      const mg_solver coarse_solver,
-                      const unsigned int n_unknowns_coarse_level,
-                      const double relaxation_parameter = 1.2);
 
       /**
-       * The flag is true if the default parameters of Paralution should be
-       * used.
+       * Constructor. @p relaxation_parameter is used only with smoothed
+       * aggregation and is disregarded otherwise. @p over_interpolation is used
+       * only with aggregation and is disregared otherwise.
        */
-      bool default_parameters;
+      AdditionalData (const unsigned int verbose = 0,
+                      const mg_solver coarse_solver = gmres,
+                      const unsigned int n_unknowns_coarse_level = 300,
+                      const unsigned int max_iter = 2,
+                      const mg_cycle cycle = V_cycle,
+                      const mg_interpolation interpolation = smoothed_aggregation,
+                      const double coupling_strength = 0.01,
+                      const double relaxation_parameter = 2./3.,
+                      const double over_interpolation = 1.5);
 
       /**
-       * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
-       * cycle, or F cycle.
+       * Verbosity: 0 =  no output, 1 = print information at the
+       * start and at the end, 3 = print residual at each iteration.
        */
-      mg_cycle cycle;
+      unsigned int verbose;
 
       /**
-       * Interpolation used by the multigrid preconditioner: smoothed
-       * aggregation or aggregation.
+       * Solver of the coarse system.
        */
-      mg_interpolation interpolation;
+      mg_solver coarse_solver;
 
       /**
-       * Preconditioner for the smoother.
+       * Number of unknowns on the coarsest level.
        */
-      mg_preconditioner preconditioner;
+      unsigned int n_unknowns_coarse_level;
 
       /**
-       * Smoother.
+       * Maximum number of iterations.
        */
-      mg_solver smoother;
+      unsigned int max_iter;
 
       /**
-       * Solver of the coarse system.
+       * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
+       * cycle, or F cycle.
        */
-      mg_solver coarse_solver;
+      mg_cycle cycle;
 
       /**
-       * Verbosity: 0 =  no output, 1 = print information at the
-       * start and at the end, 3 = print residual at each iteration.
+       * Interpolation used by the multigrid preconditioner: smoothed
+       * aggregation or aggregation.
        */
-      unsigned int verbose;
+      mg_interpolation interpolation;
 
       /**
-       * Number of unknowns on the coarsest level.
+       * Coupling strength.
        */
-      unsigned int n_unknowns_coarse_level;
+      double coupling_strength;
 
       /**
-       * Maximum number of iterations.
+       * Relaxation parameter for smoothed interpolation aggregation.
        */
-      unsigned int max_iter;
+      double relaxation_parameter;
 
       /**
-       * Relaxation parameter for smoothed interpolation aggregation.
+       * Over-interpolation parameter for aggregation.
        */
-      double relaxation_parameter;
+      double over_interpolation;
     };
 
     /**
@@ -467,25 +456,10 @@ namespace ParalutionWrappers
 
   private :
     /**
-     * Free the smoothers, the preconditioners, and the coarse level solver.
+     * Free the coarse level solver.
      */
     void clear();
 
-    /**
-     * Number of levels of the multigrid.
-     */
-    unsigned int n_levels;
-
-    /**
-     * Smoothers for each level.
-     */
-    paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>, paralution::
-    LocalVector<Number>,Number> * *smoothers;
-    /**
-     * Preconditioners of the smoothers for each level.
-     */
-    paralution::Preconditioner<paralution::LocalMatrix<Number>, paralution::
-    LocalVector<Number>,Number> * *preconditioners;
     /**
      * Coarse grid solver.
      */
index 555d0481c037fb46f7a955a7861e9acfdd26659d..9b3a738cdb3f3de4b98c10ebb639c3a580d23797 100644 (file)
@@ -36,7 +36,8 @@ DEAL_II_NAMESPACE_OPEN
 namespace ParalutionWrappers
 {
   /**
-   * Base class for solver classes using the Paralution solvers.
+   * Base class for solver classes using the Paralution solvers except AMG which 
+   * is not derived from this class.
    *
    * @ingroup ParalutionWrappers
    * @author Bruno Turcksin 2013
@@ -68,10 +69,10 @@ namespace ParalutionWrappers
                        bool                            move_to_accelerator);
 
     /**
-     * Reference to the object that controls convergence of the iteratove
+     * Reference to the object that controls convergence of the iterative
      * solver. In fact, for these Paralution wrappers, Paralution does so
-     * itself, but we copy the data from this object solition process, and
-     * copy the data back into it afterwards.
+     * itself, but we copy the data from this object before starting solution 
+     * process, and copy the data back into it afterwards.
      */
     SolverControl &solver_control;
   };
@@ -134,7 +135,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -194,7 +195,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -249,7 +250,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -315,7 +316,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -371,7 +372,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -435,7 +436,7 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
@@ -499,7 +500,134 @@ namespace ParalutionWrappers
 
   private:
     /**
-     * store a copy of the flags for this solver.
+     * Store a copy of the flags for this solver.
+     */
+    const AdditionalData additional_data;
+  };
+
+
+
+  /**
+   * An implementation of the solver interface using the Paralution GMRES
+   * solver.
+   *
+   * @ingroup ParalutionWrappers
+   * @author Bruno Turcksin, 2014
+   */
+  class SolverAMG 
+  {
+  public:
+    /**
+     * Standardized data struct to pipe additional data to the solver.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor. @p relaxation_parameter is used only with smoothed
+       * aggregation and is disregarded otherwise. @p over_interpolation is used
+       * only with aggregation and is disregared otherwise.
+       */
+      AdditionalData (const unsigned int verbose = 0,
+                      const mg_solver coarse_solver = gmres,
+                      const unsigned int n_unknowns_coarse_level = 300,
+                      const mg_solver smoother = richardson,
+                      const mg_preconditioner = multicolored_sor,
+                      const mg_cycle cycle = V_cycle,
+                      const mg_interpolation interpolation = smoothed_aggregation,
+                      const double coupling_strength = 0.01,
+                      const double relaxation_parameter = 2./3.,
+                      const double over_interpolation = 1.5);
+
+      /**
+       * Verbosity: 0 =  no output, 1 = print information at the
+       * start and at the end, 3 = print residual at each iteration.
+       */
+      unsigned int verbose;
+
+      /**
+       * Solver of the coarse system.
+       */
+      mg_solver coarse_solver;
+
+      /**
+       * Number of unknowns on the coarsest level.
+       */
+      unsigned int n_unknowns_coarse_level;
+
+      /**
+       * Smoother (pre- and post-smoothers).
+       */
+      mg_solver smoother;
+
+      /**
+       * Preconditioner of the smoother.
+       */
+      mg_preconditioner preconditioner;
+
+      /**
+       * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
+       * cycle, or F cycle.
+       */
+      mg_cycle cycle;
+
+      /**
+       * Interpolation used by the multigrid preconditioner: smoothed
+       * aggregation or aggregation.
+       */
+      mg_interpolation interpolation;
+
+      /**
+       * Coupling strength.
+       */
+      double coupling_strength;
+
+      /**
+       * Relaxation parameter for smoothed interpolation aggregation.
+       */
+      double relaxation_parameter;
+
+      /**
+       * Over-interpolation parameter for aggregation.
+       */
+      double over_interpolation;
+    };
+
+    /**
+     * Constructor. AdditionalData is a structure that contains additional
+     * flags for tuning this particular solver.
+     */
+    SolverAMG (SolverControl        &cn,
+              const AdditionalData &data = AdditionalData());
+
+    /**
+     * Solve the linear system <tt>Ax=b</tt> using the AMG solver of
+     * Paralution. If the flag @p move_to_accelerator is set to true, the
+     * solver, the right-hand side, the matrix, and the solution vector are
+     * moved on the accelerator once the solver is built. The multigrid must be
+     * build on the CPU.
+     */
+    template <typename Number>
+    void solve (SparseMatrix<Number>     &A,
+                Vector<Number>           &x,
+                Vector<Number>           &b,
+                bool                      move_to_accelerator=false);
+
+    /**
+     * Access to object that controls convergence.
+     */
+    SolverControl &control() const;
+
+  private:
+    /**
+     * Reference to the object that controls convergence of the iterative
+     * solver. In fact, for these Paralution wrappers, Paralution does so
+     * itself, but we copy the data from this object before starting solution 
+     * process, and copy the data back into it afterwards.
+     */
+    SolverControl &solver_control;
+
+    /**
+     * Store a copy of the flags for this solver.
      */
     const AdditionalData additional_data;
   };
index c7b0b72926ea27bdc1d828d9abfbddcddab134ce..2d857a849c09dc9365523f9ff5694ffae8985baa 100644 (file)
@@ -135,7 +135,6 @@ namespace ParalutionWrappers
      * Change the dimension of the vector to @p N. The vector is filled with
      * zeros.
      */
-    //TODO look to add fast
     void reinit(const size_type N);
 
     /**
@@ -444,6 +443,18 @@ namespace ParalutionWrappers
      */
     void move_to_host_async();
 
+    /**
+     * Print to a stream. @p precision denotes the desired precision with
+     * which values shall be printed, @p scientific whether scientific
+     * notation shall be used. If @p across is @p true then the vector is
+     * printed in a line, while if @p false then the elements are printed on a
+     * separate line each.
+     */
+    void print (std::ostream      &out,         
+                const unsigned int precision,   
+                const bool         scientific,  
+                const bool         across) const;
+
     /**
      * Synchronize the code when move_to_host_async or move_to_accelerator_async
      * is used.
index 3729452c039a64c06d2c72eb83fb182fd932f986..2919a9c9c86c11d118a997bcab619fb187e0f395 100644 (file)
@@ -219,38 +219,26 @@ namespace ParalutionWrappers
 
   /* -------------------------- PreconditionAMG --------------------------- */
 
-  template <typename Number>
-  PreconditionAMG<Number>::AdditionalData::AdditionalData(const unsigned int verbose, 
-                                                          const unsigned int max_iter)
-    :
-    default_parameters(true),
-    verbose(verbose),
-    max_iter(max_iter)
-  {}
-
-
-
   template <typename Number>
   PreconditionAMG<Number>::AdditionalData::AdditionalData(const unsigned int verbose,
-                                                          const unsigned int max_iter,
-                                                          const mg_cycle cycle,
-                                                          const mg_interpolation interpolation,
-                                                          const mg_preconditioner preconditioner,
-                                                          const mg_solver smoother,
                                                           const mg_solver coarse_solver,
                                                           const unsigned int n_unknowns_coarse_level,
-                                                          const double relaxation_parameter)
+                                                          const unsigned int max_iter,
+                                                          const mg_cycle cycle,                           
+                                                          const mg_interpolation interpolation,
+                                                          const double coupling_strength,
+                                                          const double relaxation_parameter,
+                                                          const double over_interpolation)                     
     :
-    default_parameters(false),
-    cycle(cycle),
-    interpolation(interpolation),
-    preconditioner(preconditioner),
-    smoother(smoother),
-    coarse_solver(coarse_solver),
     verbose(verbose),
+    coarse_solver(coarse_solver),
     n_unknowns_coarse_level(n_unknowns_coarse_level),
     max_iter(max_iter),
-    relaxation_parameter(relaxation_parameter)
+    cycle(cycle),
+    interpolation(interpolation),
+    coupling_strength(coupling_strength),
+    relaxation_parameter(relaxation_parameter),
+    over_interpolation(over_interpolation)
   {}
 
 
@@ -258,9 +246,6 @@ namespace ParalutionWrappers
   template <typename Number>
   PreconditionAMG<Number>::PreconditionAMG(const AdditionalData &additional_data)
     :
-    n_levels(1),
-    smoothers(NULL),
-    preconditioners(NULL),
     coarse_solver(NULL)
   {
     this->preconditioner.reset(new paralution::AMG<paralution::LocalMatrix<Number>,
@@ -284,318 +269,137 @@ namespace ParalutionWrappers
   {
     // Downcast the preconditioner pointer
     paralution::AMG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
-               Number>* downcasted_ptr = static_cast<paralution::AMG<paralution::
-                                         LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
+      Number>* downcasted_ptr = static_cast<paralution::AMG<paralution::
+        LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
 
     // Set the maximum number of iterations.
     downcasted_ptr->InitMaxIter(additional_data.max_iter);
-    
-    // Set the verbosity
+
+    // Set the verbosity of the multigrid.
     downcasted_ptr->Verbose(additional_data.verbose);
 
-    if (additional_data.default_parameters==false)
-      {
-        // Free the smoothers, the preconditioners, and the coarse level solver,
-        // if necessary.
-        clear();
-
-        // Set the number of unknowns on the coarsest level.
-        downcasted_ptr->SetCoarsestLevel(additional_data.n_unknowns_coarse_level);
-
-        // Set the interpolation type.
-        switch (additional_data.interpolation)
-          {
-          case aggregation :
-          {
-            downcasted_ptr->SetInterpolation(paralution::Aggregation);
-            break;
-          }
-          case smoothed_aggregation :
-          {
-            downcasted_ptr->SetInterpolation(paralution::SmoothedAggregation);
-            downcasted_ptr->SetInterpRelax(additional_data.relaxation_parameter);
-            break;
-          }
-          default :
-          {
-            AssertThrow(false,ExcMessage("Unknown interpolation type for PreconditionAMG."));
-          }
-          }
-
-        // Set the type of cycle
-        switch (additional_data.cycle)
-          {
-          case V_cycle :
-          {
-            downcasted_ptr->SetCycle(paralution::Vcycle);
-            break;
-          }
-          case  W_cycle :
-          {
-            downcasted_ptr->SetCycle(paralution::Wcycle);
-            break;
-          }
-          case K_cycle :
-          {
-            downcasted_ptr->SetCycle(paralution::Kcycle);
-            break;
-          }
-          case F_cycle :
-          {
-            downcasted_ptr->SetCycle(paralution::Fcycle);
-          }
-          default :
-          {
-            AssertThrow(false,ExcMessage("Unknown cycle type for PreconditionAMG."));
-          }
-          }
-
-    // Use manual smoothers.
-    downcasted_ptr->SetManualSmoothers(true);
+    // Free the coarse level solver if necessary.
+    clear();
+
+    // Set the number of unknowns on the coarsest level.
+    downcasted_ptr->SetCoarsestLevel(additional_data.n_unknowns_coarse_level);
+
+    // Set the interpolation type.
+    switch (additional_data.interpolation)
+    {
+      case aggregation :
+        {
+          downcasted_ptr->SetInterpolation(paralution::Aggregation);
+          downcasted_ptr->SetOverInterp(additional_data.over_interpolation);
+          break;
+        }
+      case smoothed_aggregation :
+        {
+          downcasted_ptr->SetInterpolation(paralution::SmoothedAggregation);
+          downcasted_ptr->SetInterpRelax(additional_data.relaxation_parameter);
+          break;
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown interpolation type for PreconditionAMG."));
+        }
+    }
+
+    // Set the coupling strength
+    downcasted_ptr->SetCouplingStrength(additional_data.coupling_strength);
+
+    // Set the type of cycle
+    switch (additional_data.cycle)
+    {
+      case V_cycle :
+        {
+          downcasted_ptr->SetCycle(paralution::Vcycle);
+          break;
+        }
+      case  W_cycle :
+        {
+          downcasted_ptr->SetCycle(paralution::Wcycle);
+          break;
+        }
+      case K_cycle :
+        {
+          downcasted_ptr->SetCycle(paralution::Kcycle);
+          break;
+        }
+      case F_cycle :
+        {
+          downcasted_ptr->SetCycle(paralution::Fcycle);
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown cycle type for PreconditionAMG."));
+        }
+    }
 
     // Use manual coarse grid solver.
     downcasted_ptr->SetManualSolver(true);
 
-    // Build the hierarchy of the grids.
-    downcasted_ptr->BuildHierarchy();
-
-    // Get the number of grids/levels.
-    n_levels = downcasted_ptr->GetNumLevels();
 
-    // Smoothers for each level.
-    smoothers = new paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
-    paralution::LocalVector<Number>, Number>* [n_levels-1];
-    switch (additional_data.smoother)
-      {
+    // Coarse solver
+    switch (additional_data.coarse_solver)
+    {
       case richardson :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
+        {
+          paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       case cg :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::CG<paralution::LocalMatrix<Number>,
+        {
+          paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::CG<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       case cr :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::CR<paralution::LocalMatrix<Number>,
+        {
+          paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::CR<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       case bicgstab :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
+        {
+          paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs =  new paralution::BiCGStab<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       case gmres :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::GMRES<paralution::LocalMatrix<Number>,
+        {
+          paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::GMRES<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       case fgmres :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *sm = new paralution::FGMRES<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            smoothers[i] = sm;
-          }
-        break;
-      }
-      default :
-      {
-        AssertThrow(false,ExcMessage("Unknown smoother for PreconditionAMG."));
-      }
-      }
-
-    // Preconditioners for each level.
-    preconditioners = new paralution::Preconditioner<paralution::LocalMatrix<Number>,
-    paralution::LocalVector<Number>, Number>*[n_levels-1];
-    switch (additional_data.preconditioner)
-      {
-      case jacobi :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::Jacobi<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *precond = new paralution::Jacobi<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case sgs :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::SGS<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *precond = new paralution::SGS<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case multicolored_sgs :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
-                       paralution::LocalVector<Number>, Number> *precond =
-                         new paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case multicolored_sor :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
-                       paralution::LocalVector<Number>, Number> *precond =
-                         new paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case ilu :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::ILU<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *precond = new paralution::ILU<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case ilut :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::ILUT<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                       Number> *precond = new paralution::ILUT<paralution::LocalMatrix<Number>,
-            paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
-      case multicolored_ilu :
-      {
-        for (unsigned int i=0; i<n_levels-1; ++i)
-          {
-            paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
-                       paralution::LocalVector<Number>, Number> *precond =
-                         new paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
+        {
+          paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::FGMRES<paralution::LocalMatrix<Number>,
             paralution::LocalVector<Number>, Number>;
-            preconditioners[i] = precond;
-          }
-        break;
-      }
+          coarse_solver = cs;
+          break;
+        }
       default :
-      {
-        AssertThrow(false,ExcMessage("Unknown preconditioner for the smoother for PreconditionAMG."));
-      }
-      }
-
-    for (unsigned int i=0; i<n_levels-1; ++i)
-      {
-        smoothers[i]->SetPreconditioner(*preconditioners[i]);
-        smoothers[i]->Verbose(0);
-      }
-    downcasted_ptr->SetSmoother(smoothers);
-
-    // Coarse solver
-    switch (additional_data.coarse_solver)
-      {
-      case richardson :
-      {
-        paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      case cg :
-      {
-        paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs = new paralution::CG<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      case cr :
-      {
-        paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs = new paralution::CR<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      case bicgstab :
-      {
-        paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs =  new paralution::BiCGStab<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      case gmres :
-      {
-        paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs = new paralution::GMRES<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      case fgmres :
-      {
-        paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
-                   Number> *cs = new paralution::FGMRES<paralution::LocalMatrix<Number>,
-        paralution::LocalVector<Number>, Number>;
-        coarse_solver = cs;
-        break;
-      }
-      default :
-      {
-        AssertThrow(false,ExcMessage("Unknown smoother for PreconditionAMG."));
-      }
-      }
+        {
+          AssertThrow(false,ExcMessage("Unknown coarse solver for PreconditionAMG."));
+        }
+    }
 
+    // Set the coarse solver and the verbosity of the coarse solver.
+    coarse_solver->Verbose(additional_data.verbose);
     downcasted_ptr->SetSolver(*coarse_solver);
-      }
   }
 
 
@@ -609,36 +413,6 @@ namespace ParalutionWrappers
         delete coarse_solver;
         coarse_solver = NULL;
       }
-
-    // Free the preconditioners
-    for (unsigned int i=0; i<n_levels-1; ++i)
-      {
-        if (preconditioners[i]!=NULL)
-          {
-            delete preconditioners[i];
-            preconditioners[i] = NULL;
-          }
-      }
-    if (preconditioners!=NULL)
-      {
-        delete [] preconditioners;
-        preconditioners = NULL;
-      }
-
-    // Free the smoothers.
-    for (unsigned int i=0; i<n_levels-1; ++i)
-      {
-        if (smoothers[i]!=NULL)
-          {
-            delete smoothers[i];
-            smoothers[i] = NULL;
-          }
-      }
-    if (smoothers!=NULL)
-      {
-        delete [] smoothers;
-        smoothers = NULL;
-      }
   }
 }
 
index 20f15b08c9a00646ae1079b6fc6ad81170ef875b..aadf732b87e747d34bfb3381c679cae1b5b7a0fc 100644 (file)
@@ -29,7 +29,7 @@ namespace ParalutionWrappers
 
 
 
-  SolverControl &SolverBase::control() const
+  SolverControlSolverBase::control() const
   {
     return solver_control;
   }
@@ -357,6 +357,468 @@ namespace ParalutionWrappers
 
     this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
   }
+
+
+
+  /* ---------------------- SolverAMG------------------------- */
+
+  SolverAMG::AdditionalData::AdditionalData(const unsigned int verbose,
+                                            const mg_solver coarse_solver,
+                                            const unsigned int n_unknowns_coarse_level,
+                                            const mg_solver smoother,
+                                            const mg_preconditioner preconditioner,
+                                            const mg_cycle cycle,
+                                            const mg_interpolation interpolation,
+                                            const double coupling_strength,
+                                            const double relaxation_parameter,
+                                            const double over_interpolation)
+    :
+    verbose(verbose),
+    coarse_solver(coarse_solver),
+    n_unknowns_coarse_level(n_unknowns_coarse_level),
+    smoother(smoother),
+    preconditioner(preconditioner),
+    cycle(cycle),
+    interpolation(interpolation),
+    coupling_strength(coupling_strength),
+    relaxation_parameter(relaxation_parameter),
+    over_interpolation(over_interpolation)
+  {}
+
+
+
+  SolverAMG::SolverAMG (SolverControl        &cn,
+                        const AdditionalData &data)
+    :
+    solver_control (cn),
+    additional_data (data)
+  {}
+
+
+
+  template <typename Number>
+  void SolverAMG::solve(SparseMatrix<Number>           &A,
+                        Vector<Number>                 &x,
+                        Vector<Number>                 &b,
+                        bool                            move_to_accelerator)
+  {
+    paralution::AMG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,Number> solver;
+
+    // Set the system to solve.
+    solver.SetOperator(A.paralution_matrix());
+    solver.SetOperatorFormat(CSR);
+    
+    // Set the verbosity of the multigrid.
+    solver.Verbose(additional_data.verbose);
+    
+    // Set number of unknowns on coarsest level.
+    solver.SetCoarsestLevel(additional_data.n_unknowns_coarse_level);
+
+    // Set the interpolation type.
+    switch (additional_data.interpolation)
+    {
+      case aggregation :
+        {
+          solver.SetInterpolation(paralution::Aggregation);
+          solver.SetOverInterp(additional_data.over_interpolation);
+          break;
+        }
+      case smoothed_aggregation :
+        {
+          solver.SetInterpolation(paralution::SmoothedAggregation);
+          solver.SetInterpRelax(additional_data.relaxation_parameter);
+          break;
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown interpolation type for PreconditionAMG."));
+        }
+    }
+
+    // Set the coupling strength.
+    solver.SetCouplingStrength(additional_data.coupling_strength);
+
+    // Set the type of cycle
+    switch (additional_data.cycle)
+    {
+      case V_cycle :
+        {
+          solver.SetCycle(paralution::Vcycle);
+          break;
+        }
+      case  W_cycle :
+        {
+          solver.SetCycle(paralution::Wcycle);
+          break;
+        }
+      case K_cycle :
+        {
+          solver.SetCycle(paralution::Kcycle);
+          break;
+        }
+      case F_cycle :
+        {
+          solver.SetCycle(paralution::Fcycle);
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown cycle type for PreconditionAMG."));
+        }
+    }
+    
+    // Manual smoothers
+    solver.SetManualSmoothers(true);
+    
+    // Manual course grid solver
+    solver.SetManualSolver(true);
+    
+    // Set grid transfer scaling
+    solver.SetScaling(true);
+
+    // Build the grid hierarchy.
+    solver.BuildHierarchy();
+
+    unsigned int n_levels = solver.GetNumLevels();
+
+    // Smoother for each level
+    paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+      paralution::LocalVector<Number>,Number> **smoothers = NULL;
+    smoothers = new paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+              paralution::LocalVector<Number>,Number>* [n_levels-1];
+
+    // Preconditioner for each smoother
+    paralution::Preconditioner<paralution::LocalMatrix<Number>,
+      paralution::LocalVector<Number>,Number> **preconds = NULL;
+    preconds = new paralution::Preconditioner<paralution::LocalMatrix<Number>,
+             paralution::LocalVector<Number>,Number>* [n_levels-1];
+
+    // Coarse Grid Solver
+    paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+      paralution::LocalVector<Number>,Number> *coarse_solver;
+
+    // Smoother
+    switch (additional_data.smoother)
+    {
+      case richardson :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::FixedPoint<paralution::LocalMatrix<Number>,
+              paralution::LocalVector<Number>,Number> *richardson = NULL;
+            richardson = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
+                       paralution::LocalVector<Number>,Number>;
+            smoothers[i] = richardson;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      case cg :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::CG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *cg = NULL;
+            cg = new paralution::CG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+               Number>;
+            smoothers[i] = cg;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      case cr :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::CR<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *cr = NULL;
+            cr = new paralution::CR<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+               Number>;
+            smoothers[i] = cr;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      case bicgstab :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::BiCGStab<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *bicgstab = NULL;
+            bicgstab = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
+                     paralution::LocalVector<Number>,Number>;
+            smoothers[i] = bicgstab;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      case gmres :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::GMRES<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *gmres = NULL;
+            gmres = new paralution::GMRES<paralution::LocalMatrix<Number>,
+                  paralution::LocalVector<Number>,Number>;
+            smoothers[i] = gmres;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      case fgmres :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::FGMRES<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *fgmres = NULL;
+            fgmres = new paralution::FGMRES<paralution::LocalMatrix<Number>,
+                  paralution::LocalVector<Number>,Number>;
+            smoothers[i] = fgmres;
+            smoothers[i]->Verbose(additional_data.verbose);
+          }
+
+          break;
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown smoother for AMG."));
+        }
+    }
+
+    // Preconditioner
+    switch (additional_data.preconditioner)
+    {
+      case jacobi :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::Jacobi<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *jacobi = NULL;
+            jacobi = new paralution::Jacobi<paralution::LocalMatrix<Number>,
+                   paralution::LocalVector<Number>,Number>;
+            preconds[i] = jacobi;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case sgs :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::SGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *sgs = NULL;
+            sgs = new paralution::SGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+                Number>;
+            preconds[i] = sgs;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case multicolored_sgs :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *multicolored_sgs = NULL;
+            multicolored_sgs = new paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
+                             paralution::LocalVector<Number>,Number>;
+            multicolored_sgs->SetPrecondMatrixFormat(ELL);
+            preconds[i] = multicolored_sgs;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case multicolored_sor :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::MultiColoredGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *multicolored_sor = NULL;
+            multicolored_sor = new paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
+                     paralution::LocalVector<Number>,Number>;
+            multicolored_sor->SetPrecondMatrixFormat(ELL);
+            preconds[i] = multicolored_sor;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case ilu :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::ILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *ilu = NULL;
+            ilu = new paralution::ILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+                Number>;
+            preconds[i] = ilu;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case ilut :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::ILUT<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *ilut = NULL;
+            ilut = new paralution::ILUT<paralution::LocalMatrix<Number>,
+                  paralution::LocalVector<Number>,Number>;
+            preconds[i] = ilut;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      case multicolored_ilu :
+        {
+          for (unsigned int i=0; i<n_levels-1; ++i)
+          {
+            paralution::MultiColoredILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+              Number> *multicolored_ilu = NULL;
+            multicolored_ilu = new paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
+                     paralution::LocalVector<Number>,Number>;
+            multicolored_ilu->SetPrecondMatrixFormat(ELL);
+            preconds[i] = multicolored_ilu;
+            smoothers[i]->SetPreconditioner(*preconds[i]);
+          }
+
+          break;
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown preconditioner for the smoother for AMG."));
+        }
+    }
+
+    // Set the smoothers.
+    solver.SetSmoother(smoothers);
+    solver.SetSmootherPreIter(1);
+    solver.SetSmootherPostIter(2);
+
+    // Coarse solver
+    switch (additional_data.coarse_solver)
+    {
+      case richardson :
+        {
+          paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      case cg :
+        {
+          paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::CG<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      case cr :
+        {
+          paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::CR<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      case bicgstab :
+        {
+          paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs =  new paralution::BiCGStab<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      case gmres :
+        {
+          paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::GMRES<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      case fgmres :
+        {
+          paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+            Number> *cs = new paralution::FGMRES<paralution::LocalMatrix<Number>,
+            paralution::LocalVector<Number>, Number>;
+          coarse_solver = cs;
+          break;
+        }
+      default :
+        {
+          AssertThrow(false,ExcMessage("Unknown coarse solver for PreconditionAMG."));
+        }
+    }
+    solver.SetSolver(*coarse_solver);
+
+    // Set absolute tolerance, relative tolerance, divergence tolerance,
+    // maximum number of iterations.
+    solver.Init(solver_control.tolerance(),0.,1.e100,solver_control.max_steps());
+
+    // Build the solver.
+    solver.Build();
+
+    // Move the solver to the accelerator if necessary.
+    if (move_to_accelerator==true)
+      {
+        A.move_to_accelerator();
+        x.move_to_accelerator();
+        b.move_to_accelerator();
+        solver.MoveToAccelerator();
+      }
+    solver.Solve(b.paralution_vector(),&(x.paralution_vector()));
+
+    // Let the deal.II SolverControl object know what has happened. If the solve
+    // succeeded, the status of the solver control will turn into
+    // SolverControl::success. solver_control only exists on the host so the solver
+    // needs to be on the most. 
+    solver.MoveToHost();
+
+    solver_control.check (solver.GetIterationCount(), solver.GetCurrentResidual());
+
+    if (solver_control.last_check() != SolverControl::success)
+      AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(),
+                                                       solver_control.last_value()));
+
+    // Free the coarse solver, the preconditioners and the smoothers.
+    if (coarse_solver!=NULL)
+    {
+      delete coarse_solver;
+      coarse_solver = NULL;
+    }
+    if (preconds!=NULL)
+    {
+      for (unsigned int i=0; i<n_levels-1; ++i)
+        delete [] preconds[i];
+      delete [] preconds;
+      preconds = NULL;
+    }
+    if (smoothers!=NULL)
+    {
+      for (unsigned int i=0; i<n_levels-1; ++i)
+        delete [] smoothers[i];
+      delete [] smoothers;
+      smoothers = NULL;
+    }
+  }
+
+
+
+  SolverControl& SolverAMG::control() const
+  {
+    return solver_control;
+  }
 }
 
 
@@ -435,6 +897,16 @@ namespace ParalutionWrappers
                                             Vector<double>           &b,
                                             PreconditionBase<double> &preconditioner,
                                             bool                      move_to_accelerator);
+
+  template void SolverAMG::solve<float>(SparseMatrix<float> &A,
+                                        Vector<float>       &x,
+                                        Vector<float>       &b,
+                                        bool                 move_to_accelerator);
+
+  template void SolverAMG::solve<double>(SparseMatrix<double> &A,
+                                         Vector<double>       &x,
+                                         Vector<double>       &b,
+                                         bool                  move_to_accelerator);
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 87b7217ab01de381ec3065cc89e3f58379975c5d..bd14accea92f6651d98d7c0004fe386b0ad17dc8 100644 (file)
@@ -240,6 +240,36 @@ namespace ParalutionWrappers
     local_vector.ScaleAdd2(0.,u.paralution_vector(),a,v.paralution_vector(),b);
     local_vector.AddScale(w.paralution_vector(),c);
   }
+
+
+
+  template <typename Number>
+  void Vector<Number>::print (std::ostream      &out,
+                              const unsigned int precision,
+                              const bool         scientific,
+                              const bool         across) const
+  {
+    AssertThrow (out, ExcIO());
+
+    // get a representation of the vector and loop over all the elements 
+    out.precision (precision);
+    if (scientific)
+      out.setf (std::ios::scientific, std::ios::floatfield);
+    else
+      out.setf (std::ios::fixed, std::ios::floatfield);
+
+    if (across)
+      for (size_type i=0; i<size(); ++i)
+        out << static_cast<double>(local_vector[i]) << ' ';
+    else
+      for (size_type i=0; i<size(); ++i)
+        out << static_cast<double>(local_vector[i]) << std::endl;
+    out << std::endl;
+
+    // restore the representation
+    // of the vector
+    AssertThrow (out, ExcIO());
+  }
 }
 
 // Explicit instantiations

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.