PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[shanglong zhang]

Dear all,

I modify the step-17 so that it uses PETScWrappers::SparseDirectMUMPS as solver. However, when I active the option of set_symmetric_mode, the solution has discrepancy from CG or MUMPS without the symmetry option. (The solutions from CG and MUMPS without symmetry option are identical though.) Any suggestion?

And it seems to me that  PETScWrappers::SparseDirectMUMPS doesn't have the functionality of reusing the factorization of system matrix for different right-hand side vector. Do I miss anything?

Thank you very much.

Best,

Shanglong

[Denis Davydov]

Hi,

I can't comment on the first part of your question, but as for the second part -- you are right, 

currently you can't reuse factorization. I you want it, it should be relatively easy to add.

See Direct solver in TrilinosWrappers which was recently modified to support this: 

https://www.dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1SolverDirect.html#a6c0fecc921bc93f2c005e265cfbdd882

Regards,

Denis.

[Alexander]

I modify the step-17 so that it uses PETScWrappers::SparseDirectMUMPS as solver. However, when I active the option of set_symmetric_mode, the solution has discrepancy from CG or MUMPS without the symmetry option. (The solutions from CG and MUMPS without symmetry option are identical though.) Any suggestion?

Which version of PETSc, MUMPS and deal.II do you use? What is the matrix size you work with?
 

And it seems to me that  PETScWrappers::SparseDirectMUMPS doesn't have the functionality of reusing the factorization of system matrix for different right-hand side vector. Do I miss anything?

Yes, it does. The first time you call solve it performs factorization, but the following calls will reuse it (since you don't call KSPSetOperators).
 

[shanglong zhang]

Hi Denis,

Thank you for your information.

Best,

Shanglong

在 2016年7月25日星期一 UTC-4下午12:02:09，Denis Davydov写道：

[shanglong zhang]

Hi Alexander,

The versions of library I am using are listed below:

deal.II 8.4.0

petsc 3.3.0.3

MUMPS 4.10.0-p3

I run the first two mesh refining loops of the step-17. So the matrix size is 162*162 and 302*302 respectively. The symmetry direct solver doesn't give the 'right solution' for both cases.

Sorry, I am not familiar with the petsc library. You meant if I call solver.solve(system_matrix,solution,system_rhs) twice, internally, it would use the factorization from the first time for the second solving process? (of course the system matrix is the same.) Do i understand this correctly? Thanks.

Best,

Shanglong

在 2016年7月25日星期一 UTC-4下午3:24:25，Alexander写道：

[Denis Davydov]

Inside SparseDirectMUMPS::solve  there is https://www.dealii.org/developer/doxygen/deal.II/petsc__solver_8cc_source.html#l00703

 if (solver_data.get() == 0)
  {

 solver_data.reset (new SolverDataMUMPS ());

…

ierr = KSPSetOperators (solver_data->ksp, A, A);

…

  }

ierr = KSPSolve (solver_data->ksp, b, x);

so i think Alexander is right that it will reuse factorisation.

However, this is wrong because you might be calling the same solver consequently

but with different matrices and RHS vectors. In other words, the second time you call solve() 

matrix is ignored altogether.

So i would say that it would still be good to refactor this solver similar to the director solver in Trilinos 

to make it (a) more user friendly (b) have less surprises when one calls solve() for different matrices consequently.

Kind regards,

Denis 

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/Sw0wHKEFuiA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Alexander]

deal.II 8.4.0

petsc 3.3.0.3

MUMPS 4.10.0-p3

This PETSc and MUMPS > 4 years old. Update to the latest versions and try again. If problem persists, provide the relevant snippet of the code.
 

Sorry, I am not familiar with the petsc library. You meant if I call solver.solve(system_matrix,solution,system_rhs) twice, internally, it would use the factorization from the first time for the second solving process? (of course the system matrix is the same.) Do i understand this correctly? Thanks.

Yes, correct.

@Denis: I agree, the logic is somewhat deficient. I also think documentation and interface itself (e.g. I have no idea why does method set_solver_type exist and exposed to the public) could be improved.

[Denis Davydov]

Hi Alexander,

On 26 Jul 2016, at 09:40, Alexander <agra...@gmail.com> wrote:

@Denis: I agree, the logic is somewhat deficient. I also think documentation and interface itself (e.g. I have no idea why does method set_solver_type exist and exposed to the public) could be improved.

