Enhancements on IterativeSolver, ComplexPetscParMatrix, PetscParVector

[ghost]

Dear developers,

Here I want to propose a few enhancements:

For IterativeSolver:
(1) member function to judge if the IterativeSolver has converged to the preset tolerance
(2) member function which returns the absolute (or absolute) tolerance for each linear solve
These two functions will be very useful when the users don't want the solver to print any information to the standard output stream.

(3) a new data structure ComplexPetscParMatrix
Currently only ComplexHypreParMatrix is implemented in the complex version of mfem. There is no member function SetOperatorType(Operator::PETSC_MATAIJ) for the class ParSesquilinearForm either, so we can't obtain a PETSc format matrix after calling FormLinearSystem. This new matrix format will introduce more solvers from PETSc (like MUMPS) for complex mfem.

(4) more flexible PetscParVector
I noticed the following code:
/// Creates a PetscParVector from a Vector (data is not copied) PetscParVector(MPI_Comm comm, const Vector &_x);
When I want to get a PetscParVector from a Vector (or vice versa), I have to copy the data manually. Maybe more functions (like operator= for copying Vector and PetscParVector) are expected.

Thanks for your work and consideration.

Best,
Changkai

[tzanio]

Thanks for the suggestions, Changkai!

(1) member function to judge if the IterativeSolver has converged to the preset tolerance

Can you clarify why IterativeSolver:: GetConverged() is not sufficient for your needs?

(2) member function which returns the absolute (or absolute) tolerance for each linear solve

Are you thinking here of NewtonSolver? Can your use case be addressed if we add a methods IterativeSolver::GetRelTol() and IterativeSolver::GetAbsTol()?

(3) a new data structure ComplexPetscParMatrix
(4) more flexible PetscParVector

@stefanozampini is probably the best person to address these.

[ghost]

Hi Tzanio,
Thanks for your attention.

(1) IterativeSolver::GetConverged() has public attribute in the base class IterativeSolver, while its's not visible for CGSolver and GMRESSolver. (Hope I understand the code correctly. I am not so proficient in C++). This may be not so important if we have IterativeSolver::GetRelTol() or IterativeSolver::GetAbsTol().

(2) The second point is also for IterativeSolver. I mean for each linear system solve (but not for NewtonSolver ). IterativeSolver::GetRelTol() and/or IterativeSolver::GetAbsTol() would be very useful. By the way, the HyprePCG solver prints the relative tolerance (||r||_C/||b||_C), while the fgmres solver prints the absolute tolerance (||r||). The inconsistent behaviors may make the users confused.

Changkai

[ghost]

Update:
GetConverged() is visible for all IterativeSolver. I am sorry. I just mix CGSolver : public IterativeSolver with HyprePCG.

It seems that there are no functions like GetAbsTol or GetRelTol for HyprePCG either.

[tzanio]

Hi Changkai,

Are you interested in proposing a pull request to add the GetRel/AbsTol() functions you need?

Tzanio

[ghost]

Hi Tzanio, I am glad to make any contributions (if I could). It may take some time because I have no experience with git/github and I have to read the source code solvers and hypre.

[stefanozampini]

@qiuchangkai

(3) a new data structure ComplexPetscParMatrix
Currently only ComplexHypreParMatrix is implemented in the complex version of mfem. There is no member function SetOperatorType(Operator::PETSC_MATAIJ) for the class ParSesquilinearForm either, so we can't obtain a PETSc format matrix after calling FormLinearSystem. This new matrix format will introduce more solvers from PETSc (like MUMPS) for complex mfem.

Where is this complex version of MFEM? AFAIK MFEM fakes complex operators using special real block operators. You can always extract the blocks that you want and convert them into a PetscParMatrix. Can you be more specific about your needs? PETSc can be configured with support for complex numbers, i.e. MATAIJ with complex entries. Is this what you would like to have to use MUMPS?

4. more flexible PetscParVector
   I noticed the following code:
   /// Creates a PetscParVector from a Vector (data is not copied) PetscParVector(MPI_Comm comm, const Vector &_x);
   When I want to get a PetscParVector from a Vector (or vice versa), I have to copy the data manually. Maybe more functions (like operator= for copying Vector and PetscParVector) are expected.

