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.
