Merge branch 'release/0.1'

Author: dpo

examples/bmark.py

@@ -3,12 +3,13 @@
 # The test matrix may be obtained from http://math.nist.gov/MatrixMarket
 
 import numpy as np
+from pykrylov.linop import PysparseLinearOperator
 from pykrylov.cgs import CGS
 from pykrylov.tfqmr import TFQMR
 from pykrylov.bicgstab import BiCGSTAB
 from pysparse import spmatrix
 from pysparse.sparse.pysparseMatrix import PysparseMatrix as sp
-from math import sqrt
+import sys
 
 class DiagonalPrec:
 
@@ -30,9 +31,9 @@ def __call__(self, y, **kwargs):
     print hdr
     print '-' * len(hdr)
 
-    #AA = spmatrix.ll_mat_from_mtx('mcca.mtx')
-    AA = spmatrix.ll_mat_from_mtx('jpwh_991.mtx')
+    AA = spmatrix.ll_mat_from_mtx(sys.argv[1])
     A = sp(matrix=AA)
+    op = PysparseLinearOperator(A)
 
     # Create diagonal preconditioner
     dp = DiagonalPrec(A)
@@ -42,13 +43,13 @@ def __call__(self, y, **kwargs):
     rhs = A*e
 
     for KSolver in [CGS, TFQMR, BiCGSTAB]:
