]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Patch by Uwe Koecher: select among Trilinos solvers
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 9 Jul 2014 08:01:15 +0000 (08:01 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 9 Jul 2014 08:01:15 +0000 (08:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@33117 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/trilinos_precondition.h
deal.II/include/deal.II/lac/trilinos_solver.h
deal.II/source/lac/trilinos_solver.cc
tests/trilinos/direct_solver.cc [new file with mode: 0644]
tests/trilinos/direct_solver.output [new file with mode: 0644]

index 974eef62d323188401c808bec1949847fa722075..635384ed356d8c3f027bc12f1f1ea14a371ae90a 100644 (file)
@@ -148,6 +148,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Improved: The Trilinos direct solver, TrilinosWrappers::SolverDirect,
+  now takes a string to select among the available Amesos solvers. Moreover,
+  the solver now also supports deal.II's distributed vectors.
+  <br>
+  (Uwe K&ouml;cher, Martin Kronbichler, 2014/07/09)
+  </li>
+
   <li> New: There are now three new preconditioner classes
   TrilinosWrappers::PreconditionBlockJacobi,
   TrilinosWrappers::PreconditionBlockSSOR, and
index 01806b9646e1d37faa75beb2f592e6244ede2b18..c1ceb62471c8dd6dca6efc879d11055827adec5a 100644 (file)
@@ -1416,34 +1416,36 @@ namespace TrilinosWrappers
       /**
        * Determines which smoother to use for the AMG cycle. Possibilities
        * for smoother_type are the following:
-       *     "Aztec"
-       *     "IFPACK"
-       *     "Jacobi"
-       *     "ML symmetric Gauss-Seidel"
-       *     "symmetric Gauss-Seidel"
-       *     "ML Gauss-Seidel"
-       *     "Gauss-Seidel"
-       *     "block Gauss-Seidel"
-       *     "symmetric block Gauss-Seidel"
-       *     "Chebyshev"
-       *     "MLS"
-       *     "Hiptmair"
-       *     "Amesos-KLU"
-       *     "Amesos-Superlu"
-       *     "Amesos-UMFPACK"
-       *     "Amesos-Superludist"
-       *     "Amesos-MUMPS"
-       *     "user-defined"
-       *     "SuperLU"
-       *     "IFPACK-Chebyshev"
-       *     "self"
-       *     "do-nothing"
-       *     "IC"
-       *     "ICT"
-       *     "ILU"
-       *     "ILUT"
-       *     "Block Chebyshev"
-       *     "IFPACK-Block Chebyshev"
+       * <ol>
+       *   <li>  "Aztec" </li>
+       *   <li>  "IFPACK" </li>
+       *   <li>  "Jacobi" </li>
+       *   <li>  "ML symmetric Gauss-Seidel" </li>
+       *   <li>  "symmetric Gauss-Seidel" </li>
+       *   <li>  "ML Gauss-Seidel" </li>
+       *   <li>  "Gauss-Seidel" </li>
+       *   <li>  "block Gauss-Seidel" </li>
+       *   <li>  "symmetric block Gauss-Seidel" </li>
+       *   <li>  "Chebyshev" </li>
+       *   <li>  "MLS" </li>
+       *   <li>  "Hiptmair" </li>
+       *   <li>  "Amesos-KLU" </li>
+       *   <li>  "Amesos-Superlu" </li>
+       *   <li>  "Amesos-UMFPACK" </li>
+       *   <li>  "Amesos-Superludist" </li>
+       *   <li>  "Amesos-MUMPS" </li>
+       *   <li>  "user-defined" </li>
+       *   <li>  "SuperLU" </li>
+       *   <li>  "IFPACK-Chebyshev" </li>
+       *   <li>  "self" </li>
+       *   <li>  "do-nothing" </li>
+       *   <li>  "IC" </li>
+       *   <li>  "ICT" </li>
+       *   <li>  "ILU" </li>
+       *   <li>  "ILUT" </li>
+       *   <li>  "Block Chebyshev" </li>
+       *   <li>  "IFPACK-Block Chebyshev" </li>
+       * </ol>
        */
       const char* smoother_type;
 
index 1c4ead78a83c8f93fe3ac26c81ba5791058711dd..9ccc910391269df3052ce8668fc446b390a72ce1 100644 (file)
@@ -26,6 +26,7 @@
 #  include <deal.II/lac/exceptions.h>
 #  include <deal.II/lac/solver_control.h>
 #  include <deal.II/lac/vector.h>
+#  include <deal.II/lac/parallel_vector.h>
 
 #  include <Epetra_LinearProblem.h>
 #  include <AztecOO.h>
@@ -67,18 +68,11 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Enumeration object that is
-     * set in the constructor of
-     * the derived classes and
-     * tells Trilinos which solver
-     * to use. This option can also
-     * be set in the user program,
-     * so one might use this base
-     * class instead of one of the
-     * specialized derived classes
-     * when the solver should be
-     * set at runtime. Currently
-     * enabled options are:
+     * Enumeration object that is set in the constructor of the derived
+     * classes and tells Trilinos which solver to use. This option can also be
+     * set in the user program, so one might use this base class instead of
+     * one of the specialized derived classes when the solver should be set at
+     * runtime. Currently enabled options are:
      */
     enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
 
@@ -91,49 +85,36 @@ namespace TrilinosWrappers
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format and puts
-       * the restart parameter in case the
-       * derived class is GMRES.
+       * Sets the additional data field to the desired output format and puts
+       * the restart parameter in case the derived class is GMRES.
        *
-       * TODO: Find a better way for
-       * setting the GMRES restart
-       * parameter since it is quite
-       * inelegant to set a specific option
-       * of one solver in the base class
-       * for all solvers.
+       * TODO: Find a better way for setting the GMRES restart parameter since
+       * it is quite inelegant to set a specific option of one solver in the
+       * base class for all solvers.
        */
       AdditionalData (const bool         output_solver_details   = false,
                       const unsigned int gmres_restart_parameter = 30);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       const bool output_solver_details;
 
       /**
-       * Restart parameter for GMRES
-       * solver.
+       * Restart parameter for GMRES solver.
        */
       const unsigned int gmres_restart_parameter;
     };
 
     /**
-     * Constructor. Takes the
-     * solver control object and
-     * creates the solver.
+     * Constructor. Takes the solver control object and creates the solver.
      */
     SolverBase (SolverControl  &cn);
 
     /**
-     * Second constructor. This
-     * constructor takes an enum
-     * object that specifies the
-     * solver name and sets the
-     * appropriate Krylov
-     * method.
+     * Second constructor. This constructor takes an enum object that
+     * specifies the solver name and sets the appropriate Krylov method.
      */
     SolverBase (const enum SolverName  solver_name,
                 SolverControl         &cn);
@@ -144,15 +125,9 @@ namespace TrilinosWrappers
     virtual ~SolverBase ();
 
     /**
-     * Solve the linear system
-     * <tt>Ax=b</tt>. 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.
+     * Solve the linear system <tt>Ax=b</tt>. 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 (const SparseMatrix     &A,
@@ -161,18 +136,11 @@ namespace TrilinosWrappers
            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.
+     * 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,
@@ -181,22 +149,14 @@ namespace TrilinosWrappers
            const PreconditionBase &preconditioner);
 
     /**
-     * Solve the linear system
-     * <tt>Ax=b</tt>. 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.
+     * Solve the linear system <tt>Ax=b</tt>. 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 (const SparseMatrix           &A,
@@ -205,24 +165,15 @@ namespace TrilinosWrappers
            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.
+     * Solve the linear system <tt>Ax=b</tt> where <tt>A</tt> is an
+     * operator. This function can be used for matrix free
+     * computations. 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,
@@ -231,8 +182,33 @@ namespace TrilinosWrappers
            const PreconditionBase       &preconditioner);
 
     /**
-     * Access to object that controls
-     * convergence.
+     * Solve the linear system <tt>Ax=b</tt> for deal.II's parallel
+     * distributed vectors. 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 (const SparseMatrix                                  &A,
+           dealii::parallel::distributed::Vector<double>       &x,
+           const dealii::parallel::distributed::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 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,
+           dealii::parallel::distributed::Vector<double>       &x,
+           const dealii::parallel::distributed::Vector<double> &b,
+           const PreconditionBase                              &preconditioner);
+
+
+    /**
+     * Access to object that controls convergence.
      */
     SolverControl &control() const;
 
@@ -247,16 +223,10 @@ namespace TrilinosWrappers
   protected:
 
     /**
-     * Reference to the object that
-     * controls convergence of the
-     * iterative solver. In fact,
-     * for these Trilinos wrappers,
-     * Trilinos does so itself, but
-     * we copy the data from this
-     * object before starting the
-     * solution process, and copy
-     * the data back into it
-     * afterwards.
+     * Reference to the object that controls convergence of the iterative
+     * solver. In fact, for these Trilinos wrappers, Trilinos does so itself,
+     * but we copy the data from this object before starting the solution
+     * process, and copy the data back into it afterwards.
      */
     SolverControl &solver_control;
 
@@ -266,28 +236,23 @@ namespace TrilinosWrappers
      * 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);