There are some enhancements here https://github.com/mfem/mfem/tree/stefanozampini/small-additions . I have plans to open a PR. Tell me if they fit your needs.

[ghost]

Hi @stefanozampini ,

Thanks for your attention. Please accept my apology for the ambiguous description.

I know the support of complex operators are handled by block real operators in mfem-complex-mfem-dev.

For real-value problems, OperatorHandle provides a high level common interface for operators, so users can pass a OperatorHandle to ParBilinearForm->AssembleSystemMatrix(ess_tdof_list, OperatorHandle);, and then transform the class type from OperatorHandle to expected matrix type for HyreSolver, MFEM solver, Petsc solver via dynamic_cast.

Below are my needs:

While, for complex-value problems (in the real block forms), only ComplexSparseMatrix and ComplexHypreParMatrix are implemented in ComplexOperator.hpp. If I pass a OperatorHandle to
ParSesquilinearForm->FormLinearSystem(ess_tdof_list, OperatorHandle), this OperatorHandle can only be used for Hypre solvers, IterativeSolver, but it doesn't work for Petsc KSP.

The following code
ierr = MatSetOption(dynamic_cast<PetscParMatrix&>(*OperatorHandle.Ptr()), MAT_SYMMETRIC , PETSC_TRUE);
will report
terminate called after throwing an instance of 'std::bad_cast' what(): std::bad_cast

So I expect a data structure which stores the sparse matrix arising from complex block problems in the Petsc Mat format, like that implemented for the real-value problem.

Hope I have make my needs clear. Thank you very much!

Best wishes,
Changkai

[stefanozampini]

the following code
ierr = MatSetOption(dynamic_cast<PetscParMatrix&>(*OperatorHandle.Ptr()), MAT_SYMMETRIC , PETSC_TRUE);
will report
terminate called after throwing an instance of 'std::bad_cast' what(): std::bad_cast

There should be no need for you to explicitly call PETSc functions.
If you need some functionality in PETSc which is not exposed through MFEM, let me know.

Once the complex branch will get merged to master I'll see if it is reasonable to add the ComplexPetscParMatrix class.

For the time being, if you just need the action of the matrix in the Krylov solver, you don't need to pass a PetscParMatrix to a PetscLinearSolver; you can use any mfem::Operator with the SetOperator() method and pass the ComplexHypreParMatrix object ( which I presume is derived from the mfem::Operator class). The code will then use the operator's Mult() method to solve the linear system. You can also create your own Hypre based preconditioner (deriving from Solver) and call the SetPreconditioner method of the PetscLinearSolver class to use it (again, through the Mult() method) inside a KSP.

[ghost]

Hi @stefanozampini ,

Thanks for your kind suggestions about the PetscLinearSolver. I tried to use MUMPS through PetscLinearSolver today with the following code:

     PetscLinearSolver petscsolver(MPI_COMM_WORLD);

     KSP ksp;
     PC pc;
     Mat F;
     PetscBool flg_mumps_lu = PETSC_TRUE;
     PetscBool flg_mumps_ldlt = PETSC_FALSE;

     petscsolver.SetOperator(*A.Ptr(), *A.Ptr());

     ksp = (KSP)(petscsolver);

     PetscErrorCode ierr;

     ierr = KSPSetType(ksp, KSPPREONLY);

     ierr = KSPGetPC(ksp, &pc);

     if (flg_mumps_lu)
     {
        ierr = PCSetType(pc, PCLU);
     }
     else if (flg_mumps_ldlt)
     {
        ierr = PCSetType(pc, PCCHOLESKY);
     }

     ierr = PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
     ierr = PCFactorSetUpMatSolverType(pc);
     ierr = PCFactorGetMatrix(pc, &F);

     B = 1.0e0;

     petscsolver.Mult(B,U);

     for (int i = 0; i < U.Size(); i++)
     {
        cout << U[i] << endl;
     }

