From c33c9bee55e65485f36b8c116ff38369342bd9cb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 30 Oct 2013 22:42:25 -0600 Subject: [PATCH] Added forward and backwards solvers implemented in python. Added tests for direct solvers. --- SimPEG/Solver.py | 39 +++++++++++-- SimPEG/tests/test_Solver.py | 111 ++++++++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 4 deletions(-) create mode 100644 SimPEG/tests/test_Solver.py diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index cbd4a3ea..a9105c84 100644 --- a/SimPEG/Solver.py +++ b/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) + diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py new file mode 100644 index 00000000..9b5cc0e7 --- /dev/null +++ b/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()