+    void do_solve(const PreconditionBase &preconditioner);
 
     /**
-     * A structure that collects
-     * the Trilinos sparse matrix,
-     * the right hand side vector
-     * and the solution vector,
-     * which is passed down to the
+     * A structure that collects the Trilinos sparse matrix, the right hand
+     * side vector and the solution vector, which is passed down to the
      * Trilinos solver.
      */
     std_cxx1x::shared_ptr<Epetra_LinearProblem> linear_problem;
 
     /**
-     * A structure that contains
-     * the Trilinos solver and
-     * preconditioner objects.
+     * A structure that contains the Trilinos solver and preconditioner
+     * objects.
      */
     AztecOO solver;
 
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
 
@@ -306,35 +271,28 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
 
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format.
+       * Sets the additional data field to the desired output format.
        */
       AdditionalData (const bool output_solver_details = false);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
     };
 
     /**
-     * Constructor. In contrast to
-     * deal.II's own solvers, there is no
-     * need to give a vector memory
-     * object.
+     * Constructor. In contrast to deal.II's own solvers, there is no need to
+     * give a vector memory object.
      *
-     * The last argument takes a structure
-     * with additional, solver dependent
+     * The last argument takes a structure with additional, solver dependent
      * flags for tuning.
      */
     SolverCG (SolverControl        &cn,
@@ -342,8 +300,7 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
   };
@@ -361,34 +318,27 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format.
+       * Sets the additional data field to the desired output format.
        */
       AdditionalData (const bool output_solver_details = false);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
     };
 
     /**
-     * Constructor. In contrast to
-     * deal.II's own solvers, there is no
-     * need to give a vector memory
-     * object.
+     * Constructor. In contrast to deal.II's own solvers, there is no need to
+     * give a vector memory object.
      *
-     * The last argument takes a structure
-     * with additional, solver dependent
+     * The last argument takes a structure with additional, solver dependent
      * flags for tuning.
      */
     SolverCGS (SolverControl        &cn,
@@ -396,8 +346,7 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
   };
@@ -414,43 +363,34 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
     struct AdditionalData
     {
       /**
-       * Constructor. By default, set the
-       * number of temporary vectors to
-       * 30, i.e. do a restart every 30
-       * iterations.
+       * Constructor. By default, set the number of temporary vectors to 30,
+       * i.e. do a restart every 30 iterations.
        */
       AdditionalData (const bool         output_solver_details = false,
                       const unsigned int restart_parameter = 30);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
 
       /**
-       * Maximum number of
-       * tmp vectors.
+       * Maximum number of tmp vectors.
        */
       unsigned int restart_parameter;
     };
 
     /**
-     * Constructor. In contrast to
-     * deal.II's own solvers, there is no
-     * need to give a vector memory
-     * object.
+     * Constructor. In contrast to deal.II's own solvers, there is no need to
+     * give a vector memory object.
      *
-     * The last argument takes a structure
-     * with additional, solver dependent
+     * The last argument takes a structure with additional, solver dependent
      * flags for tuning.
      */
     SolverGMRES (SolverControl        &cn,
@@ -458,8 +398,7 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
   };
@@ -477,34 +416,27 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format.
+       * Sets the additional data field to the desired output format.
        */
       AdditionalData (const bool output_solver_details = false);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
     };
 
     /**
-     * Constructor. In contrast to
-     * deal.II's own solvers, there is no
-     * need to give a vector memory
-     * object.
+     * Constructor. In contrast to deal.II's own solvers, there is no need to
+     * give a vector memory object.
      *
-     * The last argument takes a structure
-     * with additional, solver dependent
+     * The last argument takes a structure with additional, solver dependent
      * flags for tuning.
      */
     SolverBicgstab (SolverControl        &cn,
@@ -512,8 +444,7 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
   };
