From: Martin Kronbichler Date: Wed, 9 Jul 2014 08:01:15 +0000 (+0000) Subject: Patch by Uwe Koecher: select among Trilinos solvers X-Git-Tag: v8.2.0-rc1~331 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9f3a0d2f9d71075d8dedd13d2cbdb92193b26eaf;p=dealii.git Patch by Uwe Koecher: select among Trilinos solvers git-svn-id: https://svn.dealii.org/trunk@33117 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 974eef62d3..635384ed35 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,13 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Uwe Köcher, Martin Kronbichler, 2014/07/09) +
    2. +
  2. New: There are now three new preconditioner classes TrilinosWrappers::PreconditionBlockJacobi, TrilinosWrappers::PreconditionBlockSSOR, and diff --git a/deal.II/include/deal.II/lac/trilinos_precondition.h b/deal.II/include/deal.II/lac/trilinos_precondition.h index 01806b9646..c1ceb62471 100644 --- a/deal.II/include/deal.II/lac/trilinos_precondition.h +++ b/deal.II/include/deal.II/lac/trilinos_precondition.h @@ -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" + *
      + *
    1. "Aztec"
      2. + *
    2. "IFPACK"
      3. + *
    3. "Jacobi"
      4. + *
    4. "ML symmetric Gauss-Seidel"
      5. + *
    5. "symmetric Gauss-Seidel"
      6. + *
    6. "ML Gauss-Seidel"
      7. + *
    7. "Gauss-Seidel"
      8. + *
    8. "block Gauss-Seidel"
      9. + *
    9. "symmetric block Gauss-Seidel"
      10. + *
    10. "Chebyshev"
      11. + *
    11. "MLS"
      12. + *
    12. "Hiptmair"
      13. + *
    13. "Amesos-KLU"
      14. + *
    14. "Amesos-Superlu"
      15. + *
    15. "Amesos-UMFPACK"
      16. + *
    16. "Amesos-Superludist"
      17. + *
    17. "Amesos-MUMPS"
      18. + *
    18. "user-defined"
      19. + *
    19. "SuperLU"
      20. + *
    20. "IFPACK-Chebyshev"
      21. + *
    21. "self"
      22. + *
    22. "do-nothing"
      23. + *
    23. "IC"
      24. + *
    24. "ICT"
      25. + *
    25. "ILU"
      26. + *
    26. "ILUT"
      27. + *
    27. "Block Chebyshev"
      28. + *
    28. "IFPACK-Block Chebyshev"
      29. + *
    */ const char* smoother_type; diff --git a/deal.II/include/deal.II/lac/trilinos_solver.h b/deal.II/include/deal.II/lac/trilinos_solver.h index 1c4ead78a8..9ccc910391 100644 --- a/deal.II/include/deal.II/lac/trilinos_solver.h +++ b/deal.II/include/deal.II/lac/trilinos_solver.h @@ -26,6 +26,7 @@ # include # include # include +# include # include # include @@ -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 - * Ax=b. 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 Ax=b. 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 - * Ax=b where A - * 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 Ax=b where A 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 - * Ax=b. 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 Ax=b. 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 - * Ax=b where A - * 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 Ax=b where A 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 Ax=b 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 &x, + const dealii::parallel::distributed::Vector &b, + const PreconditionBase &preconditioner); + + /** + * Solve the linear system Ax=b where A 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 &x, + const dealii::parallel::distributed::Vector &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 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: + *
      + *
    1. "Amesos_Lapack"
      2. + *
    2. "Amesos_Scalapack"
      3. + *
    3. "Amesos_Klu"
      4. + *
    4. "Amesos_Umfpack"
      5. + *
    5. "Amesos_Pardiso"
      6. + *
    6. "Amesos_Taucs"
      7. + *
    7. "Amesos_Superlu"
      8. + *
    8. "Amesos_Superludist"
      9. + *
    9. "Amesos_Dscpack"
      10. + *
    10. "Amesos_Mumps"
      11. + *
    + * 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 - * Ax=b. 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 Ax=b. 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 - * Ax=b. 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 Ax=b. 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 &x, const dealii::Vector &b); + /** + * Solve the linear system Ax=b 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 &x, + const dealii::parallel::distributed::Vector &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 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 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; - }; } diff --git a/deal.II/source/lac/trilinos_solver.cc b/deal.II/source/lac/trilinos_solver.cc index 7d206765fd..f7a567fe8e 100644 --- a/deal.II/source/lac/trilinos_solver.cc +++ b/deal.II/source/lac/trilinos_solver.cc @@ -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(&A.trilinos_matrix()), &x.trilinos_vector(), const_cast(&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(&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(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(&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(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 &x, + const dealii::parallel::distributed::Vector &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(x.local_size()), + A.domain_partitioner().NumMyElements()); + AssertDimension (static_cast(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(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(&A.trilinos_matrix()), + &ep_x, &ep_b)); + + do_solve(preconditioner); + } + + + + void + SolverBase::solve (Epetra_Operator &A, + dealii::parallel::distributed::Vector &x, + const dealii::parallel::distributed::Vector &b, + const PreconditionBase &preconditioner) + { + linear_problem.reset(); + + AssertDimension (static_cast(x.local_size()), + A.OperatorDomainMap().NumMyElements()); + AssertDimension (static_cast(b.local_size()), + A.OperatorRangeMap().NumMyElements()); + + Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.begin()); + Epetra_Vector ep_b (View, A.OperatorRangeMap(), const_cast(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(&A.trilinos_matrix()), - &x.trilinos_vector(), - const_cast(&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(&A.trilinos_matrix()), + &x.trilinos_vector(), + const_cast(&b.trilinos_vector()))); + + do_solve(); + } + + void SolverDirect::solve (const SparseMatrix &A, dealii::Vector &x, const dealii::Vector &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(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(&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 &x, + const dealii::parallel::distributed::Vector &b) + { + AssertDimension (static_cast(x.local_size()), + A.domain_partitioner().NumMyElements()); + AssertDimension (static_cast(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(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(&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 index 0000000000..1c243190b2 --- /dev/null +++ b/tests/trilinos/direct_solver.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + + +template +double RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4 () + : + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (6); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + 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 +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + FEValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::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_pointget_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 +void Step4::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 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 +void Step4::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 index 0000000000..9b38539243 --- /dev/null +++ b/tests/trilinos/direct_solver.output @@ -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 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 but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use.