Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean up linear solvers, add StencilDiagonalMatrix and diagonal(), require pc to be LinearOperator #360

Merged
merged 37 commits into from
Jan 31, 2024

Conversation

spossann
Copy link
Collaborator

@spossann spossann commented Dec 6, 2023

Reduce code duplication in the linear solvers by adding the method __init__ to the common base class InverseLinearOperator. Add new class StencilDiagonalMatrix to be used in combination with the new method diagonal of StencilMatrix and BlockLinearOperator -- this fixes #369. Require all preconditioners to be LinearOperator objects -- this fixes #343.

Detailed list of changes

All solvers are sub-classes of InverseLinearOperator. As such, they now call now super().__init__ which takes the base arguments for any iterative solver. The method _check_options is moved to the base class and checks the common parameters x0, tol, maxiter and verbose. Additional arguments must be checked in the individual solvers.

The methods get_options, set_options and transpose are also moved to the base class. In particular get_options has been improved as it accepts now a specific option name. The property options has been removed because unsafe (it gave direct access to the underlying dictionary without performing any check). The method set_options must be used instead.

The base class now has a setter method for linop, which allows to change the matrix A of the solver.

Moreover, changed the handling of preconditioners in InverseLinearOperators. The option pc must now take as value a LinearOperator, which is the preconditioner. If pc=None, an IdentityOperator of the domain space is created automatically.

The new class StencilDiagonalMatrix has been added. This represents a linear operator acting between two StencilVectorSpaces, and which corresponds to a diagonal matrix. This linear operator is completely local.

The method diagonal has been added to the classes StencilDiagonalMatrix and BlockLinearOperator (for which it returns an object of the same class), and to the class StencilMatrix (for which it returns a StencilDiagonalMatrix object). As well as an out optional argument, this method accepts the inverse boolean flag which is useful for creating a Jacobi preconditioner.

The new solver PBiConjugateGradientStabilized (preconditioned BiCGSTAB) has been added from Struphy (does not work for complex).

test_solvers.py has been updated, as well as several API unit tests.

Notes

When passing a preconditioner to DiscreteEquation.set_solver(), this also has to be a LinearOperator now. This means that in the case of a Jacobi preconditioner one needs first to assemble the linear system, then extract its diagonal, and finally pass the new LinearOperator as the pc solver option.

The option `pc` must now take as value a `LinearOperator`, which is the preconditioner.
If `pc=None`, an `IdentityOperator` of the domain space is created automatically.
The class `InverseLinearOperators` could be deprecated, since the Jacobi preconditiner should be implemented as a `LinearOperator`.

The new solver `PBiConjugateGradientStabilized` (preconditioned) has been added from Struphy.

`test_solvers.py` has been updated. The Jacobi preconditioner is not tested at the moment.
…e base class `InverseLinearOperator`.

All solvers call now `super().__init__` which takes the base arguments for any iterative solver.
Additional arguments must be checked in the individual solvers.

The base class now has the setter method for `linop`, which allows to change tha matrix A of the solver.
Moreover, there is the abstract property `solver` which is the name string of the solver.
@spossann spossann changed the title Moved _check_options() and transpose() to the base class InverseLinearOperator Moved _check_options() and transpose() to InverseLinearOperator, require pc to be LinearOperator Dec 6, 2023
psydac/linalg/basic.py Show resolved Hide resolved
psydac/linalg/basic.py Outdated Show resolved Hide resolved
@yguclu yguclu changed the title Moved _check_options() and transpose() to InverseLinearOperator, require pc to be LinearOperator Clean up linear solvers, add StencilDiagonalMatrix and diagonal(), require pc to be LinearOperator Jan 23, 2024
The text clarifies that replacing the original linear operator through
the setter is a dangerous operation which should only be made in extreme
cases.
@yguclu yguclu requested review from jowezarek, e-moral-sanchez and vcarlier and removed request for jowezarek and e-moral-sanchez January 30, 2024 13:10
psydac/linalg/tests/test_linalg.py Outdated Show resolved Hide resolved
psydac/linalg/tests/test_linalg.py Outdated Show resolved Hide resolved
psydac/linalg/solvers.py Outdated Show resolved Hide resolved
psydac/linalg/block.py Outdated Show resolved Hide resolved
psydac/linalg/basic.py Outdated Show resolved Hide resolved
psydac/linalg/solvers.py Outdated Show resolved Hide resolved
psydac/linalg/basic.py Show resolved Hide resolved
@e-moral-sanchez
Copy link
Contributor

e-moral-sanchez commented Jan 31, 2024

The pc flag can only be passed when the solver is either 'pcg' or 'pbicgstab', otherwise it returns an error.

In the case of solvers like GMRES, we could easily pass a (left) preconditioner P, and then instead of solving A x = b, we solve P @ A x = P.dot(b) with a ComposedLinearOperator.

@yguclu
Copy link
Member

yguclu commented Jan 31, 2024

The pc flag can only be passed when the solver is either 'pcg' or 'pbicgstab', otherwise it returns an error.

In the case of solvers like GMRES, we could easily pass a (left) preconditioner P, and then instead of solving A x = b, we solve P @ A x = P.dot(b) with a ComposedLinearOperator.

AFAIK this was the original behavior prior to this PR. I suggest to defer any further improvements or fixes to a future PR. Please feel free to open new issues.

@yguclu yguclu merged commit e4b5580 into devel Jan 31, 2024
6 checks passed
@yguclu yguclu deleted the 343-pc-as-linop branch January 31, 2024 15:40
yguclu pushed a commit that referenced this pull request Mar 13, 2024
Fix regression bug in file `psydac/examples/maxwell_2d_multi_patch.py`:

Pass a `LinearOperator` object as preconditioner to the PCG solver,
instead of just a string with the method name 'jacobi' (which is not possible
since PR #360). To achieve this we follow the approach in
`psydac/api/tests/test_2d_multipatch_mapping_maxwell.py`:

1. Assemble the linear system explicitly;
2. Extract the diagonal linear operator from the matrix;
3. Pass the linear operator above as preconditioner to the solver.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add class DiagonalStencilMatrix Require a preconditioner to be a LinearOperator
5 participants