@@ -531,34 +462,27 @@ namespace TrilinosWrappers
   {
   public:
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format.
+       * Sets the additional data field to the desired output format.
        */
       AdditionalData (const bool output_solver_details = false);
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
     };
 
     /**
-     * Constructor. In contrast to
-     * deal.II's own solvers, there is no
-     * need to give a vector memory
-     * object.
+     * Constructor. In contrast to deal.II's own solvers, there is no need to
+     * give a vector memory object.
      *
-     * The last argument takes a structure
-     * with additional, solver dependent
+     * The last argument takes a structure with additional, solver dependent
      * flags for tuning.
      */
     SolverTFQMR (SolverControl        &cn,
@@ -566,8 +490,7 @@ namespace TrilinosWrappers
 
   protected:
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
   };
@@ -586,31 +509,46 @@ namespace TrilinosWrappers
   public:
 
     /**
-     * Standardized data struct to
-     * pipe additional data to the
-     * solver.
+     * Standardized data struct to pipe additional data to the solver.
      */
 
     struct AdditionalData
     {
       /**
-       * Sets the additional data field to
-       * the desired output format.
+       * Sets the additional data field to the desired output format.
        */
-      AdditionalData (const bool output_solver_details = false);
+      AdditionalData (const bool output_solver_details = false,
+                      const std::string &solver_type = "Amesos_Klu");
 
       /**
-       * Enables/disables the output of
-       * solver details (residual in each
+       * Enables/disables the output of solver details (residual in each
        * iterations etc.).
        */
       bool output_solver_details;
+
+      /**
+       * Set the solver type (for third party solver support of Trilinos
+       * Amesos package). Possibilities are:
+       * <ol>
+       *   <li>  "Amesos_Lapack" </li>
+       *   <li>  "Amesos_Scalapack" </li>
+       *   <li>  "Amesos_Klu" </li>
+       *   <li>  "Amesos_Umfpack" </li>
+       *   <li>  "Amesos_Pardiso" </li>
+       *   <li>  "Amesos_Taucs" </li>
+       *   <li>  "Amesos_Superlu" </li>
+       *   <li>  "Amesos_Superludist" </li>
+       *   <li>  "Amesos_Dscpack" </li>
+       *   <li>  "Amesos_Mumps" </li>
+       * </ol>
+       * Note that the availability of these solvers in deal.II depends on
+       * which solvers were set when configuring Trilinos.
+       */
+      std::string solver_type;
     };
 
     /**
-     * Constructor. Takes the
-     * solver control object and
-     * creates the solver.
+     * Constructor. Takes the solver control object and creates the solver.
      */
     SolverDirect (SolverControl  &cn,
                   const AdditionalData &data = AdditionalData());
@@ -621,12 +559,10 @@ namespace TrilinosWrappers
     virtual ~SolverDirect ();
 
     /**
-     * Solve the linear system
-     * <tt>Ax=b</tt>. Creates a KLU
-     * factorization of the matrix and
-     * performs the solve. Note that
-     * there is no need for a
-     * preconditioner here.
+     * Solve the linear system <tt>Ax=b</tt>. Creates a factorization of the
+     * matrix with the package chosen from the additional data structure and
+     * performs the solve. Note that there is no need for a preconditioner
+     * here.
      */
     void
     solve (const SparseMatrix     &A,
@@ -634,28 +570,28 @@ namespace TrilinosWrappers
            const VectorBase       &b);
 
     /**
-     * Solve the linear system
-     * <tt>Ax=b</tt>. 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.
+     * Solve the linear system <tt>Ax=b</tt>. This class works with Trilinos
+     * matrices, but takes deal.II serial vectors as argument. Since these
+     * vectors are not distributed, this function does only what you expect in
+     * case the matrix is serial (i.e., locally owned). Otherwise, an
+     * exception will be thrown.
      */
     void
     solve (const SparseMatrix           &A,
            dealii::Vector<double>       &x,
            const dealii::Vector<double> &b);
 