I created a GitHub issue to track this https://github.com/dealii/dealii/issues/2869

As per set_solver_type, it is protected so it is not exposed to users. Should be ok, AFAICT.

Regards,

Denis.

[shanglong zhang]

Hi Alexander,

I switch to :

dealii 8.4.1

petsc 3.6.4

mumps 5.0.1

The result is still not correct.

The modification I did for the code is replacing the following code in step-17.cc:

    SolverControl           solver_control (solution.size(),

                                            1e-8*system_rhs.l2_norm());

    PETScWrappers::SolverCG cg (solver_control,

                                mpi_communicator);

    PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);

    cg.solve (system_matrix, solution, system_rhs,

              preconditioner);

with

 SolverControl solver_control (solution.size(), 1e-12);

 PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator);

 solver.set_symmetric_mode(true);

 solver.solve(system_matrix, solution, system_rhs);

Thank you very much.

Regards,

Shanglong

在 2016年7月26日星期二 UTC-4上午3:40:09，Alexander写道：

[Alexander]

Shanglong,
I can reproduce this behaviour, but this isn't a bug. I output the 162x162 system matrix at first cycle and as far as I can see it fails symmetry check. You can see this already from first few non-zero entries:

(0,0) 1.33333
(1,1) 1.33333
(2,2) 1.33333
(3,3) 1.33333
(4,4) 1.33333
(5,5) 1.33333
(6,0) -0.666667

See, there is no (0, 6) and symmetry is violated.

The matrix is still positive-definite, that is probably why CG still gets something out of it.

Alexander

[shanglong zhang]

Hi Alexander,

Thank you for your explanation. I should have checked the symmetricity first.

However, since the step-17 is isotropic linear elasticity problem, should the system matrix be symmetry?

Thank you very much for your help and patience.

Regards,

Shanglong

在 2016年7月27日星期三 UTC-4上午4:49:53，Alexander写道：

[Alexander]

Shanglong,
I'm not an expert in this field, but I assume it should be.

Can you please check if system matrix in step-8 (which was used to implement step-17) is symmetric?

Alexander

[shanglong zhang]

Alexander,

I run step-8 and the system matrix is symmetry at least for the first cycle. I attach the output txt of the system matrix in this reply. Thanks.

Best,
Shanglong

在 2016年7月27日星期三 UTC-4上午10:53:46，Alexander写道：

[Alexander]

So why is it not 162x162?

step-17 claims to solve same problem as step-8. If one produces symmetry matrix and the other does not or both do not -- this is suspicious.
Anyway, I recommend you creating an issue on github and put short description with the link to this thread. People who created step-17 will more likely see it there than go into details of this thread.

Best,
Alexander

[Wolfgang Bangerth]

On 07/27/2016 01:43 PM, Alexander wrote:
>
> step-17 claims to solve same problem as step-8. If one produces symmetry
> matrix and the other does not or both do not -- this is suspicious.

The reason why the matrix in step-17 is not symmetric is actually
explained in the code:

// The last argument to the call to
// MatrixTools::apply_boundary_values() below allows for some
// optimizations. It controls whether we should also delete
// entries (i.e., set them to zero) in the matrix columns
// corresponding to boundary nodes, or to keep them (and passing
// <code>true</code> means: yes, do eliminate the columns). If we
// do eliminate columns, then the resulting matrix will be
// symmetric again if it was before; if we don't, then it
// won't. The solution of the resulting system should be the same,
// though. The only reason why we may want to make the system
// symmetric again is that we would like to use the CG method,
// which only works with symmetric matrices. The reason why we may
// <i>not</i> want to make the matrix symmetric is because this
// would require us to write into column entries that actually
// reside on other processes, i.e., it involves communicating
// data. This is always expensive.
//
// Experience tells us that CG also works (and works almost as
// well) if we don't remove the columns associated with boundary
// nodes, which can be explained by the special structure of this
// particular non-symmetry. To avoid the expense of communication,
// we therefore do not eliminate the entries in the affected
// columns.
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(dim),
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs,
false);

In other words, it is expected that the matrix is not symmetric, and
that consequently when you run a direct solver, you can't use the
symmetric setting.

Best
Wolfgang

[shanglong zhang]

Hi Wolfgang,

Thank you for the explanation. Then everything make sense.

Best,

Shanglong

在 2016年7月27日星期三 UTC-4下午6:13:16，bangerth写道：