wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-12 07:12:45 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fortran Solvers. Faster than Matlab!!

Browse Source 
Need to implement multiple RHSs.
This commit is contained in:
Rowan Cockett
2013-11-10 17:38:23 -08:00
parent 059ccd40b3
commit ea5dc21517
3 changed files with 136 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/tests/test_Solver.py
+37 -4
View File
@@ -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={})
SimPEG/utils/Solver.py
+45 -15
View File
@@ -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
SimPEG/utils/TriSolve.f
+54
View File
@@ -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