+    /**
+     * Solve the linear system <tt>Ax=b</tt> for deal.II's own parallel
+     * vectors. Creates a factorization of the matrix with the package chosen
+     * from the additional data structure and performs the solve. Note that
+     * there is no need for a preconditioner here.
+     */
+    void
+    solve (const SparseMatrix                                  &A,
+           dealii::parallel::distributed::Vector<double>       &x,
+           const dealii::parallel::distributed::Vector<double> &b);
+
     /**
      * Access to object that controls
      * convergence.
@@ -671,44 +607,37 @@ namespace TrilinosWrappers
                     << " occurred while calling a Trilinos function");
 
   private:
+    /**
+     * Actually performs the operations for solving the linear system,
+     * including the factorization and forward and backward substitution.
+     */
+    void do_solve();
 
     /**
-     * Reference to the object that
-     * controls convergence of the
-     * iterative solver. In fact,
-     * for these Trilinos wrappers,
-     * Trilinos does so itself, but
-     * we copy the data from this
-     * object before starting the
-     * solution process, and copy
-     * the data back into it
-     * afterwards.
+     * Reference to the object that controls convergence of the iterative
+     * solver. In fact, for these Trilinos wrappers, Trilinos does so itself,
+     * but we copy the data from this object before starting the solution
+     * process, and copy the data back into it afterwards.
      */
     SolverControl &solver_control;
 
     /**
-     * A structure that collects
-     * the Trilinos sparse matrix,
-     * the right hand side vector
-     * and the solution vector,
-     * which is passed down to the
+     * A structure that collects the Trilinos sparse matrix, the right hand
+     * side vector and the solution vector, which is passed down to the
      * Trilinos solver.
      */
     std_cxx1x::shared_ptr<Epetra_LinearProblem> linear_problem;
 
     /**
-     * A structure that contains
-     * the Trilinos solver and
-     * preconditioner objects.
+     * A structure that contains the Trilinos solver and preconditioner
+     * objects.
      */
     std_cxx1x::shared_ptr<Amesos_BaseSolver> solver;
 
     /**
-     * Store a copy of the flags for this
-     * particular solver.
+     * Store a copy of the flags for this particular solver.
      */
     const AdditionalData additional_data;
-
   };
 
 }
index 7d206765fdaee6151d3f35e8dc3a9d8b714dd202..f7a567fe8e52d96bec82e4c5ac78321ffd7a2682 100644 (file)
@@ -77,17 +77,14 @@ namespace TrilinosWrappers
   {
     linear_problem.reset();
 
-    // We need an
-    // Epetra_LinearProblem object
-    // to let the AztecOO solver
-    // know about the matrix and
-    // vectors.
+    // We need an Epetra_LinearProblem object to let the AztecOO solver know
+    // about the matrix and vectors.
     linear_problem.reset
     (new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
                               &x.trilinos_vector(),
                               const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
 
-    execute_solve(preconditioner);
+    do_solve(preconditioner);
   }
 
 
@@ -100,17 +97,14 @@ namespace TrilinosWrappers
   {
     linear_problem.reset();
 
-    // We need an
-    // Epetra_LinearProblem object
-    // to let the AztecOO solver
-    // know about the matrix and
-    // vectors.
+    // 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())));
 
-    execute_solve(preconditioner);
+    do_solve(preconditioner);
   }
 
 
@@ -123,9 +117,8 @@ namespace TrilinosWrappers
   {
     linear_problem.reset();
 
-    // In case we call the solver with
-    // deal.II vectors, we create views
-    // of the vectors in Epetra format.
+    // In case we call the solver with deal.II vectors, we create views of the
+    // vectors in Epetra format.
     Assert (x.size() == A.n(),
             ExcDimensionMismatch(x.size(), A.n()));
     Assert (b.size() == A.m(),
@@ -138,16 +131,13 @@ namespace TrilinosWrappers
     Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
     Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
 
-    // We need an
-    // Epetra_LinearProblem object
-    // to let the AztecOO solver
-    // know about the matrix and
-    // vectors.
+    // We need an Epetra_LinearProblem object to let the AztecOO solver know
+    // about the matrix and vectors.
     linear_problem.reset (new Epetra_LinearProblem
                           (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
                            &ep_x, &ep_b));
 
-    execute_solve(preconditioner);
+    do_solve(preconditioner);
   }
 
 
@@ -163,29 +153,78 @@ namespace TrilinosWrappers
     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.
+    // 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);
+    do_solve(preconditioner);
   }
 
 
 
   void
-  SolverBase::execute_solve(const PreconditionBase &preconditioner)
+  SolverBase::solve (const SparseMatrix                                  &A,
+                     dealii::parallel::distributed::Vector<double>       &x,
+                     const dealii::parallel::distributed::Vector<double> &b,
+                     const PreconditionBase                              &preconditioner)
+  {
+    linear_problem.reset();
+
+    // In case we call the solver with deal.II vectors, we create views of the
+    // vectors in Epetra format.
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+                     A.domain_partitioner().NumMyElements());
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+                     A.range_partitioner().NumMyElements());
+
+    Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
+    Epetra_Vector ep_b (View, A.range_partitioner(), 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
+                          (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+                           &ep_x, &ep_b));
+
+    do_solve(preconditioner);
+  }
+
+
+
+  void
+  SolverBase::solve (Epetra_Operator                                     &A,
+                     dealii::parallel::distributed::Vector<double>       &x,
+                     const dealii::parallel::distributed::Vector<double> &b,
+                     const PreconditionBase                              &preconditioner)
+  {
+    linear_problem.reset();
+
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+                     A.OperatorDomainMap().NumMyElements());
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+                     A.OperatorRangeMap().NumMyElements());
+
+    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));
+
+    do_solve(preconditioner);
+  }
+
+
+
+  void
+  SolverBase::do_solve(const PreconditionBase &preconditioner)
   {
     int ierr;
 
-    // Next we can allocate the
-    // AztecOO solver...
+    // Next we can allocate the AztecOO solver...
     solver.SetProblem(*linear_problem);
 
-    // ... and we can specify the
-    // solver to be used.
+    // ... and we can specify the solver to be used.
     switch (solver_name)
       {
       case cg:
@@ -208,8 +247,7 @@ namespace TrilinosWrappers
         Assert (false, ExcNotImplemented());
       }
 
-    // Introduce the preconditioner,
-    // if the identity preconditioner is used,
+    // Introduce the preconditioner, if the identity preconditioner is used,
     // the precondioner is set to none, ...
     if (preconditioner.preconditioner.use_count()!=0)
       {
@@ -229,13 +267,9 @@ namespace TrilinosWrappers
     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
+    // 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:
@@ -249,18 +283,14 @@ namespace TrilinosWrappers
                                        "loss of precision"));
       case -4:
         AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: "
