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