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.

Reply via email to