-                                       "GMRES hessenberg ill-conditioned"));
+                                       "GMRES Hessenberg ill-conditioned"));
       default:
         AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
       }
 
-    // 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.
+    // 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());
 
     if (solver_control.last_check() != SolverControl::success)
@@ -381,9 +411,11 @@ namespace TrilinosWrappers
   /* ---------------------- SolverDirect ------------------------ */
 
   SolverDirect::AdditionalData::
-  AdditionalData (const bool output_solver_details)
+  AdditionalData (const bool output_solver_details,
+                  const std::string &solver_type)
     :
-    output_solver_details (output_solver_details)
+    output_solver_details (output_solver_details),
+    solver_type(solver_type)
   {}
 
 
@@ -393,7 +425,7 @@ namespace TrilinosWrappers
                               const AdditionalData &data)
     :
     solver_control (cn),
-    additional_data (data.output_solver_details)
+    additional_data (data.output_solver_details,data.solver_type)
   {}
 
 
@@ -412,49 +444,51 @@ namespace TrilinosWrappers
 
 
   void
-  SolverDirect::solve (const SparseMatrix     &A,
-                       VectorBase             &x,
-                       const VectorBase       &b)
+  SolverDirect::do_solve()
   {
-    // First set whether we want to print
-    // the solver information to screen
-    // or not.
+    // Fetch return value of Amesos Solver functions
+    int ierr;
+
+    // First set whether we want to print the solver information to screen or
+    // not.
     ConditionalOStream  verbose_cout (std::cout,
                                       additional_data.output_solver_details);
 
-    linear_problem.reset();
     solver.reset();
 
-    // We need an
-    // Epetra_LinearProblem object
-    // to let the AztecOO solver
-    // know about the matrix and
-    // vectors.
-    linear_problem.reset
-    (new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
-                              &x.trilinos_vector(),
-                              const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
+    // Next allocate the Amesos solver, this is done in two steps, first we
+    // create a solver Factory and and generate with that the concrete Amesos
+    // solver, if possible.
+    Amesos Factory;
+
+    AssertThrow(
+      Factory.Query(additional_data.solver_type.c_str()),
+      ExcMessage (std::string ("You tried to select the solver type <") +
+                      additional_data.solver_type +
+                      "> but this solver is not supported by Trilinos either "
+                      "because it does not exist, or because Trilinos was not "
+                      "configured for its use.")
+    );
 
-    // Next we can allocate the
-    // AztecOO solver...
-    solver.reset (Amesos().Create("Amesos_Klu", *linear_problem));
+    solver.reset (
+      Factory.Create(additional_data.solver_type.c_str(), *linear_problem)
+    );
 
     verbose_cout << "Starting symbolic factorization" << std::endl;
-    solver->SymbolicFactorization();
+    ierr = solver->SymbolicFactorization();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
     verbose_cout << "Starting numeric factorization" << std::endl;
-    solver->NumericFactorization();
+    ierr = solver->NumericFactorization();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
     verbose_cout << "Starting solve" << std::endl;
-    solver->Solve();
-
-    // 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.
+    ierr = solver->Solve();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    // 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 (0, 0);
 
     if (solver_control.last_check() != SolverControl::success)
@@ -463,25 +497,31 @@ namespace TrilinosWrappers
   }
 
 
+  void
+  SolverDirect::solve (const SparseMatrix     &A,
+                       VectorBase             &x,
+                       const VectorBase       &b)
+  {
+    // We need an Epetra_LinearProblem object to let the Amesos solver know
+    // about the matrix and vectors.
+    linear_problem.reset
+    (new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+                              &x.trilinos_vector(),
+                              const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
+
+    do_solve();
+  }
+
+
 
   void
   SolverDirect::solve (const SparseMatrix           &A,
                        dealii::Vector<double>       &x,
                        const dealii::Vector<double> &b)
   {
-    // First set whether we want to print
-    // the solver information to screen
-    // or not.
-    ConditionalOStream  verbose_cout (std::cout,
-                                      additional_data.output_solver_details);
-
-    linear_problem.reset();
-    solver.reset();
-
 
-    // In case we call the solver with
-    // deal.II vectors, we create views
-    // of the vectors in Epetra format.
+    // In case we call the solver with deal.II vectors, we create views of the
+    // vectors in Epetra format.
     Assert (x.size() == A.n(),
             ExcDimensionMismatch(x.size(), A.n()));
     Assert (b.size() == A.m(),
@@ -491,44 +531,37 @@ namespace TrilinosWrappers
     Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
     Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
 
-    // We need an
-    // Epetra_LinearProblem object
-    // to let the AztecOO solver
-    // know about the matrix and
-    // vectors.
+    // We need an Epetra_LinearProblem object to let the Amesos solver know
+    // about the matrix and vectors.
     linear_problem.reset (new Epetra_LinearProblem
                           (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
                            &ep_x, &ep_b));
 
-    // Next we can allocate the
-    // AztecOO solver...
-    solver.reset (Amesos().Create("Amesos_Klu", *linear_problem));
+    do_solve();
+  }
 
-    verbose_cout << "Starting symbolic factorization" << std::endl;
-    solver->SymbolicFactorization();
 
-    verbose_cout << "Starting numeric factorization" << std::endl;
-    solver->NumericFactorization();
-
-    verbose_cout << "Starting solve" << std::endl;
-    solver->Solve();
-
-    // 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 (0, 0);
-
-    if (solver_control.last_check() != SolverControl::success)
-      AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(),
-                                                       solver_control.last_value()));
-  }
 