-        ks = KSolver( lambda v: A*v,
+        ks = KSolver( op,
                       #precon = dp,
                       #verbose=False,
                       reltol = 1.0e-8
                       )
-        ks.solve(rhs, guess = 1+np.arange(n, dtype=np.float), matvec_max=2*n)
+        ks.solve(rhs, guess = 1+np.arange(n, dtype=op.dtype), matvec_max=2*n)
 
-        err = np.linalg.norm(ks.bestSolution-e)/sqrt(n)
+        err = np.linalg.norm(ks.bestSolution-e)/np.sqrt(n)
         print fmt % (ks.acronym, ks.nMatvec, ks.residNorm0, ks.residNorm, err)
 

pykrylov/bicgstab/bicgstab.py

@@ -2,7 +2,7 @@
 __docformat__ = 'restructuredtext'
 
 import numpy as np
-from math import sqrt
+
 
 from pykrylov.generic import KrylovMethod
 
@@ -14,9 +14,9 @@ class BiCGSTAB( KrylovMethod ):
 
         A x = b
 
-    where the matrix A is unsymmetric and nonsingular.
+    where the operator A is unsymmetric and nonsingular.
 
-    Bi-CGSTAB requires 2 matrix-vector products, 6 dot products and 6 daxpys
+    Bi-CGSTAB requires 2 operator-vector products, 6 dot products and 6 daxpys
     per iteration.
 
     In addition, if a preconditioner is supplied, it needs to solve 2
@@ -33,8 +33,8 @@ class BiCGSTAB( KrylovMethod ):
                    Computing **13** (2), pp. 631--644, 1992.
     """
 
-    def __init__(self, matvec, **kwargs):
-        KrylovMethod.__init__(self, matvec, **kwargs)
+    def __init__(self, op, **kwargs):
+        KrylovMethod.__init__(self, op, **kwargs)
 
         self.name = 'Bi-Conjugate Gradient Stabilized'
         self.acronym = 'Bi-CGSTAB'
@@ -53,19 +53,20 @@ def solve(self, rhs, **kwargs):
         nMatvec = 0
 
         # Initial guess is zero unless one is supplied
+        result_type = np.result_type(self.op.dtype, rhs.dtype)
         guess_supplied = 'guess' in kwargs.keys()
-        x = kwargs.get('guess', np.zeros(n))
+        x = kwargs.get('guess', np.zeros(n)).astype(result_type)
         matvec_max = kwargs.get('matvec_max', 2*n)
 
         # Initial residual is the fixed vector
         r0 = rhs
         if guess_supplied:
-            r0 = rhs - self.matvec(x)
+            r0 = rhs - self.op * x
             nMatvec += 1
 
         rho = alpha = omega = 1.0
         rho_next = np.dot(r0,r0)
-        residNorm = self.residNorm0 = sqrt(rho_next)
+        residNorm = self.residNorm0 = np.abs(np.sqrt(rho_next))
         threshold = max( self.abstol, self.reltol * self.residNorm0 )
 
         finished = (residNorm <= threshold or nMatvec >= matvec_max)
@@ -78,8 +79,8 @@ def solve(self, rhs, **kwargs):
 
         if not finished:
             r = r0.copy()
-            p = np.zeros(n)
-            v = np.zeros(n)
+            p = np.zeros(n, dtype=result_type)
+            v = np.zeros(n, dtype=result_type)
 
         while not finished:
 
@@ -93,17 +94,17 @@ def solve(self, rhs, **kwargs):
 
             # Compute preconditioned search direction
             if self.precon is not None:
-                q = self.precon(p)
+                q = self.precon * p
             else:
                 q = p
 
-            v = self.matvec(q) ; nMatvec += 1
+            v = self.op * q ; nMatvec += 1
 
             alpha = rho/np.dot(r0, v)
             s = r - alpha * v
 
             # Check for CGS termination
-            residNorm = sqrt(np.dot(s,s))
+            residNorm = np.linalg.norm(s)
 
             self.logger.info('%6d  %8.2e' % (nMatvec, residNorm))
 
@@ -117,11 +118,11 @@ def solve(self, rhs, **kwargs):
                 continue
 
             if self.precon is not None:
-                z = self.precon(s)
+                z = self.precon * s
             else:
                 z = s
 
-            t = self.matvec(z) ; nMatvec += 1
+            t = self.op * z ; nMatvec += 1
             omega = np.dot(t,s)/np.dot(t,t)
             rho_next = -omega * np.dot(r0,t)
 
@@ -135,7 +136,7 @@ def solve(self, rhs, **kwargs):
             x += z
             x += alpha * q
 
-            residNorm = sqrt(np.dot(r,r))
+            residNorm = np.linalg.norm(r)
 
             self.logger.info('%6d  %8.2e' % (nMatvec, residNorm))
 

pykrylov/cg/cg.py

@@ -2,7 +2,6 @@
 __docformat__ = 'restructuredtext'
 
 import numpy as np
-from math import sqrt
 
 from pykrylov.tools.utils import check_symmetric
 from pykrylov.generic import KrylovMethod
@@ -15,24 +14,24 @@ class CG( KrylovMethod ):
 
         A x = b
 
-    where the matrix A is square, symmetric and positive definite. This is
+    where the operator A is square, symmetric and positive definite. This is
     equivalent to solving the unconstrained convex quadratic optimization
     problem
 
         minimize    -<b,x> + 1/2 <x, Ax>
 
     in the variable x.
 
-    CG performs 1 matrix-vector product, 2 dot products and 3 daxpys per
+    CG performs 1 operator-vector product, 2 dot products and 3 daxpys per
     iteration.
 
     If a preconditioner is supplied, it needs to solve one preconditioning
     system per iteration. Our implementation is standard and follows [Kelley]_
     and [Templates]_.
     """
 
-    def __init__(self, matvec, **kwargs):
-        KrylovMethod.__init__(self, matvec, **kwargs)
+    def __init__(self, op, **kwargs):
+        KrylovMethod.__init__(self, op, **kwargs)
 
         self.name = 'Conjugate Gradient'
         self.acronym = 'CG'
@@ -52,9 +51,9 @@ def solve(self, rhs, **kwargs):
         :Keywords:
 
            :guess:           Initial guess (Numpy array). Default: 0.
-           :matvec_max:      Max. number of matrix-vector produts. Default: 2n.
-           :check_symmetric: Ensure matrix is symmetric. Default: False.
-           :check_curvature: Ensure matrix is positive definite. Default: True.
+           :matvec_max:      Max. number of operator-vector produts. Default: 2n.
+           :check_symmetric: Ensure operator is symmetric. Default: False.
+           :check_curvature: Ensure operator is positive definite. Default: True.
            :store_resids:    Store full residual vector history. Default: False.
            :store_iterates:  Store full iterate history. Default: False.
 
@@ -68,13 +67,14 @@ def solve(self, rhs, **kwargs):
         store_iterates = kwargs.get('store_iterates', False)
 
         if check_sym:
-            if not check_symmetric(self.matvec):
-                self.logger.error('Coefficient matrix is not symmetric')
+            if not check_symmetric(self.op):
+                self.logger.error('Coefficient operator is not symmetric')
                 return
 
         # Initial guess
+        result_type = np.result_type(self.op.dtype, rhs.dtype)
         guess_supplied = 'guess' in kwargs.keys()
-        x = kwargs.get('guess', np.zeros(n))
+        x = kwargs.get('guess', np.zeros(n)).astype(result_type)
 
         if store_iterates:
             self.iterates.append(x.copy())
@@ -84,20 +84,20 @@ def solve(self, rhs, **kwargs):
         # Initial residual vector
         r = -rhs
         if guess_supplied:
-            r += self.matvec(x)
+            r += self.op * x
             nMatvec += 1
 
         # Initial preconditioned residual vector
         if self.precon is not None:
-            y = self.precon(r)
+            y = self.precon * r
         else:
             y = r
 
         if store_resids:
             self.resids.append(y.copy())
 
         ry = np.dot(r,y)
-        self.residNorm0 = residNorm = sqrt(ry)
+        self.residNorm0 = residNorm = np.abs(np.sqrt(ry))
         self.residHistory.append(self.residNorm0)
         threshold = max(self.abstol, self.reltol * self.residNorm0)
 
@@ -112,13 +112,13 @@ def solve(self, rhs, **kwargs):
 
         while residNorm > threshold and nMatvec < matvec_max and definite:
 
-            Ap  = self.matvec(p)
+            Ap  = self.op * p
             nMatvec += 1
             pAp = np.dot(p, Ap)
 
             if check_curvature:
-                if pAp <= 0:
-                    self.logger.error('Coefficient matrix is not positive definite')
+                if np.imag(pAp) > 1.0e-8 * np.abs(pAp) or np.real(pAp) <= 0:
+                    self.logger.error('Coefficient operator is not positive definite')
                     self.infiniteDescent = p
                     definite = False
                     continue
@@ -135,7 +135,7 @@ def solve(self, rhs, **kwargs):
 
             # Compute preconditioned residual
             if self.precon is not None:
-                y = self.precon(r)
+                y = self.precon * r
             else:
                 y = r
 
@@ -151,10 +151,10 @@ def solve(self, rhs, **kwargs):
             p -= r
 
             ry = ry_next
-            residNorm = sqrt(ry)
+            residNorm = np.abs(np.sqrt(ry))
             self.residHistory.append(residNorm)
 
-            info = '%6d  %7.1e  %8.1e' % (nMatvec, residNorm, pAp)
+            info = '%6d  %7.1e  %8.1e' % (nMatvec, residNorm, np.real(pAp))
             self.logger.info(info)
 
 

pykrylov/cgs/cgs.py

@@ -2,7 +2,6 @@
 __docformat__ = 'restructuredtext'
 
 import numpy as np
-from math import sqrt
 
 from pykrylov.generic import KrylovMethod
 
@@ -14,10 +13,10 @@ class CGS( KrylovMethod ):
 
         A x = b
 
-    where the matrix A may be unsymmetric.
+    where the operator A may be unsymmetric.
 
-    CGS requires 2 matrix-vector products with A, 3 dot products and 7 daxpys
-    per iteration. It does not require products with the transpose of A.
+    CGS requires 2 operator-vector products with A, 3 dot products and 7 daxpys
+    per iteration. It does not require products with the adjoint of A.
 
     If a preconditioner is supplied, CGS needs to solve two preconditioning
     systems per iteration. The original description appears in [Sonn89]_, which
@@ -31,8 +30,8 @@ class CGS( KrylovMethod ):
                 Computing **10** (1), pp. 36--52, 1989.
     """
 
-    def __init__(self, matvec, **kwargs):
-        KrylovMethod.__init__(self, matvec, **kwargs)
+    def __init__(self, op, **kwargs):
+        KrylovMethod.__init__(self, op, **kwargs)
 
         self.name = 'Conjugate Gradient Squared'
         self.acronym = 'CGS'
@@ -51,16 +50,17 @@ def solve(self, rhs, **kwargs):
         nMatvec = 0
 
         # Initial guess is zero unless one is supplied
+        result_type = np.result_type(self.op.dtype, rhs.dtype)
         guess_supplied = 'guess' in kwargs.keys()
-        x = kwargs.get('guess', np.zeros(n))
+        x = kwargs.get('guess', np.zeros(n)).astype(result_type)
         matvec_max = kwargs.get('matvec_max', 2*n)
 
         r0 = rhs  # Fixed vector throughout
         if guess_supplied:
-            r0 = rhs - self.matvec(x)
+            r0 = rhs - self.op * x
 
         rho = np.dot(r0,r0)
-        residNorm = sqrt(rho)
+        residNorm = np.abs(np.sqrt(rho))
         self.residNorm0 = residNorm
         threshold = max( self.abstol, self.reltol * self.residNorm0 )
         self.logger.info('Initial residual = %8.2e\n' % self.residNorm0)
@@ -76,27 +76,27 @@ def solve(self, rhs, **kwargs):
         while not finished:
 
             if self.precon is not None:
-                y = self.precon(p)
+                y = self.precon * p
             else:
                 y = p
 
-            v = self.matvec(y) ; nMatvec += 1
+            v = self.op * y ; nMatvec += 1
             sigma = np.dot(r0,v)
             alpha = rho/sigma
             q = u - alpha * v
 
             if self.precon is not None:
-                z = self.precon(u+q)
+                z = self.precon * (u+q)
             else:
                 z = u+q
 
             # Update solution and residual
             x += alpha * z
-            Az = self.matvec(z) ; nMatvec += 1
+            Az = self.op * z ; nMatvec += 1
             r -= alpha * Az
 
             # Update residual norm and check convergence
-            residNorm = sqrt(np.dot(r,r))
+            residNorm = np.linalg.norm(r)
 
             if residNorm <= threshold or nMatvec >= matvec_max:
                 finished = True