TrilinosWrappers::SolverBase::Solve support Epetra_Operator for matrix free calculation.

git-svn-id: https://svn.dealii.org/trunk@29811 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/lac/petsc_matrix_free.h b/deal.II/include/deal.II/lac/petsc_matrix_free.h
index c7217178d3066b5c0bc0a107e259d5b0516fdd53..f5580a0cce0860f68e87c3e0347083ee5857e644 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_matrix_free.h
+++ b/deal.II/include/deal.II/lac/petsc_matrix_free.h
@@ -34,7 +34,7 @@ namespace PETScWrappers
    * Implementation of a parallel matrix class based on PETSc <tt>MatShell</tt> matrix-type.
    * This base class implements only the interface to the PETSc matrix object,
    * while all the functionality is contained in the matrix-vector
-   * multiplication which must be reimplmented in derived classes.
+   * multiplication which must be reimplemented in derived classes.
    *
    * This interface is an addition to the dealii::MatrixFree class to realize
    * user-defined matrix-classes together with PETSc solvers and functionalities.

diff --git a/deal.II/include/deal.II/lac/trilinos_solver.h b/deal.II/include/deal.II/lac/trilinos_solver.h
index a0ed04362361110d832982b8f1bdd66b46d7bfa1..54a72beaf1a518a4b6bfa5be061ff0f38c4372d8 100644 (file)
--- a/deal.II/include/deal.II/lac/trilinos_solver.h
+++ b/deal.II/include/deal.II/lac/trilinos_solver.h
@@ -155,6 +155,26 @@ namespace TrilinosWrappers
            const VectorBase       &b,
            const PreconditionBase &preconditioner);
 