+  void
+  SolverDirect::solve (const SparseMatrix                                  &A,
+                       dealii::parallel::distributed::Vector<double>       &x,
+                       const dealii::parallel::distributed::Vector<double> &b)
+  {
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+                     A.domain_partitioner().NumMyElements());
+    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+                     A.range_partitioner().NumMyElements());
+    Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
+    Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
 
+    // We need an Epetra_LinearProblem object to let the Amesos solver know
+    // about the matrix and vectors.
+    linear_problem.reset (new Epetra_LinearProblem
+                          (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+                           &ep_x, &ep_b));
 
+    do_solve();
+  }
 
 }
 
diff --git a/tests/trilinos/direct_solver.cc b/tests/trilinos/direct_solver.cc
new file mode 100644 (file)
index 0000000..1c24319
--- /dev/null
@@ -0,0 +1,332 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests Trilinos direct solvers on a 2D Poisson equation for linear elements
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/base/function.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4 ();
+  void run ();
+
+private:
+  void make_grid ();
+  void setup_system();
+  void assemble_system ();
+  void solve ();
+
+  Triangulation<dim>   triangulation;
+  FE_Q<dim>            fe;
+  DoFHandler<dim>      dof_handler;
+
+  ConstraintMatrix     constraints;
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                  const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    return_value += 4*std::pow(p(i), 4);
+
+  return return_value;
+}
+
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                   const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4 ()
+  :
+  fe (1),
+  dof_handler (triangulation)
+{}
+
+
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (6);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  constraints.clear();
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close();
+
+  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false);
+  system_matrix.reinit (c_sparsity);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void Step4<dim>::assemble_system ()
+{
+  QGauss<dim>  quadrature_formula(fe.degree+1);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
+                                   fe_values.shape_grad (j, q_point) *
+                                   fe_values.JxW (q_point));
+
+            cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                            right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                            fe_values.JxW (q_point));
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+                                             local_dof_indices,
+                                             system_matrix, system_rhs);
+    }
+  system_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void Step4<dim>::solve ()
+{
+  // Compute 'reference' solution with CG solver and SSOR preconditioner
+  TrilinosWrappers::PreconditionSSOR preconditioner;
+  solution = 0;
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              solver (solver_control);
+  preconditioner.initialize(system_matrix);
+  solver.solve (system_matrix, solution, system_rhs,
+                preconditioner);
+
+  Vector<double> output(solution);
+  {
+    deallog.push("DirectKLU");
+    TrilinosWrappers::SolverDirect::AdditionalData data;
+    data.solver_type = "Amesos_Klu";
+    SolverControl           solver_control (1000, 1e-10);
+    TrilinosWrappers::SolverDirect solver(solver_control, data);
+    solver.solve (system_matrix, output, system_rhs);
+    output -= solution;
+    deallog << "Norm of error in direct solve: " << output.l2_norm()
+            << std::endl;
+    deallog.pop();
+  }
+
+  {
+    deallog.push("DirectLAPACK");
+    TrilinosWrappers::SolverDirect::AdditionalData data;
+    data.solver_type = "Amesos_Lapack";
+    SolverControl           solver_control (1000, 1e-10);
+    TrilinosWrappers::SolverDirect solver(solver_control, data);
+    solver.solve (system_matrix, output, system_rhs);
+    output -= solution;
+    deallog << "Norm of error in direct solve: " << output.l2_norm()
+            << std::endl;
+    deallog.pop();
+  }
+
+  {
+    deallog.push("DirectDscpack");
+    TrilinosWrappers::SolverDirect::AdditionalData data;
+    data.solver_type = "Amesos_Dscpack";
+    SolverControl           solver_control (1000, 1e-10);
+    TrilinosWrappers::SolverDirect solver(solver_control, data);
+    try
+      {
+        solver.solve (system_matrix, output, system_rhs);
+      }
+    catch (dealii::ExceptionBase &exc)
+      {
+        deallog << "Error: " << std::endl;
+        exc.print_info(deallog.get_file_stream());
+      }
+    deallog.pop();
+  }
+
+  {
+    deallog.push("Dummy");
+    TrilinosWrappers::SolverDirect::AdditionalData data;
+    data.solver_type = "DummySolver";
+    SolverControl           solver_control (1000, 1e-10);
+    TrilinosWrappers::SolverDirect solver(solver_control, data);
+    try
+      {
+        solver.solve (system_matrix, output, system_rhs);
+      }
+    catch (dealii::ExceptionBase &exc)
+      {
+        deallog << "Error: " << std::endl;
+        exc.print_info(deallog.get_file_stream());
+      }
+    deallog.pop();
+  }
+}
+
+
+
+template <int dim>
+void Step4<dim>::run()
+{
+  make_grid();
+  setup_system();
+  assemble_system();
+  solve();
+}
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  try
+    {
+      Step4<2> test;
+      test.run();
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/direct_solver.output b/tests/trilinos/direct_solver.output
new file mode 100644 (file)
index 0000000..9b38539
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:cg::Starting value 21.8299
+DEAL:cg::Convergence step 86 value 0
+DEAL:DirectKLU::Starting value 0
+DEAL:DirectKLU::Convergence step 0 value 0
+DEAL:DirectKLU::Norm of error in direct solve: 0
+DEAL:DirectLAPACK::Starting value 0
+DEAL:DirectLAPACK::Convergence step 0 value 0
+DEAL:DirectLAPACK::Norm of error in direct solve: 0
+DEAL:DirectDscpack::Error: 
+You tried to select the solver type <Amesos_Dscpack> but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use.
+DEAL:Dummy::Error: 
+You tried to select the solver type <DummySolver> but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use.

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.