Solver updates. Remove .solve() replace with a multiplication.

SimPEG/Mesh/TreeMesh.py

@@ -1127,7 +1127,7 @@ if __name__ == '__main__':
     q = np.zeros(M.nC)
     q[208] = -1.0
     q[291] = 1.0
-    b = Solver(-DIV*Mf*DIV.T).solve(q)
+    b = Solver(-DIV*Mf*DIV.T) * (q)
     plt.figure()
     M.plotImage(b)
     # plt.gca().invert_yaxis()

SimPEG/Mesh/View.py

@@ -186,7 +186,7 @@ class TensorView(object):
             q[[4,4],[4,4],[2,6]]=[-1,1]
             q = Utils.mkvc(q)
             A = M.faceDiv*M.cellGrad
-            b = Solver(A).solve(q)
+            b = Solver(A) * (q)
             M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
 
         """
@@ -628,7 +628,7 @@ if __name__ == '__main__':
     q[[4,4],[4,4],[2,6]]=[-1,1]
     q = Utils.mkvc(q)
     A = M.faceDiv*M.cellGrad
-    b = Solver(A).solve(q)
+    b = Solver(A) * (q)
 
     M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8})
     M2 = Mesh.TensorMesh([10,20],x0=[10,5])

SimPEG/Optimization.py

@@ -1,5 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp
-from Solver import Solver, SolverCG
+from Utils.SolverUtils import Solver, SolverCG
 norm = np.linalg.norm
 
 
@@ -692,7 +692,7 @@ class BFGS(Minimize, Remember):
     def bfgsrec(self,k,n,nn,S,Y,d):
         """BFGS recursion"""
         if k < 0:
-            d = self.bfgsH0.solve(d)
+            d = self.bfgsH0 * (d)
         else:
             khat    = 0 if nn is 0 else np.mod(n-nn+k,nn)
             gamma   = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])
@@ -731,7 +731,7 @@ class GaussNewton(Minimize, Remember):
 
     @Utils.timeIt
     def findSearchDirection(self):
-        return Solver(self.H).solve(-self.g)
+        return Solver(self.H) * (-self.g)
 
 
 class InexactGaussNewton(BFGS, Minimize, Remember):
@@ -779,7 +779,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
     @Utils.timeIt
     def findSearchDirection(self):
         Hinv = SolverCG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG)
-        p = Hinv.solve(-self.g)
+        p = Hinv * (-self.g)
         return p
 
 
@@ -846,7 +846,7 @@ class NewtonRoot(object):
             r, J = fun(x, return_g=True)
 
             Jinv = self.Solver(J, **self.solverOpts)
-            dh   = - Jinv.solve(r)
+            dh   = - (Jinv * r)
 
             muLS = 1.
             LScnt  = 1

SimPEG/Solver.py

@@ -1,6 +0,0 @@
-import scipy.sparse as sp
-from SimPEG.Utils import SolverUtils
-
-Solver   = SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False)
-SolverLU = SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
-SolverCG = SolverUtils.ISolverWrap(sp.linalg.cg)

SimPEG/Tests/test_boundaryPoisson.py

@@ -55,12 +55,12 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
         elif self.myTest == 'xc':
             #TODO: fix the null space
             solver = SolverCG(A, maxiter=1000)
-            xc = solver.solve(rhs)
+            xc = solver * (rhs)
             print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
             err = np.linalg.norm((xc-xc_anal), np.inf)
         elif self.myTest == 'xcJ':
             #TODO: fix the null space
-            xc = Solver(A).solve(rhs)
+            xc = Solver(A) * (rhs)
             print np.linalg.norm(Utils.mkvc(A*xc) - rhs)
             j = McI*(G*xc + P*phi_bc)
             err = np.linalg.norm((j-j_anal), np.inf)
@@ -140,10 +140,10 @@ class Test2D_InhomogeneousDirichlet(OrderTest):
         elif self.myTest == 'q':
             err = np.linalg.norm((q-q_anal), np.inf)
         elif self.myTest == 'xc':
-            xc = Solver(A).solve(rhs)
+            xc = Solver(A) * (rhs)
             err = np.linalg.norm((xc-xc_anal), np.inf)
         elif self.myTest == 'xcJ':
-            xc = Solver(A).solve(rhs)
+            xc = Solver(A) * (rhs)
             j = McI*(G*xc + P*bc)
             err = np.linalg.norm((j-j_anal), np.inf)
 

SimPEG/Tests/test_tensorMesh.py

@@ -79,7 +79,7 @@ class TestPoissonEqn(OrderTest):
             err = np.linalg.norm((sA - sN), np.inf)
         else:
             fA = fun(self.M.gridCC)
-            fN = Solver(D*G).solve(sol(self.M.gridCC))
+            fN = Solver(D*G) * (sol(self.M.gridCC))
             err = np.linalg.norm((fA - fN), np.inf)
         return err
 

SimPEG/Utils/SolverUtils.py

@@ -1,4 +1,4 @@
-import numpy as np
+import numpy as np, scipy.sparse as sp
 from matutils import mkvc
 import warnings
 
@@ -13,7 +13,16 @@ def _checkAccuracy(A, b, X, accuracyTol):
         warnings.warn(msg, RuntimeWarning)
 
 
-def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
+def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
+    """
+    Wraps a direct Solver.
+
+    ::
+
+        Solver   = SolverUtils.SolverWrapD(sp.linalg.spsolve, factorize=False)
+        SolverLU = SolverUtils.SolverWrapD(sp.linalg.splu, factorize=True)
+
+    """
 
     def __init__(self, A, **kwargs):
         self.A = A.tocsc()
@@ -21,7 +30,10 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
         if factorize:
             self.solver = fun(self.A, **kwargs)
 
-    def solve(self, b):
+    def __mul__(self, b):
+        if type(b) is not np.ndarray:
+            raise TypeError('Can only multiply by a numpy array.')
+
         if len(b.shape) == 1 or b.shape[1] == 1:
             b = b.flatten()
             # Just one RHS
@@ -45,22 +57,28 @@ def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
         if factorize and hasattr(self.solver, 'clean'):
             return self.solver.clean()
 
-    def __mul__(self, val):
-        if type(val) is np.ndarray:
-            return self.solve(val)
-        raise TypeError('Can only multiply by a numpy array.')
-
-    return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean, "__mul__": __mul__})
+    return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
 
 
-def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5):
+def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
+    """
+    Wraps an iterative Solver.
+
+    ::
+
+        SolverCG = SolverUtils.SolverWrapI(sp.linalg.cg)
+
+    """
 
     def __init__(self, A, **kwargs):
         self.A = A
         self.kwargs = kwargs
 
-    def solve(self, b):
+    def __mul__(self, b):
+        if type(b) is not np.ndarray:
+            raise TypeError('Can only multiply by a numpy array.')
+
         if len(b.shape) == 1 or b.shape[1] == 1:
             b = b.flatten()
             # Just one RHS
@@ -90,9 +108,9 @@ def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5):
         if hasattr(self.solver, 'clean'):
             return self.solver.clean()
 
-    def __mul__(self, val):
-        if type(val) is np.ndarray:
-            return self.solve(val)
-        raise TypeError('Can only multiply by a numpy array.')
+    return type(fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
-    return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean, "__mul__": __mul__})
+
+Solver   = SolverWrapD(sp.linalg.spsolve, factorize=False)
+SolverLU = SolverWrapD(sp.linalg.splu, factorize=True)
+SolverCG = SolverWrapI(sp.linalg.cg)

SimPEG/__init__.py

@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.sparse as sp
 import Utils
-from Solver import *
+from Utils.SolverUtils import *
 import Mesh
 import Maps
 import Problem

docs/api_Solver.rst

@@ -4,7 +4,43 @@
 Solver
 ******
 
-.. automodule:: SimPEG.Solver
-    :members:
-    :undoc-members:
+The numerical linear algebra solver that you use will ultimately be the
+bottleneck of your large scale inversion. To be the most flexible, SimPEG
+provides wrappers rather than a comprehensive set of solvers (i.e. BYOS).
+
+The interface is as follows::
+
+    A # Where A is a sparse matrix (or linear operator)
+    Ainv = Solver(A, **solverOpts) # Create a solver object with key word arguments
+    x = Ainv * b # Where b is a numpy array of shape (n,) or (n,*)
+    Ainv.clean() # This cleans the memory footprint (if any)
+
+.. note::
+
+    This is somewhat an abuse of notation for solvers as we never actually
+    create A inverse. Instead we are creating an object that acts like A
+    inverse, whether that be a Krylov subspace solver or an LU decomposition.
+
+To wrap up solvers in scipy.sparse.linalg it takes one line of code::
+
+
+    Solver   = SolverWrapD(sp.linalg.spsolve, factorize=False)
+    SolverLU = SolverWrapD(sp.linalg.splu, factorize=True)
+    SolverCG = SolverWrapI(sp.linalg.cg)
+
+.. note::
+
+    The above solvers are loaded into the base name space of SimPEG.
+
+.. seealso::
+
+    https://github.com/rowanc1/pymatsolver
+
+
+The API
+=======
+
+.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
+
+.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI