From ea5dc21517b02d30b8a8c6e2f13876711ce7e8ec Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 10 Nov 2013 17:38:23 -0800 Subject: [PATCH 1/2] Fortran Solvers. Faster than Matlab!! Need to implement multiple RHSs. --- SimPEG/tests/test_Solver.py | 41 ++++++++++++++++++++++--- SimPEG/utils/Solver.py | 60 +++++++++++++++++++++++++++---------- SimPEG/utils/TriSolve.f | 54 +++++++++++++++++++++++++++++++++ 3 files changed, 136 insertions(+), 19 deletions(-) create mode 100644 SimPEG/utils/TriSolve.f diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py index 9b5cc0e7..bfa57761 100644 --- a/SimPEG/tests/test_Solver.py +++ b/SimPEG/tests/test_Solver.py @@ -58,17 +58,33 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directLower_1(self): + def test_directLower_1_fortran(self): AL = sparse.tril(self.A) - solve = Solver(AL, doDirect=True, flag='L', options={}) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) 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): + # def test_directLower_M_fortran(self): + # AL = sparse.tril(self.A) + # solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) + # 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_directLower_1_python(self): AL = sparse.tril(self.A) - solve = Solver(AL, doDirect=True, flag='L', options={}) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'}) + 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_python(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'}) e = np.ones((self.M.nC,numRHS)) rhs = AL.dot(e) x = solve.solve(rhs) @@ -90,6 +106,23 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_1_fortran(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) + 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_fortran(self): + # AU = sparse.triu(self.A) + # solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) + # 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={}) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index d9ac4a1d..1cc2c743 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -2,6 +2,13 @@ import numpy as np import scipy.sparse as sparse import scipy.sparse.linalg as linalg +try: + import TriSolve +except Exception, e: + import os + # Note: this may not work from SublimeText, if that is the case, just run the command in your shell. + os.system('f2py -c TriSolve.f -m TriSolve') + import TriSolve class Solver(object): """ @@ -55,11 +62,11 @@ class Solver(object): elif self.flag is None and not self.doDirect: return self.solveIter(b, **self.options) elif self.flag == 'U': - return self.solveBackward(b) + return self.solveBackward(b, **self.options) elif self.flag == 'L': - return self.solveForward(b) + return self.solveForward(b, **self.options) elif self.flag == 'D': - return self.solveDiagonal(b) + return self.solveDiagonal(b, **self.options) else: raise Exception('Unknown flag.') pass @@ -120,12 +127,15 @@ class Solver(object): 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] + if backend == 'fortran': + x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size) + elif backend == 'python': + 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, backend='python'): @@ -144,12 +154,15 @@ class Solver(object): 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] + if backend == 'fortran': + x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size) + elif backend == 'python': + 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, backend='python'): @@ -205,3 +218,20 @@ if __name__ == '__main__': print np.linalg.norm(e-x,np.inf) + n = 6000 + A_dense = np.random.random((n,n)) + L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned. + e = np.ones(n) + b = np.dot(L, e) + + A = sparse.csr_matrix(L) + pSolve = Solver(A,flag='L',options={'backend':'python'}); + fSolve = Solver(A,flag='L',options={'backend':'fortran'}) + tic = time() + x = pSolve.solve(b) + toc = time() - tic + print 'Error Forward Python = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc + tic = time() + x = fSolve.solve(b) + toc = time() - tic + print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc diff --git a/SimPEG/utils/TriSolve.f b/SimPEG/utils/TriSolve.f new file mode 100644 index 00000000..6d15d013 --- /dev/null +++ b/SimPEG/utils/TriSolve.f @@ -0,0 +1,54 @@ +c File TriSolve.f + subroutine forward(al, ial, jal, b, nv, n, x) + double precision al(nv) + integer ial(n+1) + integer jal(nv) + double precision b(n) + double precision x(n) + integer nv + integer n +cf2py intent(in) :: al +cf2py intent(in) :: ial +cf2py intent(in) :: jal +cf2py intent(in) :: b +cf2py intent(in) :: nv +cf2py intent(in) :: n +cf2py intent(out) :: x + real ( kind = 8 ) t + + do k = 1, n + t = b(k) + do j = ial(k)+1, ial(k+1) + t = t - al(j) * x(jal(j)+1) + end do + x(k) = t/al(ial(k+1)) + end do + end subroutine forward + + + subroutine backward(au,iau, jau, b, nv, n, x) + double precision au(nv) + integer iau(n+1) + integer jau(nv) + double precision b(n) + double precision x(n) + integer nv + integer n +cf2py intent(in) :: au +cf2py intent(in) :: iau +cf2py intent(in) :: jau +cf2py intent(in) :: b +cf2py intent(in) :: nv +cf2py intent(in) :: n +cf2py intent(out) :: x + real ( kind = 8 ) t + + do k = n, 1, -1 + t = b(k) + do j = iau(k)+1, iau(k+1) + t = t - au(j) * x(jau(j)+1) + end do + x(k) = t/au(iau(k)+1) + end do + + end subroutine backward From d3f38047e40a18899e0be98e70bd537d85beecbc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 12 Nov 2013 10:36:20 -0800 Subject: [PATCH 2/2] Multiple RHSs on solvers in Fortran. ~2x speed up on matlab implementation for a single RHS. for multiple RHS there are still some problems. Someone with some knowledge of how fortran works should look at this code. Added a setup.py script that complies things. f2py should work on most computers, because it is included in the numpy distribution. --- .gitignore | 3 +- SimPEG/setup.py | 8 + SimPEG/tests/api_TestResults.rst | 355 ------------------------------- SimPEG/tests/test_Solver.py | 50 ++--- SimPEG/utils/Solver.py | 42 +++- SimPEG/utils/TriSolve.f | 42 ++-- SimPEG/utils/__init__.py | 4 +- 7 files changed, 95 insertions(+), 409 deletions(-) create mode 100644 SimPEG/setup.py delete mode 100644 SimPEG/tests/api_TestResults.rst diff --git a/.gitignore b/.gitignore index bcfb3d73..abbec73c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.pyc SimPEG.sublime-project SimPEG.sublime-workspace -docs/_build/ \ No newline at end of file +docs/_build/ +myNotebooks/* diff --git a/SimPEG/setup.py b/SimPEG/setup.py new file mode 100644 index 00000000..c421cb4b --- /dev/null +++ b/SimPEG/setup.py @@ -0,0 +1,8 @@ +import os +print 'Compiling TriSolve.' +os.system('f2py -c utils/TriSolve.f -m TriSolve') +print 'TriSolve Compiled! yay.' +print 'Moving TriSolve into Utils.' +os.system('mv TriSolve.so utils/TriSolve.so') +print 'Thats it. Well Done Computer.' + diff --git a/SimPEG/tests/api_TestResults.rst b/SimPEG/tests/api_TestResults.rst deleted file mode 100644 index 0d63a328..00000000 --- a/SimPEG/tests/api_TestResults.rst +++ /dev/null @@ -1,355 +0,0 @@ -.. _api_TestResults: - -.. raw:: html - - - - - -
-