+    /**
+     * Solve the linear system
+     * <tt>Ax=b</tt> where <tt>A</tt> 
+     * is an operator. This function 
+     * can be used for matrix free 
+     * computation. Depending on
+     * the information provided by
+     * derived classes and the
+     * object passed as a
+     * preconditioner, one of the
+     * linear solvers and
+     * preconditioners of Trilinos
+     * is chosen.
+     */
+    void
+    solve (Epetra_Operator        &A,
+           VectorBase             &x,
+           const VectorBase       &b,
+           const PreconditionBase &preconditioner);
+
     /**
      * Solve the linear system
      * <tt>Ax=b</tt>. Depending on the
@@ -179,6 +199,32 @@ namespace TrilinosWrappers
            const dealii::Vector<double> &b,
            const PreconditionBase       &preconditioner);
 
+    /**
+     * Solve the linear system
+     * <tt>Ax=b</tt> where <tt>A</tt>
+     * is an operator. This function can 
+     * be used for matric free. Depending on the
+     * information provided by derived
+     * classes and the object passed as a
+     * preconditioner, one of the linear
+     * solvers and preconditioners of
+     * Trilinos is chosen. This class
+     * works with matrices according to
+     * the TrilinosWrappers format, but
+     * can take deal.II vectors as
+     * argument. Since deal.II are serial
+     * vectors (not distributed), this
+     * function does only what you expect
+     * in case the matrix is locally
+     * owned. Otherwise, an exception
+     * will be thrown.
+     */
+    void
+    solve (Epetra_Operator              &A,
+           dealii::Vector<double>       &x,
+           const dealii::Vector<double> &b,
+           const PreconditionBase       &preconditioner);
+
     /**
      * Access to object that controls
      * convergence.
@@ -211,6 +257,12 @@ namespace TrilinosWrappers
 
   private:
 
+    /**
+     * The solve function is used to set properly the Epetra_LinearProblem,
+     * once it is done this function solves the linear problem.
+     */
+    void execute_solve(const PreconditionBase &preconditioner);
+
     /**
      * A structure that collects
      * the Trilinos sparse matrix,

diff --git a/deal.II/source/lac/trilinos_solver.cc b/deal.II/source/lac/trilinos_solver.cc
index 60ae66c5de613c98378c27bfa10f9e96c3183073..c315f191397f81ccabc3a5abdd7888f101a93d1c 100644 (file)
--- a/deal.II/source/lac/trilinos_solver.cc
+++ b/deal.II/source/lac/trilinos_solver.cc
@@ -73,8 +73,6 @@ namespace TrilinosWrappers
                      const VectorBase       &b,
                      const PreconditionBase &preconditioner)
   {
-    int ierr;
-
     linear_problem.reset();
 
     // We need an
@@ -87,86 +85,30 @@ namespace TrilinosWrappers
                               &x.trilinos_vector(),
                               const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
 
-    // Next we can allocate the
-    // AztecOO solver...
-    solver.SetProblem(*linear_problem);
+    execute_solve(preconditioner);
+  }
 
-    // ... and we can specify the
-    // solver to be used.
-    switch (solver_name)
-      {
-      case cg:
-        solver.SetAztecOption(AZ_solver, AZ_cg);
-        break;
-      case cgs:
-        solver.SetAztecOption(AZ_solver, AZ_cgs);
-        break;
-      case gmres:
-        solver.SetAztecOption(AZ_solver, AZ_gmres);
-        solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter);
-        break;
-      case bicgstab:
-        solver.SetAztecOption(AZ_solver, AZ_bicgstab);
-        break;
-      case tfqmr:
-        solver.SetAztecOption(AZ_solver, AZ_tfqmr);
-        break;
-      default:
-        Assert (false, ExcNotImplemented());
-      }
 
-    // Introduce the
-    // preconditioner, ...
-    ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
-                                   (preconditioner.preconditioner.get()));
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    // ... set some options, ...
-    solver.SetAztecOption (AZ_output, additional_data.output_solver_details ?
-                           AZ_all : AZ_none);
-    solver.SetAztecOption (AZ_conv, AZ_noscaled);
-
-    // ... and then solve!
-    ierr = solver.Iterate (solver_control.max_steps(),
-                           solver_control.tolerance());
-
-    // report errors in more detail
-    // than just by checking whether
-    // the return status is zero or
-    // greater. the error strings are
-    // taken from the implementation
-    // of the AztecOO::Iterate
-    // function
-    switch (ierr)
-      {
-      case -1:
-        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -1: "
-                                       "option not implemented"));
-      case -2:
-        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -2: "
-                                       "numerical breakdown"));
-      case -3:
-        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -3: "
-                                       "loss of precision"));
-      case -4:
-        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: "
-                                       "GMRES hessenberg ill-conditioned"));
-      default:
-        AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
-      }
+  void
+  SolverBase::solve (Epetra_Operator        &A,
+                     VectorBase             &x,
+                     const VectorBase       &b,
+                     const PreconditionBase &preconditioner)
+  {
+    linear_problem.reset();
 
-    // Finally, 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.check (solver.NumIters(), solver.TrueResidual());
+    // We need an
+    // Epetra_LinearProblem object
+    // to let the AztecOO solver
+    // know about the matrix and
+    // vectors.
+    linear_problem.reset
+    (new Epetra_LinearProblem(&A,
+                              &x.trilinos_vector(),
+                              const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
 
-    if (solver_control.last_check() != SolverControl::success)
-      throw SolverControl::NoConvergence (solver_control.last_step(),
-                                          solver_control.last_value());
+    execute_solve(preconditioner);
   }
 
 
@@ -177,8 +119,6 @@ namespace TrilinosWrappers
                      const dealii::Vector<double> &b,
                      const PreconditionBase       &preconditioner)
   {
-    int ierr;
-
     linear_problem.reset();
 
     // In case we call the solver with
@@ -205,6 +145,39 @@ namespace TrilinosWrappers
                           (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
                            &ep_x, &ep_b));
 
+    execute_solve(preconditioner);
+  }
+
+
+
+  void
+  SolverBase::solve (Epetra_Operator              &A,
+                     dealii::Vector<double>       &x,
+                     const dealii::Vector<double> &b,
+                     const PreconditionBase       &preconditioner)
+  {
+    linear_problem.reset();
+
+    Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.begin());
+    Epetra_Vector ep_b (View, A.OperatorRangeMap(), const_cast<double *>(b.begin()));
+
+    // We need an
+    // Epetra_LinearProblem object
+    // to let the AztecOO solver
+    // know about the matrix and
+    // vectors.
+    linear_problem.reset (new Epetra_LinearProblem(&A,&ep_x, &ep_b));
+
+    execute_solve(preconditioner);
+  }
+
+
+
+  void
+  SolverBase::execute_solve(const PreconditionBase &preconditioner)
+  {
+    int ierr;
+
     // Next we can allocate the
     // AztecOO solver...
     solver.SetProblem(*linear_problem);
@@ -247,7 +220,31 @@ namespace TrilinosWrappers
     // ... and then solve!
     ierr = solver.Iterate (solver_control.max_steps(),
                            solver_control.tolerance());
-    AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+
+    // report errors in more detail
+    // than just by checking whether
+    // the return status is zero or
+    // greater. the error strings are
+    // taken from the implementation
+    // of the AztecOO::Iterate
+    // function
+    switch (ierr)
+      {
+      case -1:
+        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -1: "
+                                       "option not implemented"));
+      case -2:
+        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -2: "
+                                       "numerical breakdown"));
+      case -3:
+        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -3: "
+                                       "loss of precision"));
+      case -4:
+        AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: "
+                                       "GMRES hessenberg ill-conditioned"));
+      default:
+        AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+      }
 
     // Finally, let the deal.II
     // SolverControl object know