Matrix-Free Linear Operator (#423)

Provide a `MatrixFreeLinearOperator` class which is a subclass of the
abstract class `LinearOperator`.

The new class allows creating a general matrix-free linear operator. The
constructor only requires the domain, codomain and a callable `dot`
method.

Notes:
- The provided `dot` method may or may not take an `out` argument.
- A `transpose_dot` method may also be provided (mandatory to
instantiate the `transpose()` linear operator).
- Unit tests have been added.

Additional changes: 
- We stop passing the deprecated `tol` argument in calls to SciPy's
minres and we use `rtol` instead, which requires SciPy >= 1.12. This
fixes #419.
- We stop supporting Python 3.8 because of SciPy 1.12. Python 3.8 is
close to end-of-life anyway, see https://devguide.python.org/versions.
- We also verify the initial convergence in the minres linear solver to
avoid iterating when the initial solution is good enough (which may
cause a division-by-zero error).
- We also revert the changes to files *psydac/api/settings.py* and
*pyproject.toml* which were erroneously made in PR #320

---------

Co-authored-by: e-moral-sanchez <[email protected]>
Co-authored-by: Martin Campos Pinto <[email protected]>

.github/workflows/continuous-integration.yml

@@ -16,11 +16,10 @@ jobs:
       fail-fast: false
       matrix:
         os: [ ubuntu-latest ]
-        python-version: [ 3.8, 3.9, '3.10', '3.11' ]
+        python-version: [ 3.9, '3.10', '3.11' ]
         isMerge:
           - ${{ github.event_name == 'push' && github.ref == 'refs/heads/devel' }}
         exclude:
-          - { isMerge: false, python-version: 3.9 }
           - { isMerge: false, python-version: '3.10' }
         include:
           - os: macos-latest

psydac/api/settings.py

@@ -67,4 +67,4 @@
     'pyccel-intel' : PSYDAC_BACKEND_IPYCCEL,
     'pyccel-pgi'   : PSYDAC_BACKEND_PGPYCCEL,
     'pyccel-nvidia': PSYDAC_BACKEND_NVPYCCEL,
-}
+}

psydac/api/tests/test_api_2d_compatible_spaces.py

