mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-18 12:40:30 +08:00
Merge branch 'master' of https://bitbucket.org/rcockett/simpeg into BoundConstraint
Conflicts: .gitignore SimPEG/utils/__init__.py
This commit is contained in:
Rowan Cockett
+71
-19
@@ -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
|
||||
|
||||
@@ -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
|
||||
@@ -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."""
|
||||
|
||||
Reference in New Issue
Block a user