From 6b9b80f6e74508a9d574e87c5ea45cbeff43b560 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 5 Jan 2009 07:38:24 +0000 Subject: [PATCH] Provide some more Trilinos preconditioners and solvers. git-svn-id: https://svn.dealii.org/trunk@18069 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_block_vector.h | 75 ++--- .../lac/include/lac/trilinos_precondition.h | 182 ++++++++++ deal.II/lac/include/lac/trilinos_solver.h | 279 +++++++++++++++- .../lac/include/lac/trilinos_sparse_matrix.h | 6 +- .../lac/include/lac/trilinos_vector_base.h | 12 +- deal.II/lac/source/trilinos_block_vector.cc | 2 +- deal.II/lac/source/trilinos_precondition.cc | 54 +++ deal.II/lac/source/trilinos_solver.cc | 316 +++++++++++++++++- 8 files changed, 843 insertions(+), 83 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index 13e27772f5..cf313ad4fe 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -237,46 +237,41 @@ namespace TrilinosWrappers */ void reinit (const unsigned int num_blocks); - /** - * This reinit function is - * meant to be used for - * parallel calculations where - * some non-local data has to - * be used. The typical - * situation where one needs - * this function is the call of - * the - * FEValues::get_function_values - * function (or of some - * derivatives) in - * parallel. Since it is - * usually faster to retrieve - * the data in advance, this - * function can be called - * before the assembly forks - * out to the different - * processors. What this - * function does is the - * following: It takes the - * information in the columns - * of the given matrix and - * looks which data couples - * between the different - * processors. That data is - * then queried from the input - * vector. Note that you should - * not write to the resulting - * vector any more, since the - * some data can be stored - * several times on different - * processors, leading to - * unpredictable results. In - * particular, such a vector - * cannot be used for - * matrix-vector products as - * for example done during the - * solution of linear systems. - */ + /** + * This reinit function is meant to + * be used for parallel + * calculations where some + * non-local data has to be + * used. The typical situation + * where one needs this function is + * the call of the + * FEValues::get_function_values + * function (or of some + * derivatives) in parallel. Since + * it is usually faster to retrieve + * the data in advance, this + * function can be called before + * the assembly forks out to the + * different processors. What this + * function does is the following: + * It takes the information in the + * columns of the given matrix and + * looks which data couples between + * the different processors. That + * data is then queried from the + * input vector. Note that you + * should not write to the + * resulting vector any more, since + * the some data can be stored + * several times on different + * processors, leading to + * unpredictable results. In + * particular, such a vector cannot + * be used for matrix-vector + * products as for example done + * during the solution of linear + * systems. + */ void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m, const BlockVector &v); diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index 04101161ad..dde9daeff0 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -882,6 +882,188 @@ namespace TrilinosWrappers + + + +/** + * A wrapper class for a thresholded incomplete LU factorization (ILU-T) + * preconditioner for Trilinos matrices. This preconditioner works both in + * serial and in parallel, depending on the matrix it is based on. In + * general, an incomplete factorization does not take all fill-in elements + * that would appear in a full factorization (that is the basis for a direct + * solve). For the ILU-T precondtioner, the parameter ilut_drop + * lets the user specify which elements should be dropped (i.e., should not + * be part of the incomplete decomposition). Trilinos calculates first the + * complete factorization for one row, and then skips those elements that + * are lower than the threshold. This is the main difference to the + * non-thresholded ILU preconditioner, where the parameter + * ilut_fill governs the incomplete factorization structure. This + * parameter is available here as well, but provides only some extra + * information here. + * + * The AdditionalData data structure allows to set preconditioner + * options. Besides the fill-in arguments, these options are some options + * for perturbations (see the documentation of the AdditionalData structure + * for details), and a parameter overlap that determines if and how + * much overlap there should be between the matrix partitions on the various + * MPI processes. The default settings are 0 for the additional fill-in, 0 + * for the absolute augmentation tolerance, 1 for the relative augmentation + * tolerance, 0 for the overlap. + * + * Note that a parallel applicatoin of the ILU-T preconditioner is + * actually a block-Jacobi preconditioner with block size equal to the + * local matrix size. Spoken more technically, this parallel operation + * is an additive + * Schwarz method with an ILU approximate solve as inner + * solver, based on the (outer) parallel partitioning. + * + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2009 + */ + class PreconditionILUT : public PreconditionBase + { + public: + /** + * Standardized data struct to pipe + * additional parameters to the + * preconditioner. The Trilinos ILU-T + * decomposition allows for some + * fill-in, so it actually is a + * threshold incomplete LU + * factorization. The amount of + * fill-in, and hence, the amount of + * memory used by this + * preconditioner, is controlled by + * the parameters ilut_drop + * and ilut_fill, which + * specifies a threshold about which + * values should form the incomplete + * factorization and the level of + * additional fill-in. When forming + * the preconditioner, for certain + * problems bad conditioning (or just + * bad luck) can cause the + * preconditioner to be very poorly + * conditioned. Hence it can help to + * add diagonal perturbations to the + * original matrix and form the + * preconditioner for this slightly + * better matrix. ilut_atol + * is an absolute perturbation that + * is added to the diagonal before + * forming the prec, and + * ilu_rtol is a scaling + * factor $rtol \geq 1$. The last + * parameter specifies the overlap of + * the partitions when the + * preconditioner runs in parallel. + */ + struct AdditionalData + { + /** + * Constructor. By default, no + * element will be dropped, the level + * of extra fill-ins is set to be + * zero (just use the matrix + * structure, do not generate any + * additional fill-in except the one + * that results from non-dropping + * large elements), the tolerance + * level are 0 and 1, respectively, + * and the overlap in case of a + * parallel execution is zero. This + * overlap in a block-application of + * the ILU in the parallel case makes + * the preconditioner a so-called + * additive Schwarz preconditioner. + */ + AdditionalData (const double ilut_drop = 0., + const unsigned int ilut_fill = 0, + const double ilut_atol = 0., + const double ilut_rtol = 1., + const unsigned int overlap = 0); + + /** + * This specifies the relative size + * of elements which should be + * dropped when forming an incomplete + * LU decomposition with threshold. + */ + double ilut_drop; + + /** + * This specifies the amount of + * additional fill-in elements + * besides the sparse matrix + * structure. When + * ilu_fill is large, + * this means that many + * fill-ins will be added, so + * that the ILU preconditioner + * comes closer to a (direct) + * sparse LU + * decomposition. Note, + * however, that this will + * drastically increase the + * memory requirement, + * especially when the + * preconditioner is used in + * 3D. + */ + unsigned int ilut_fill; + + /** + * This specifies the amount of + * an absolute perturbation + * that will be added to the + * diagonal of the matrix, + * which sometimes can help to + * get better preconditioners. + */ + double ilut_atol; + + /** + * This specifies the factor by + * which the diagonal of the + * matrix will be scaled, which + * sometimes can help to get + * better preconditioners. + */ + double ilut_rtol; + + /** + * This determines how large + * the overlap of the local + * matrix portions on each + * processor in a parallel + * application should be. + */ + unsigned int overlap; + }; + + /** + * Initialize function. Takes + * the matrix which is used to + * form the preconditioner, and + * additional flags if there + * are any. + */ + void initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data = AdditionalData()); + + private: + /** + * This is a pointer to the + * Ifpack data contained in + * this preconditioner. + */ + Teuchos::RCP ifpack; + }; + + + /** * A wrapper class for a sparse direct LU decomposition on parallel * blocks for Trilinos matrices. When run in serial, this corresponds diff --git a/deal.II/lac/include/lac/trilinos_solver.h b/deal.II/lac/include/lac/trilinos_solver.h index f7a01c17f1..0ce1af8b7d 100644 --- a/deal.II/lac/include/lac/trilinos_solver.h +++ b/deal.II/lac/include/lac/trilinos_solver.h @@ -18,16 +18,17 @@ #include #include +#include #ifdef DEAL_II_USE_TRILINOS #include #include #include +#include DEAL_II_NAMESPACE_OPEN - namespace TrilinosWrappers { // forward declarations @@ -54,7 +55,7 @@ namespace TrilinosWrappers * deal.II class SolverSelector. * * @ingroup TrilinosWrappers - * @author Martin Kronbichler, 2008 + * @author Martin Kronbichler, 2008, 2009 */ class SolverBase { @@ -81,18 +82,36 @@ namespace TrilinosWrappers * pipe additional data to the * solver. */ + struct AdditionalData { /** - * 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 - * doing this. + * 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 + * iterations etc.). */ - AdditionalData (const unsigned int gmres_restart_parameter = 30); + const bool output_solver_details; + /** + * Restart parameter for GMRES + * solver. + */ const unsigned int gmres_restart_parameter; }; @@ -136,6 +155,30 @@ namespace TrilinosWrappers const VectorBase &b, 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. + */ + void + solve (const SparseMatrix &A, + dealii::Vector &x, + const dealii::Vector &b, + const PreconditionBase &preconditioner); + /** * Access to object that controls * convergence. @@ -210,9 +253,23 @@ namespace TrilinosWrappers * pipe additional data to the * solver. */ + struct AdditionalData - {}; - + { + /** + * 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 + * iterations etc.). + */ + bool output_solver_details; + }; + /** * Constructor. In contrast to * deal.II's own solvers, there is no @@ -252,7 +309,20 @@ namespace TrilinosWrappers * solver. */ struct AdditionalData - {}; + { + /** + * 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 + * iterations etc.). + */ + bool output_solver_details; + }; /** * Constructor. In contrast to @@ -299,12 +369,20 @@ namespace TrilinosWrappers * 30, i.e. do a restart every 30 * iterations. */ - AdditionalData (const unsigned int restart_parameter = 30); - - /** - * Maximum number of - * tmp vectors. - */ + AdditionalData (const bool output_solver_details = false, + const unsigned int restart_parameter = 30); + + /** + * Enables/disables the output of + * solver details (residual in each + * iterations etc.). + */ + bool output_solver_details; + + /** + * Maximum number of + * tmp vectors. + */ unsigned int restart_parameter; }; @@ -347,7 +425,20 @@ namespace TrilinosWrappers * solver. */ struct AdditionalData - {}; + { + /** + * 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 + * iterations etc.). + */ + bool output_solver_details; + }; /** * Constructor. In contrast to @@ -388,7 +479,20 @@ namespace TrilinosWrappers * solver. */ struct AdditionalData - {}; + { + /** + * 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 + * iterations etc.). + */ + bool output_solver_details; + }; /** * Constructor. In contrast to @@ -411,6 +515,145 @@ namespace TrilinosWrappers const AdditionalData additional_data; }; + + +/** + * An implementation of the Trilinos KLU direct solver (using the Amesos + * package). + * + * @ingroup TrilinosWrappers + * @author Martin Kronbichler, 2009 + */ + class SolverDirect + { + public: + + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + + struct AdditionalData + { + /** + * 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 + * iterations etc.). + */ + bool output_solver_details; + }; + + /** + * Constructor. Takes the + * solver control object and + * creates the solver. + */ + SolverDirect (SolverControl &cn, + const AdditionalData &data = AdditionalData()); + + /** + * Destructor. + */ + 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. + */ + void + solve (const SparseMatrix &A, + VectorBase &x, + 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. + */ + void + solve (const SparseMatrix &A, + dealii::Vector &x, + const dealii::Vector &b); + + /** + * Access to object that controls + * convergence. + */ + SolverControl & control() const; + + /** + * Exception + */ + DeclException1 (ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occured while calling a Trilinos function"); + + private: + + /** + * 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 + * Trilinos solver. + */ + std::auto_ptr linear_problem; + + /** + * A structure that contains + * the Trilinos solver and + * preconditioner objects. + */ + std::auto_ptr solver; + + /** + * Store a copy of the flags for this + * particular solver. + */ + const AdditionalData additional_data; + + }; + } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index fe8f684854..de6f2bbd2f 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -2316,7 +2316,7 @@ namespace TrilinosWrappers Assert (values.m() == values.n(), ExcNotQuadratic()); for (unsigned int i=0; icomponents.resize(this->n_blocks()); for (unsigned int i=0;in_blocks();++i) - block(i).clear(); + components[i].clear(); collect_sizes(); } diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index cdf1e1e914..5df4665db0 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -371,6 +371,60 @@ namespace TrilinosWrappers +/* -------------------------- PreconditionILUT -------------------------- */ + + PreconditionILUT::AdditionalData:: + AdditionalData (const double ilut_drop, + const unsigned int ilut_fill, + const double ilut_atol, + const double ilut_rtol, + const unsigned int overlap) + : + ilut_drop (ilut_drop), + ilut_fill (ilut_fill), + ilut_atol (ilut_atol), + ilut_rtol (ilut_rtol), + overlap (overlap) + {} + + + + void + PreconditionILUT::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + preconditioner.release(); + ifpack.release(); + + ifpack = Teuchos::rcp (Ifpack().Create ("ILUT", &*matrix.matrix, + additional_data.overlap)); + + Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("fact: drop value",additional_data.ilut_drop); + parameter_list.set ("fact: level-of-fill",(int)additional_data.ilut_fill); + parameter_list.set ("fact: absolute threshold",additional_data.ilut_atol); + parameter_list.set ("fact: relative threshold",additional_data.ilut_rtol); + parameter_list.set ("schwarz: combine mode", "Add"); + + ierr = ifpack->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = ifpack->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = ifpack->Compute(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + preconditioner = Teuchos::rcp (&*ifpack, false); + } + + + /* ---------------------- PreconditionBlockDirect --------------------- */ PreconditionBlockwiseDirect::AdditionalData:: diff --git a/deal.II/lac/source/trilinos_solver.cc b/deal.II/lac/source/trilinos_solver.cc index 4ca0f6ed8c..d22462b4b7 100644 --- a/deal.II/lac/source/trilinos_solver.cc +++ b/deal.II/lac/source/trilinos_solver.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- +#include #include #include #include @@ -21,27 +22,27 @@ #ifdef DEAL_II_USE_TRILINOS -#include - DEAL_II_NAMESPACE_OPEN namespace TrilinosWrappers { - SolverBase::AdditionalData::AdditionalData (const unsigned int gmres_restart_parameter) + SolverBase::AdditionalData::AdditionalData (const bool output_solver_details, + const unsigned int gmres_restart_parameter) : + output_solver_details (output_solver_details), gmres_restart_parameter (gmres_restart_parameter) {} - - + + SolverBase::SolverBase (SolverControl &cn) : solver_name (gmres), solver_control (cn) {} - + SolverBase::SolverBase (const enum SolverBase::SolverName solver_name, SolverControl &cn) @@ -50,12 +51,12 @@ namespace TrilinosWrappers solver_control (cn) {} - + SolverBase::~SolverBase () {} - + SolverControl & SolverBase::control() const @@ -119,7 +120,101 @@ namespace TrilinosWrappers AssertThrow (ierr == 0, ExcTrilinosError(ierr)); // ... set some options, ... - solver.SetAztecOption (AZ_output, AZ_none); + solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? + AZ_all : AZ_none); + solver.SetAztecOption (AZ_conv, AZ_noscaled); + + // ... and then solve! + ierr = solver.Iterate (solver_control.max_steps(), + solver_control.tolerance()); + 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 (solver.NumIters(), solver.TrueResidual()); + + if (solver_control.last_check() != SolverControl::success) + throw SolverControl::NoConvergence (solver_control.last_step(), + solver_control.last_value()); + } + + + + void + SolverBase::solve (const SparseMatrix &A, + dealii::Vector &x, + const dealii::Vector &b, + const PreconditionBase &preconditioner) + { + int ierr; + + linear_problem.reset(); + + // 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(), + ExcDimensionMismatch(b.size(), A.m())); + Assert (A.local_range ().second == A.m(), + ExcMessage ("Can only work in serial when using deal.II vectors.")); + Assert (A.matrix->Filled(), + ExcMessage ("Matrix is not compressed. Call compress() method.")); + + Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin()); + Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast(b.begin())); + + // We need an + // Epetra_LinearProblem object + // to let the AztecOO solver + // know about the matrix and + // vectors. + linear_problem = std::auto_ptr + (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b)); + + // Next we can allocate the + // AztecOO solver... + solver.SetProblem(*linear_problem); + + // ... and we can specify the + // solver to be used. + switch (solver_name) + { + case cg: + solver.SetAztecOption(AZ_solver, AZ_cg); + break; + case cgs: + solver.SetAztecOption(AZ_solver, AZ_cgs); + break; + case gmres: + solver.SetAztecOption(AZ_solver, AZ_gmres); + solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter); + break; + case bicgstab: + solver.SetAztecOption(AZ_solver, AZ_bicgstab); + break; + case tfqmr: + solver.SetAztecOption(AZ_solver, AZ_tfqmr); + break; + default: + Assert (false, ExcNotImplemented()); + } + + // Introduce the + // preconditioner, ... + ierr = solver.SetPrecOperator (const_cast + (&*preconditioner.preconditioner)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + // ... set some options, ... + solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? + AZ_all : AZ_none); solver.SetAztecOption (AZ_conv, AZ_noscaled); // ... and then solve! @@ -147,11 +242,19 @@ namespace TrilinosWrappers /* ---------------------- SolverCG ------------------------ */ + SolverCG::AdditionalData:: + AdditionalData (const bool output_solver_details) + : + output_solver_details (output_solver_details) + {} + + + SolverCG::SolverCG (SolverControl &cn, const AdditionalData &data) : SolverBase (cn), - additional_data (data) + additional_data (data.output_solver_details) { solver_name = cg; } @@ -160,8 +263,10 @@ namespace TrilinosWrappers /* ---------------------- SolverGMRES ------------------------ */ SolverGMRES::AdditionalData:: - AdditionalData (const unsigned int restart_parameter) + AdditionalData (const bool output_solver_details, + const unsigned int restart_parameter) : + output_solver_details (output_solver_details), restart_parameter (restart_parameter) {} @@ -171,7 +276,8 @@ namespace TrilinosWrappers const AdditionalData &data) : SolverBase (cn), - additional_data (data.restart_parameter) + additional_data (data.output_solver_details, + data.restart_parameter) { solver_name = gmres; } @@ -179,11 +285,20 @@ namespace TrilinosWrappers /* ---------------------- SolverBicgstab ------------------------ */ + SolverBicgstab::AdditionalData:: + AdditionalData (const bool output_solver_details) + : + output_solver_details (output_solver_details) + {} + + + + SolverBicgstab::SolverBicgstab (SolverControl &cn, const AdditionalData &data) : SolverBase (cn), - additional_data (data) + additional_data (data.output_solver_details) { solver_name = bicgstab; } @@ -191,11 +306,20 @@ namespace TrilinosWrappers /* ---------------------- SolverCGS ------------------------ */ + SolverCGS::AdditionalData:: + AdditionalData (const bool output_solver_details) + : + output_solver_details (output_solver_details) + {} + + + + SolverCGS::SolverCGS (SolverControl &cn, const AdditionalData &data) : SolverBase (cn), - additional_data (data) + additional_data (data.output_solver_details) { solver_name = cgs; } @@ -203,15 +327,177 @@ namespace TrilinosWrappers /* ---------------------- SolverTFQMR ------------------------ */ + SolverTFQMR::AdditionalData:: + AdditionalData (const bool output_solver_details) + : + output_solver_details (output_solver_details) + {} + + + SolverTFQMR::SolverTFQMR (SolverControl &cn, const AdditionalData &data) : SolverBase (cn), - additional_data (data) + additional_data (data.output_solver_details) { solver_name = tfqmr; } + + +/* ---------------------- SolverDirect ------------------------ */ + + SolverDirect::AdditionalData:: + AdditionalData (const bool output_solver_details) + : + output_solver_details (output_solver_details) + {} + + + + + SolverDirect::SolverDirect (SolverControl &cn, + const AdditionalData &data) + : + solver_control (cn), + additional_data (data.output_solver_details) + {} + + + + SolverDirect::~SolverDirect () + {} + + + + SolverControl & + SolverDirect::control() const + { + return solver_control; + } + + + + void + SolverDirect::solve (const SparseMatrix &A, + VectorBase &x, + const VectorBase &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(); + + // We need an + // Epetra_LinearProblem object + // to let the AztecOO solver + // know about the matrix and + // vectors. + linear_problem = std::auto_ptr ( + new Epetra_LinearProblem(&*(A.matrix), &*x.vector, + &*b.vector)); + + // Next we can allocate the + // AztecOO solver... + solver = std::auto_ptr (Amesos().Create("Amesos_Klu", + *linear_problem)); + + 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) + throw SolverControl::NoConvergence (solver_control.last_step(), + solver_control.last_value()); + } + + + + 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. + Assert (x.size() == A.n(), + ExcDimensionMismatch(x.size(), A.n())); + Assert (b.size() == A.m(), + ExcDimensionMismatch(b.size(), A.m())); + Assert (A.local_range ().second == A.m(), + ExcMessage ("Can only work in serial when using deal.II vectors.")); + Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin()); + Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast(b.begin())); + + // We need an + // Epetra_LinearProblem object + // to let the AztecOO solver + // know about the matrix and + // vectors. + linear_problem = std::auto_ptr + (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b)); + + // Next we can allocate the + // AztecOO solver... + solver = std::auto_ptr (Amesos().Create("Amesos_Klu", + *linear_problem)); + + 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) + throw SolverControl::NoConvergence (solver_control.last_step(), + solver_control.last_value()); + } + + + + } DEAL_II_NAMESPACE_CLOSE -- 2.39.5