diff --git a/.gitignore b/.gitignore index 793b6202..5c39733b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ *.pyc +*.so *.sublime-project *.sublime-workspace 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 9b5cc0e7..f8c2964d 100644 --- a/SimPEG/tests/test_Solver.py +++ b/SimPEG/tests/test_Solver.py @@ -58,23 +58,39 @@ 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_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(self): + def test_directLower_M_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,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) @@ -82,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)) @@ -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..3efaf10b 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -1,7 +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: + 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): """ @@ -55,11 +70,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 @@ -69,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. @@ -78,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: @@ -104,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. @@ -114,21 +131,29 @@ 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) 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': + 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])): + 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'): + def solveForward(self, b, backend=None): """ Use solve instead of this interface. @@ -138,21 +163,29 @@ 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) 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': + 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]): + 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'): + def solveDiagonal(self, b, backend=None): """ Use solve instead of this interface. @@ -162,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 @@ -205,3 +240,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..ab9048d6 --- /dev/null +++ b/SimPEG/utils/TriSolve.f @@ -0,0 +1,64 @@ +c File TriSolve.f + 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,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 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 + end do + end subroutine forward + + + 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,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 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 + end do + + end subroutine backward diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index c54949ad..3f82ba6d 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,14 +3,13 @@ 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 from ipythonUtils import easyAnimate as animate - +import Solver +from Solver import Solver def setKwargs(obj, **kwargs): """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""