@@ -130,7 +130,7 @@ def run_stokes_2d_dir(domain, f, ue, pe, *, homogeneous, ncells, degree, scipy=F
     # ... solve linear system using scipy.sparse.linalg or psydac
     if scipy:
 
-        tol = 1e-11
+        rtol = 1e-11
         equation_h.assemble()
         A0 = equation_h.linear_system.lhs.tosparse()
         b0 = equation_h.linear_system.rhs.toarray()
@@ -145,17 +145,17 @@ def run_stokes_2d_dir(domain, f, ue, pe, *, homogeneous, ncells, degree, scipy=F
             A1 = a1_h.assemble().tosparse()
             b1 = l1_h.assemble().toarray()
 
-            x1, info = sp_minres(A1, b1, tol=tol)
+            x1, info = sp_minres(A1, b1, rtol=rtol)
             print('Boundary solution with scipy.sparse: success = {}'.format(info == 0))
 
-            x0, info = sp_minres(A0, b0 - A0.dot(x1), tol=tol)
+            x0, info = sp_minres(A0, b0 - A0.dot(x1), rtol=rtol)
             print('Interior solution with scipy.sparse: success = {}'.format(info == 0))
 
             # Solution is sum of boundary and interior contributions
             x = x0 + x1
 
         else:
-            x, info = sp_minres(A0, b0, tol=tol)
+            x, info = sp_minres(A0, b0, rtol=rtol)
             print('Solution with scipy.sparse: success = {}'.format(info == 0))
 
         # Convert to stencil format

psydac/linalg/basic.py

@@ -8,6 +8,8 @@
 """
 
 from abc import ABC, abstractmethod
+from types import LambdaType 
+from inspect import signature
 
 import numpy as np
 from scipy.sparse import coo_matrix
@@ -25,7 +27,8 @@
     'ComposedLinearOperator',
     'PowerLinearOperator',
     'InverseLinearOperator',
-    'LinearSolver'
+    'LinearSolver',
+    'MatrixFreeLinearOperator'
 )
 
 #===============================================================================
@@ -1111,3 +1114,108 @@ def solve(self, rhs, out=None):
     @property
     def T(self):
         return self.transpose()
+
+#===============================================================================
+class MatrixFreeLinearOperator(LinearOperator):
+    """
+    General linear operator represented by a callable dot method.
+
+    Parameters
+    ----------
+    domain : VectorSpace
+        The domain of the linear operator.
+    
+    codomain : VectorSpace
+        The codomain of the linear operator.
+    
+    dot : Callable
+        The method of the linear operator, assumed to map from domain to codomain.    
+        This method can take out as an optional argument but this is not mandatory.
+
+    dot_transpose: Callable
+        The method of the transpose of the linear operator, assumed to map from codomain to domain.
+        This method can take out as an optional argument but this is not mandatory.
+
+    Examples
+    --------
+        # example 1: a matrix encapsulated as a (fake) matrix-free linear operator        
+        A_SM = StencilMatrix(V, W)
+        AT_SM = A_SM.transpose()
+        A = MatrixFreeLinearOperator(domain=V, codomain=W, dot=lambda v: A_SM @ v, dot_transpose=lambda v: AT_SM @ v)
+
+        # example 2: a truly matrix-free linear operator
+        A = MatrixFreeLinearOperator(domain=V, codomain=V, dot=lambda v: 2*v, dot_transpose=lambda v: 2*v)
+
+    """
+
+    def __init__(self, domain, codomain, dot, dot_transpose=None):
+
+        assert isinstance(domain, VectorSpace)
+        assert isinstance(codomain, VectorSpace)       
+        assert isinstance(dot, LambdaType)
+
+        self._domain = domain
+        self._codomain = codomain
+        self._dot = dot
+
+        sig =  signature(dot)
+        self._dot_takes_out_arg = ('out' in [p.name for p in sig.parameters.values() if p.kind == p.KEYWORD_ONLY])
+
+        if dot_transpose is not None:
+            assert isinstance(dot_transpose, LambdaType)
+            self._dot_transpose = dot_transpose
+            sig =  signature(dot_transpose)
+            self._dot_transpose_takes_out_arg = ('out' in [p.name for p in sig.parameters.values() if p.kind == p.KEYWORD_ONLY])
+        else:
+            self._dot_transpose = None
+            self._dot_transpose_takes_out_arg = False            
+
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def codomain(self):
+        return self._codomain
+
+    @property
+    def dtype(self):
+        return None
+
+    def dot(self, v, out=None):
+        assert isinstance(v, Vector)
+        assert v.space == self.domain
+
+        if out is not None:
+            assert isinstance(out, Vector)
+            assert out.space == self.codomain
+        else:
+            out = self.codomain.zeros()
+
+        if self._dot_takes_out_arg:
+            self._dot(v, out=out)
+        else:
+            # provided dot product does not take an out argument: we simply copy the result into out
+            self._dot(v).copy(out=out)
+
+        return out
+
+    def toarray(self):
+        raise NotImplementedError('toarray() is not defined for MatrixFreeLinearOperator.')
+
+    def tosparse(self):
+        raise NotImplementedError('tosparse() is not defined for MatrixFreeLinearOperator.')
+
+    def transpose(self, conjugate=False):
+        if self._dot_transpose is None:
+            raise NotImplementedError('no transpose dot method was given -- cannot create the transpose operator')
+
+        if conjugate:
+            if self._dot_transpose_takes_out_arg:
+                new_dot = lambda v, out=None: self._dot_transpose(v, out=out).conjugate()
+            else:
+                new_dot = lambda v: self._dot_transpose(v).conjugate()
+        else:
+            new_dot = self._dot_transpose
+
+        return MatrixFreeLinearOperator(domain=self.codomain, codomain=self.domain, dot=new_dot, dot_transpose=self._dot)

psydac/linalg/solvers.py

@@ -39,7 +39,7 @@ def inverse(A, solver, **kwargs):
     A : psydac.linalg.basic.LinearOperator
         Left-hand-side matrix A of linear system; individual entries A[i,j]
         can't be accessed, but A has 'shape' attribute and provides 'dot(p)'
-        function (i.e. matrix-vector product A*p).
+        function (e.g. a matrix-vector product A*p).
 
     solver : str
         Preferred iterative solver. Options are: 'cg', 'pcg', 'bicg',
@@ -1165,7 +1165,7 @@ def solve(self, b, out=None):
         A.dot(x, out=y)
         y -= b
         y *= -1.0
-        y.copy(out=res_old)
+        y.copy(out=res_old)   # res = b - A*x
 
         beta = sqrt(res_old.dot(res_old))
 
@@ -1193,8 +1193,15 @@ def solve(self, b, out=None):
             print( "+---------+---------------------+")
             template = "| {:7d} | {:19.2e} |"
 
+        # check whether solution is already converged:
+        if beta < tol:
+            istop = 1
+            rnorm = beta
+            if verbose:
+                print( template.format(itn, rnorm ))
 
-        for itn in range(1, maxiter + 1 ):
+        while istop == 0 and itn < maxiter:
+            itn += 1
 
             s = 1.0/beta
             y.copy(out=v)

psydac/linalg/tests/test_linalg.py

@@ -20,7 +20,7 @@ def array_equal(a, b):
 def sparse_equal(a, b):
     return (a.tosparse() != b.tosparse()).nnz == 0
 
-def is_pos_def(A):
+def assert_pos_def(A):
     assert isinstance(A, LinearOperator)
     A_array = A.toarray()
     assert np.all(np.linalg.eigvals(A_array) > 0)
@@ -50,7 +50,7 @@ def get_StencilVectorSpace(n1, n2, p1, p2, P1, P2):
     V = StencilVectorSpace(C)
     return V
 
-def get_positive_definite_stencilmatrix(V):
+def get_positive_definite_StencilMatrix(V):
 
     np.random.seed(2)
     assert isinstance(V, StencilVectorSpace)
@@ -700,9 +700,9 @@ def test_positive_definite_matrix(n1, n2, p1, p2):
     P1 = False
     P2 = False
     V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2)
-    S = get_positive_definite_stencilmatrix(V)
+    S = get_positive_definite_StencilMatrix(V)
 
-    is_pos_def(S)
+    assert_pos_def(S)
 
 #===============================================================================
 @pytest.mark.parametrize('n1', [3, 5])
@@ -745,7 +745,7 @@ def test_operator_evaluation(n1, n2, p1, p2):
     V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2)
 
     # Initiate positive definite StencilMatrices for which the cg inverse works (necessary for certain tests)
-    S = get_positive_definite_stencilmatrix(V)
+    S = get_positive_definite_StencilMatrix(V)
 
     # Initiate StencilVectors 
     v = StencilVector(V)
@@ -769,7 +769,7 @@ def test_operator_evaluation(n1, n2, p1, p2):
 
     ### 2.1 PowerLO test
     Bmat = B.toarray()
-    is_pos_def(B)
+    assert_pos_def(B)
     uarr = u.toarray()
     b0 = ( B**0 @ u ).toarray()
     b1 = ( B**1 @ u ).toarray()
@@ -799,7 +799,7 @@ def test_operator_evaluation(n1, n2, p1, p2):
     assert np.array_equal(zeros, z2)
 
     Smat = S.toarray()
-    is_pos_def(S)
+    assert_pos_def(S)
     varr = v.toarray()
     s0 = ( S**0 @ v ).toarray()
     s1 = ( S**1 @ v ).toarray()
@@ -960,8 +960,8 @@ def test_x0update(solver):
     P1 = False
     P2 = False
     V = get_StencilVectorSpace(n1, n2, p1, p2, P1, P2)
-    A = get_positive_definite_stencilmatrix(V)
-    is_pos_def(A)
+    A = get_positive_definite_StencilMatrix(V)
+    assert_pos_def(A)
     b = StencilVector(V)
     for n in range(n1):
         b[n, :] = 1.