Added forward and backwards solvers implemented in python. Added tests for direct solvers.

SimPEG/Solver.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.sparse as sparse
 import scipy.sparse.linalg as linalg
 
 
@@ -64,10 +65,36 @@ class Solver(object):
         pass
 
     def solveBackward(self, b):
-        pass
+        "Perform a backwards solve with upper triangular A in CSR format."
+        if type(self.A) is not sparse.csr.csr_matrix:
+            from scipy.sparse import csr_matrix
+            self.A = csr_matrix(self.A)
+        vals = self.A.data
+        rowptr = self.A.indptr
+        colind = self.A.indices
+        x = np.empty_like(b)   # empty() is faster than zeros().
+        for i in reversed(xrange(self.A.shape[0])):
+            ith_row = vals[rowptr[i] : rowptr[i+1]]
+            cols = colind[rowptr[i] : rowptr[i+1]]
+            x_vals = x[cols]
+            x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0]
+        return x
 
     def solveForward(self, b):
-        pass
+        "Perform a forward solve with lower triangular A in CSR format."
+        if type(self.A) is not sparse.csr.csr_matrix:
+            from scipy.sparse import csr_matrix
+            self.A = csr_matrix(self.A)
+        vals = self.A.data
+        rowptr = self.A.indptr
+        colind = self.A.indices
+        x = np.empty_like(b)   # empty() is faster than zeros().
+        for i in xrange(self.A.shape[0]):
+            ith_row = vals[rowptr[i] : rowptr[i+1]]
+            cols = colind[rowptr[i] : rowptr[i+1]]
+            x_vals = x[cols]
+            x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1]
+        return x
 
     def solveDiagonal(self, b):
         diagA = self.A.diagonal()
@@ -96,16 +123,20 @@ if __name__ == '__main__':
     G = M.cellGrad
     Msig = M.getFaceMass()
     A = D*Msig*G
+    A[0,0] *= 10 # remove the constant null space from the matrix
 
-    rhs = np.random.rand(M.nC)
-
+    e = np.ones(M.nC)
+    rhs = A.dot(e)
 
     tic = time()
     solve = Solver(A, options={'factorize':True})
     x = solve.solve(rhs)
     print 'Factorized', time() - tic
+    print np.linalg.norm(e-x,np.inf)
     tic = time()
     solve = Solver(A, options={'factorize':False})
     x = solve.solve(rhs)
     print 'spsolve', time() - tic
+    print np.linalg.norm(e-x,np.inf)
+
 

SimPEG/tests/test_Solver.py

@@ -0,0 +1,111 @@
+import unittest
+from SimPEG import Solver
+from SimPEG.mesh import TensorMesh
+from SimPEG.utils import sdiag
+import numpy as np
+import scipy.sparse as sparse
+
+TOL = 1e-10
+numRHS = 5
+
+
+class TestSolver(unittest.TestCase):
+
+    def setUp(self):
+        h1 = np.ones(10)*100.
+        h2 = np.ones(10)*100.
+        h3 = np.ones(10)*100.
+
+        h = [h1,h2,h3]
+
+        M = TensorMesh(h)
+
+        D = M.faceDiv
+        G = M.cellGrad
+        Msig = M.getFaceMass()
+        A = D*Msig*G
+        A[0,0] *= 10 # remove the constant null space from the matrix
+
+        self.A = A
+        self.M = M
+
+    def test_directFactored_1(self):
+        solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
+        e = np.ones(self.M.nC)
+        rhs = self.A.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+
+    def test_directFactored_M(self):
+        solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = self.A.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directSpsolve_1(self):
+        solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'})
+        e = np.ones(self.M.nC)
+        rhs = self.A.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directSpsolve_M(self):
+        solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'})
+        e = np.ones((self.M.nC, numRHS))
+        rhs = self.A.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directLower_1(self):
+        AL = sparse.tril(self.A)
+        solve = Solver(AL, doDirect=True, flag='L', options={})
+        e = np.ones(self.M.nC)
+        rhs = AL.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directLower_M(self):
+        AL = sparse.tril(self.A)
+        solve = Solver(AL, doDirect=True, flag='L', options={})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = AL.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directUpper_1(self):
+        AU = sparse.triu(self.A)
+        solve = Solver(AU, doDirect=True, flag='U', options={})
+        e = np.ones(self.M.nC)
+        rhs = AU.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directUpper_M(self):
+        AU = sparse.triu(self.A)
+        solve = Solver(AU, doDirect=True, flag='U', options={})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = AU.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directDiagonal_1(self):
+        AD = sdiag(self.A.diagonal())
+        solve = Solver(AD, doDirect=True, flag='D', options={})
+        e = np.ones(self.M.nC)
+        rhs = AD.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directDiagonal_M(self):
+        AD = sdiag(self.A.diagonal())
+        solve = Solver(AD, doDirect=True, flag='D', options={})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = AD.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+
+if __name__ == '__main__':
+    unittest.main()