Test Report

-

Start Time: 2013-11-05 15:24:44

-

Duration: 0:00:00.007500

-

Status: Pass 22

- -

This demonstrates the report output by Prasanna.Yelsangikar.

-
- - - -

Show - Summary - Failed - All -

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
Total222200 
- diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py index bfa57761..f8c2964d 100644 --- a/SimPEG/tests/test_Solver.py +++ b/SimPEG/tests/test_Solver.py @@ -58,22 +58,6 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directLower_1_fortran(self): - AL = sparse.tril(self.A) - solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) - 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_fortran(self): - # AL = sparse.tril(self.A) - # solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) - # 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_directLower_1_python(self): AL = sparse.tril(self.A) solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'}) @@ -88,9 +72,25 @@ class TestSolver(unittest.TestCase): e = np.ones((self.M.nC,numRHS)) rhs = AL.dot(e) x = solve.solve(rhs) + + def test_directLower_1_fortran(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) + 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_directUpper_1(self): + def test_directLower_M_fortran(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'}) + 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) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_1_python(self): AU = sparse.triu(self.A) solve = Solver(AU, doDirect=True, flag='U', options={}) e = np.ones(self.M.nC) @@ -98,7 +98,7 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - def test_directUpper_M(self): + def test_directUpper_M_python(self): AU = sparse.triu(self.A) solve = Solver(AU, doDirect=True, flag='U', options={}) e = np.ones((self.M.nC,numRHS)) @@ -115,13 +115,13 @@ class TestSolver(unittest.TestCase): x = solve.solve(rhs) self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) - # def test_directUpper_M_fortran(self): - # AU = sparse.triu(self.A) - # solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) - # 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_directUpper_M_fortran(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'}) + 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()) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 1cc2c743..3efaf10b 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -1,14 +1,22 @@ import numpy as np import scipy.sparse as sparse import scipy.sparse.linalg as linalg +from SimPEG.utils import mkvc + +DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} try: import TriSolve except Exception, e: - import os - # Note: this may not work from SublimeText, if that is the case, just run the command in your shell. - os.system('f2py -c TriSolve.f -m TriSolve') - import TriSolve + try: + import os + # Note: this may not work from SublimeText, if that is the case, just run the command in your shell. + os.system('f2py -c TriSolve.f -m TriSolve') + import TriSolve + except Exception, e: + print 'Warning: Python backend is being used for solver. Run setup.py from the command line.' + DEFAULTS['forward'] = 'python' + DEFAULTS['backward'] = 'python' class Solver(object): """ @@ -76,7 +84,7 @@ class Solver(object): del self.dsolve self.dsolve = None - def solveDirect(self, b, factorize=False, backend='scipy'): + def solveDirect(self, b, factorize=False, backend=None): """ Use solve instead of this interface. @@ -85,6 +93,8 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['scipy'] + assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' if factorize and self.dsolve is None: @@ -111,7 +121,7 @@ class Solver(object): def solveIter(self, b, M=None, iterSolver='CG'): pass - def solveBackward(self, b, backend='python'): + def solveBackward(self, b, backend=None): """ Use solve instead of this interface. @@ -121,6 +131,7 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['backward'] if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) @@ -128,7 +139,11 @@ class Solver(object): rowptr = self.A.indptr colind = self.A.indices if backend == 'fortran': - x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size) + if len(b.shape) == 1 or b.shape[1] == 1: + x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size, 1) + x = mkvc(x) + else: + x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1]) elif backend == 'python': x = np.empty_like(b) # empty() is faster than zeros(). for i in reversed(xrange(self.A.shape[0])): @@ -138,7 +153,7 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] return x - def solveForward(self, b, backend='python'): + def solveForward(self, b, backend=None): """ Use solve instead of this interface. @@ -148,6 +163,7 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['forward'] if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) @@ -155,7 +171,11 @@ class Solver(object): rowptr = self.A.indptr colind = self.A.indices if backend == 'fortran': - x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size) + if len(b.shape) == 1 or b.shape[1] == 1: + x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size, 1) + x = mkvc(x) + else: + x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1]) elif backend == 'python': x = np.empty_like(b) # empty() is faster than zeros(). for i in xrange(self.A.shape[0]): @@ -165,7 +185,7 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] return x - def solveDiagonal(self, b, backend='python'): + def solveDiagonal(self, b, backend=None): """ Use solve instead of this interface. @@ -175,6 +195,8 @@ class Solver(object): :rtype: numpy.ndarray :return: x """ + if backend is None: backend = DEFAULTS['diagonal'] + diagA = self.A.diagonal() if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS diff --git a/SimPEG/utils/TriSolve.f b/SimPEG/utils/TriSolve.f index 6d15d013..ab9048d6 100644 --- a/SimPEG/utils/TriSolve.f +++ b/SimPEG/utils/TriSolve.f @@ -1,54 +1,64 @@ c File TriSolve.f - subroutine forward(al, ial, jal, b, nv, n, x) + subroutine forward(al, ial, jal, b, nv, n, nRHS, x) double precision al(nv) integer ial(n+1) integer jal(nv) - double precision b(n) - double precision x(n) + double precision b(n,nRHS) + double precision x(n,nRHS) integer nv integer n + integer nRHS + integer rhs cf2py intent(in) :: al cf2py intent(in) :: ial cf2py intent(in) :: jal cf2py intent(in) :: b cf2py intent(in) :: nv cf2py intent(in) :: n +cf2py intent(in) :: nRHS cf2py intent(out) :: x real ( kind = 8 ) t - do k = 1, n - t = b(k) - do j = ial(k)+1, ial(k+1) - t = t - al(j) * x(jal(j)+1) + do rhs = 1, nRHS + do k = 1, n + t = b(k,rhs) + do j = ial(k)+1, ial(k+1) + t = t - al(j) * x(jal(j)+1,rhs) + end do + x(k,rhs) = t/al(ial(k+1)) end do - x(k) = t/al(ial(k+1)) end do end subroutine forward - subroutine backward(au,iau, jau, b, nv, n, x) + subroutine backward(au,iau, jau, b, nv, n, nRHS, x) double precision au(nv) integer iau(n+1) integer jau(nv) - double precision b(n) - double precision x(n) + double precision b(n,nRHS) + double precision x(n,nRHS) integer nv integer n + integer nRHS + integer rhs cf2py intent(in) :: au cf2py intent(in) :: iau cf2py intent(in) :: jau cf2py intent(in) :: b cf2py intent(in) :: nv cf2py intent(in) :: n +cf2py intent(in) :: nRHS cf2py intent(out) :: x real ( kind = 8 ) t - do k = n, 1, -1 - t = b(k) - do j = iau(k)+1, iau(k+1) - t = t - au(j) * x(jau(j)+1) + do rhs = 1, nRHS + do k = n, 1, -1 + t = b(k,rhs) + do j = iau(k)+1, iau(k+1) + t = t - au(j) * x(jau(j)+1,rhs) + end do + x(k,rhs) = t/au(iau(k)+1) end do - x(k) = t/au(iau(k)+1) end do end subroutine backward diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index ae83222e..d3955170 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,9 +3,9 @@ import sputils import lomutils import interputils import ModelBuilder -import Solver -from Solver import Solver from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird from interputils import interpmat +import Solver +from Solver import Solver