There is no problem when I tested with mfem/examples/ex3p.cpp (a ParBilinearForm), but it fails with mfem/examples/ex21p.cpp (a ParSesquilinearForm, it's from the latest complex-mfem-dev). I just simply add the above code after a->FormLinearSystem(ess_tdof_list, u, b, A, U, B). The errors are:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for possible LU and Cholesky solvers
[0]PETSC ERROR: MatSolverType mumps does not support matrix type shell
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.9.4, Sep, 11, 2018
[0]PETSC ERROR: Unknown Name on a arch-intel-opt named debian by gem Tue Nov 6 03:42:52 2018
[0]PETSC ERROR: Configure options --with-cc=mpiicc --with-cxx=mpiicpc --with-fc=mpiifort --with-blaslapack-dir=/opt/intel/mkl --with-scalapack-lib="-L/opt/intel/mkl/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64" --download-metis=yes --download-mumps=yes --download-hypre=yes
[0]PETSC ERROR: #1 MatGetFactor() line 4332 in /home/gem/Desktop/GEM/petsc-3.9.4/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 PCFactorSetUpMatSolverType_Factor() line 15 in /home/gem/Desktop/GEM/petsc-3.9.4/src/ksp/pc/impls/factor/factimpl.c
[0]PETSC ERROR: #3 PCFactorSetUpMatSolverType() line 60 in /home/gem/Desktop/GEM/petsc-3.9.4/src/ksp/pc/impls/factor/factor.c

Am I missing something in the above code?

Best wishes,
Changkai

[stefanozampini]

As I said before, since the general convertor in MFEM does not know anything about the complex class of matrices, it creates a MATSHELL which just calls the Mult method of the class when doing matvecs.
MUMPS requires knowledge of the matrix entries. To get the entries of the matrix, you should compute the operator explicitly, which requries N matvecs, with N the size of the matrix, via PetscParMatrix *convMat = new PetscParMatrix(comm,A.Ptr(),Operator::PETSC_MATAIJ);

[ghost]

Well, I don't understand 'which requries N matvecs, with N the size of the matrix'. Thanks.

[stefanozampini]

Being a shell matrix, you need N matrix vector products to compute all the entries (required by MUMPS)

[ghost]

Hi @stefanozampini ,

As you suggested, I first call

PetscParMatrix *convMat = new PetscParMatrix(MPI_COMM_WORLD,A.Ptr(),Operator::PETSC_MATAIJ);

to obtain the matrix in PetscParMatrix format.

Then, either

 petscsolver.SetOperator(*convMat, *convMat);

or

ksp = (KSP)(petscsolver);

KSPSetOperators(ksp,Mat(*convMat),Mat(*convMat));

won't work. Could you show me the pseudo code how I can solve the complex system using PETSc MUMPS? I have no experience with PETSc before, but use the Fortran version MUMPS directly.

Thanks a lot.

Changkai

[ghost]

Hi @stefanozampini ,

Could you please help me about this issue? I want to use MUMPS through PETSc to verify my solutions.

I am still working on a problem like ex21p with Dirichlet boundary condition E=E0, where E0 is the electric field on the boundary. I set up the coefficient matrix (and rhs vector) correctly.

Currently, the residual by FGMRESSolver will drop below the preset tolerance after the third iteration if I set the initial guess to be zero, and the output is like:

Pass : 1 Iteration : 0 || r || = 367.143
Pass : 1 Iteration : 1 || r || = 259.609
Pass : 1 Iteration : 2 || r || = 4.96954e-13

In this case, the final solution vector are all zero. If the initial guess is a small number like 1.0e-8, the FGMRESSolver will converge normally after several iterations, but I don't know whether the solution is correct or not (at least it's not what I expct).

@mlstowell , any suggestions about what I should do? How can I solve the system equation in ex21p with MUMPS solver? FGMRESSolver fails to solve the equation, so I want to find an alternative solver like HypreFlexGMRES from Hypre, but the BlockDiagonalPreconditioner in ex21p is not supported by HypreSolver.

Many thanks for all your help.

Changkai