Dear Timo, Thank you for your reply.
I am still having trouble with implementing my code with MUMPS. I will briefly describe the problem: I am solving 2 systems of Maxwell's equations in the following way: the systems are Ax1=b1 and Ax2=b2; A is a sparse block symmetric matrix (C -M; -M -C). I solve it with FGMRES preconditioned with Pa= (B 0; 0 B) = (C+M 0; 0 C+M). To obtain inv(Pa) I need to obtain inv(B). And this is what I would like to achieve with MUMPS. However, to save time I want to compute inv(B) only ones for both b1 and b2. So I would like to reuse LDLt the second time I solve the system with b2. With SparseDirectUMFPACK<double> I did as following. 1) I declared std_cxx11::shared_ptr<typename LinearSolvers::InnerPreconditioner<inner_solver_id>::type> matrix_B_preconditioner; in my main class, 2) at the beginning of setup_dofs I had matrix_B_preconditioner.reset(); 3) I assembled the system and rhs separately, and after assembling the system matrix I defined matrix_B_preconditioner = std_cxx11::shared_ptr< typename LinearSolvers::InnerPreconditioner<inner_solver_id>::type > ( new typename LinearSolvers::InnerPreconditioner<inner_solver_id>::type() ); matrix_B_preconditioner->initialize(matrix_B.block(0, 0), typename LinearSolvers::InnerPreconditioner<inner_solver_id>::type::AdditionalData() ); 4) when solving const LinearSolvers::Pa_Preconditioner<typename LinearSolvers::InnerPreconditioner<inner_solver_id>::type > pa_preconditioner(*matrix_B_preconditioner); Where a) InnerPreconditioner was template<int inner_solver_id> struct InnerPreconditioner; template<> struct InnerPreconditioner<0> { typedef SparseDirectUMFPACK<double> type; //typedef PETScWrappers::SparseDirectMUMPS type; }; b) Pa_Preconditioner was: template<class PreconditionerB> class Pa_Preconditioner: public Subscriptor { public: Pa_Preconditioner(const PreconditionerB &preconditioner_B); void vmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; private: const PreconditionerB &preconditioner_B; }; template<class PreconditionerB> Pa_Preconditioner<PreconditionerB>:: Pa_Preconditioner(const PreconditionerB &preconditioner_B) : preconditioner_B (&preconditioner_B) { } template<class PreconditionerB> void Pa_Preconditioner<PreconditionerB>:: vmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { preconditioner_B.vmult(dst.block(0), src.block(0)); preconditioner_B.vmult(dst.block(1), src.block(1)); } As far as I understand with MUMPS I have to change this strategy and modify Pa_Preconditioner class and have matrix B as input to it, and change Pa_Preconditioner::vmult to template<class Btype> Pa_Preconditioner<Btype>:: Pa_Preconditioner(const Btype &B,const MPI_Comm &mpi_communicator) : B (&B), mpi_communicator (&mpi_communicator) { } template<class Btype> void Pa_Preconditioner<Btype>:: vmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { SolverControl cn; PetScWrappers::SparseDirectMUMPS solver(cn,mpi_communicator); solver.set_symmetric_mode(true); solver.solve(B,dst.block(0),src.block(0)); solver.solve(B,dst.block(1),src.block(1)); } Does this look reasonable? I think, this makes program much less general than my previous version where I could choose between direct and iterative solvers with the template parameter. And also I am not sure that LDLt decomposition of B is reused in this case. Any suggestions will be highly appreciated. Thank you in advance, Anna Show trimmed content -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.