mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 05:35:13 +08:00
#69 Removed most solvers, replaced by wrappers.
This commit is contained in:
rowanc1
@@ -23,7 +23,7 @@ class LinearProblem(Problem.BaseProblem):
|
||||
return self.G.T.dot(v)
|
||||
|
||||
|
||||
def run(N):
|
||||
def run(N, plotIt=True):
|
||||
mesh = Mesh.TensorMesh([N])
|
||||
|
||||
nk = 20
|
||||
@@ -59,13 +59,14 @@ def run(N):
|
||||
|
||||
mrec = inv.run(m0)
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
plt.figure(1)
|
||||
for i in range(prob.G.shape[0]):
|
||||
plt.plot(prob.G[i,:])
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
plt.figure(1)
|
||||
for i in range(prob.G.shape[0]):
|
||||
plt.plot(prob.G[i,:])
|
||||
|
||||
plt.figure(2)
|
||||
plt.plot(M.vectorCCx, survey.mtrue, 'b-')
|
||||
plt.plot(M.vectorCCx, mrec, 'r-')
|
||||
plt.figure(2)
|
||||
plt.plot(M.vectorCCx, survey.mtrue, 'b-')
|
||||
plt.plot(M.vectorCCx, mrec, 'r-')
|
||||
|
||||
return prob, survey, mesh
|
||||
return prob, survey, mesh, mrec
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
import Utils, numpy as np, scipy.sparse as sp
|
||||
from Solver import Solver
|
||||
from Solver import Solver, SolverCG
|
||||
norm = np.linalg.norm
|
||||
|
||||
|
||||
@@ -778,7 +778,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
|
||||
|
||||
@Utils.timeIt
|
||||
def findSearchDirection(self):
|
||||
Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG})
|
||||
Hinv = SolverCG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG)
|
||||
p = Hinv.solve(-self.g)
|
||||
return p
|
||||
|
||||
|
||||
+4
-403
@@ -1,405 +1,6 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import scipy.sparse.linalg as linalg
|
||||
from Utils import mkvc, sdiag
|
||||
import warnings
|
||||
from SimPEG.Utils import SolverUtils
|
||||
|
||||
DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}
|
||||
OPTIONS = {'direct':['scipy'], 'iter':['scipy'], 'triangular':['python'], 'diagonal':['python']}
|
||||
|
||||
HELP = []
|
||||
|
||||
try:
|
||||
import Utils.TriSolve as TriSolve
|
||||
OPTIONS['triangular'].append('fortran')
|
||||
except Exception, e:
|
||||
HELP += ['Warning: Python backend is being used for solver. Run setup.py from the command line.']
|
||||
DEFAULTS['triangular'] = 'python'
|
||||
|
||||
try:
|
||||
import mumps
|
||||
OPTIONS['direct'].append('mumps')
|
||||
except Exception, e:
|
||||
HELP += ['Warning: mumps solver not available.']
|
||||
|
||||
class Solver(object):
|
||||
"""
|
||||
Solver is a light wrapper on the various types of
|
||||
linear solvers available in python.
|
||||
|
||||
:param scipy.sparse A: Matrix
|
||||
:param bool doDirect: if you want a direct solver
|
||||
:param string flag: Matrix type flag for special solves: [None, 'L', 'U', 'D']
|
||||
:param dict options: options which are passed to each sub solver, see each for details.
|
||||
:rtype: Solver
|
||||
:return: Solver
|
||||
|
||||
To use for direct solvers::
|
||||
|
||||
solve = Solver(A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
|
||||
x = solve.solve(rhs)
|
||||
|
||||
Or in one line::
|
||||
|
||||
x = Solver(A).solve(rhs)
|
||||
|
||||
The flag can be set to None, 'L', 'U', or 'D', for general, lower, upper, and diagonal matrices, respectively.
|
||||
|
||||
"""
|
||||
def __init__(self, A, doDirect=True, flag=None, options={}):
|
||||
assert type(doDirect) is bool, 'doDirect must be a boolean'
|
||||
assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
|
||||
assert type(options) is dict, 'options must be a dictionary object'
|
||||
self.A = A
|
||||
|
||||
self.dsolve = None
|
||||
self.doDirect = doDirect
|
||||
self.flag = flag
|
||||
self.options = options
|
||||
if doDirect: return
|
||||
|
||||
# Now deal with iterative stuff only
|
||||
if 'M' not in options:
|
||||
warnings.warn("You should provide a preconditioner, M.", UserWarning)
|
||||
return
|
||||
M = options['M']
|
||||
if isinstance(M, sp.linalg.LinearOperator):
|
||||
return
|
||||
PreconditionerList = ['J','GS']
|
||||
if type(M) is str:
|
||||
assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']"
|
||||
M = (M,A) # use A as the base for the preconditioner.
|
||||
if type(M) is tuple:
|
||||
assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix."
|
||||
if M[0] is 'J':
|
||||
Jacobi = sdiag(1.0/M[1].diagonal())
|
||||
options['M'] = Jacobi
|
||||
elif M[0] is 'GS':
|
||||
DD = sdiag(M[1].diagonal())
|
||||
Uinv = Solver(M[1], flag='U')
|
||||
Linv = Solver(M[1], flag='L')
|
||||
def GS(f):
|
||||
return Uinv.solve(DD*Linv.solve(f))
|
||||
options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype )
|
||||
|
||||
else:
|
||||
raise Exception('M must be a LinearOperator or a tuple')
|
||||
|
||||
|
||||
def solve(self, b):
|
||||
"""
|
||||
Solves the linear system.
|
||||
|
||||
.. math::
|
||||
|
||||
Ax=b
|
||||
|
||||
:param numpy.ndarray b: the right hand side
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if self.flag is None and self.doDirect:
|
||||
return self.solveDirect(b, **self.options)
|
||||
elif self.flag is None and not self.doDirect:
|
||||
return self.solveIter(b, **self.options)
|
||||
elif self.flag == 'U':
|
||||
return self.solveBackward(b, **self.options)
|
||||
elif self.flag == 'L':
|
||||
return self.solveForward(b, **self.options)
|
||||
elif self.flag == 'D':
|
||||
return self.solveDiagonal(b, **self.options)
|
||||
else:
|
||||
raise Exception('Unknown flag.')
|
||||
pass
|
||||
|
||||
def clean(self):
|
||||
"""Cleans up the memory"""
|
||||
if self.options.has_key('backend'):
|
||||
if self.options['backend'] == 'mumps':
|
||||
self.mctx.destroy()
|
||||
del self.dsolve
|
||||
self.dsolve = None
|
||||
|
||||
def solveDirect(self, b, factorize=False, backend=None):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
:param numpy.ndarray b: the right hand side
|
||||
:param bool factorize: if you want to factorize and store factors
|
||||
:param str backend: which backend to use. Default is scipy
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['direct']
|
||||
|
||||
assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
|
||||
|
||||
if backend == 'scipy':
|
||||
X = self.solveDirect_scipy(b, factorize)
|
||||
elif backend == 'mumps':
|
||||
X = self.solveDirect_mumps(b, factorize)
|
||||
|
||||
return X
|
||||
|
||||
def solveDirect_scipy(self, b, factorize):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
:param numpy.ndarray b: the right hand side
|
||||
:param bool factorize: if you want to factorize and store factors
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if factorize and self.dsolve is None:
|
||||
self.A = self.A.tocsc() # for efficiency
|
||||
self.dsolve = linalg.factorized(self.A)
|
||||
|
||||
if len(b.shape) == 1 or b.shape[1] == 1:
|
||||
# Just one RHS
|
||||
if factorize:
|
||||
return self.dsolve(b.flatten())
|
||||
else:
|
||||
return linalg.dsolve.spsolve(self.A, b)
|
||||
|
||||
# Multiple RHSs
|
||||
X = np.empty_like(b)
|
||||
for i in range(b.shape[1]):
|
||||
if factorize:
|
||||
X[:,i] = self.dsolve(b[:,i])
|
||||
else:
|
||||
X[:,i] = linalg.dsolve.spsolve(self.A,b[:,i])
|
||||
|
||||
return X
|
||||
|
||||
def solveDirect_mumps(self, b, factorize):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
:param numpy.ndarray b: the right hand side
|
||||
:param bool factorize: if you want to factorize and store factors
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if factorize and self.dsolve is None:
|
||||
self.mctx = mumps.DMumpsContext()
|
||||
self.mctx.set_icntl(14, 60)
|
||||
# self.mctx.set_silent()
|
||||
self.mctx.set_centralized_sparse(self.A)
|
||||
self.mctx.run(job=4)
|
||||
|
||||
def mdsolve(rhs):
|
||||
x = rhs.copy()
|
||||
self.mctx.set_rhs(x)
|
||||
self.mctx.run(job=3)
|
||||
return x
|
||||
|
||||
self.dsolve = mdsolve
|
||||
|
||||
if len(b.shape) == 1 or b.shape[1] == 1:
|
||||
# Just one RHS
|
||||
if factorize:
|
||||
X = self.dsolve(b)
|
||||
else:
|
||||
X = mumps.spsolve(self.A, b)
|
||||
|
||||
else:
|
||||
# Multiple RHSs
|
||||
X = np.empty_like(b)
|
||||
for i in range(b.shape[1]):
|
||||
if factorize:
|
||||
X[:,i] = self.dsolve(b[:,i])
|
||||
else:
|
||||
X[:,i] = mumps.spsolve(self.A,b[:,i])
|
||||
|
||||
return X
|
||||
|
||||
def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
|
||||
if backend is None: backend = DEFAULTS['iter']
|
||||
|
||||
algorithms = {'CG':sp.linalg.cg, 'QMR':sp.linalg.qmr}
|
||||
assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
|
||||
alg = algorithms[iterSolver]
|
||||
|
||||
if iterSolver == 'CG':
|
||||
opts = {'M':M}
|
||||
elif iterSolver == 'QMR':
|
||||
#TODO: make preconditioner better.
|
||||
opts = {'M1':np.sqrt(M), 'M2':np.sqrt(M)}
|
||||
|
||||
if len(b.shape) == 1 or b.shape[1] == 1:
|
||||
x, self.info = alg(self.A, b, tol=tol, maxiter=maxIter)
|
||||
else:
|
||||
x = np.empty_like(b)
|
||||
for i in range(b.shape[1]):
|
||||
x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter)
|
||||
return x
|
||||
|
||||
def solveBackward(self, b, backend=None):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
Perform a backwards solve with upper triangular A in CSR format (best, if not, it will be converted).
|
||||
|
||||
:param str backend: which backend to use. Default is python.
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['triangular']
|
||||
if backend not in OPTIONS['triangular']:
|
||||
print 'Warning: %s-backend not being used, %s-default will be used instead.'%(backend,DEFAULTS['triangular'])
|
||||
backend = DEFAULTS['triangular']
|
||||
if type(self.A) is not sp.csr.csr_matrix:
|
||||
self.A = sp.csr_matrix(self.A)
|
||||
vals = self.A.data
|
||||
rowptr = self.A.indptr
|
||||
colind = self.A.indices
|
||||
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=None):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
Perform a forward solve with lower triangular A in CSR format (best, if not, it will be converted).
|
||||
|
||||
:param str backend: which backend to use. Default is python.
|
||||
:rtype: numpy.ndarray
|
||||
:return: x
|
||||
"""
|
||||
if backend is None: backend = DEFAULTS['triangular']
|
||||
if backend not in OPTIONS['triangular']:
|
||||
print 'Warning: %s-backend not being used, %s-default will be used instead.'%(backend,DEFAULTS['triangular'])
|
||||
backend = DEFAULTS['triangular']
|
||||
if type(self.A) is not sp.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
|
||||
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=None):
|
||||
"""
|
||||
Use solve instead of this interface.
|
||||
|
||||
Perform a diagonal solve with diagonal matrix A.
|
||||
|
||||
:param str backend: which backend to use. Default is python.
|
||||
: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
|
||||
return b/diagA
|
||||
# Multiple RHSs
|
||||
X = np.empty_like(b)
|
||||
for i in range(b.shape[1]):
|
||||
X[:,i] = b[:,i]/diagA
|
||||
return X
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG.Mesh import TensorMesh
|
||||
from time import time
|
||||
h1 = np.ones(20)*100.
|
||||
h2 = np.ones(20)*100.
|
||||
h3 = np.ones(20)*100.
|
||||
|
||||
h = [h1,h2,h3]
|
||||
|
||||
M = TensorMesh(h)
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
e = np.ones(M.nC)
|
||||
rhs = A.dot(e)
|
||||
|
||||
tic = time()
|
||||
solve = Solver(A, options={'factorize':True})
|
||||
x = solve.solve(rhs)
|
||||
print 'Factorized', time() - tic
|
||||
print np.linalg.norm(e-x,np.inf)
|
||||
tic = time()
|
||||
solve = Solver(A, options={'factorize':False})
|
||||
x = solve.solve(rhs)
|
||||
print 'spsolve', time() - tic
|
||||
print np.linalg.norm(e-x,np.inf)
|
||||
|
||||
|
||||
n = 600
|
||||
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 = sp.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
|
||||
|
||||
|
||||
|
||||
A = -D*D.T
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
e = np.ones(M.nC)
|
||||
b = A.dot(e)
|
||||
|
||||
iSolve = Solver(A, doDirect=False,options={'M':('GS',A)})
|
||||
tic = time()
|
||||
x = iSolve.solve(b)
|
||||
toc = time() - tic
|
||||
print x
|
||||
print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info
|
||||
|
||||
|
||||
|
||||
A = -D*D.T
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
e = np.ones(M.nC)
|
||||
b = A.dot(e)
|
||||
|
||||
iSolve = Solver(A, doDirect=False, options={'iterSolver': 'QMR', 'M':'J'})
|
||||
tic = time()
|
||||
x = iSolve.solve(b)
|
||||
toc = time() - tic
|
||||
print x
|
||||
print 'Error QMR = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info
|
||||
Solver = SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False)
|
||||
SolverLU = SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
|
||||
SolverCG = SolverUtils.ISolverWrap(sp.linalg.cg)
|
||||
|
||||
@@ -1,28 +0,0 @@
|
||||
from mumps import DMumpsContext
|
||||
|
||||
class MumpsSolver():
|
||||
A = None
|
||||
ctx = None
|
||||
x = None
|
||||
|
||||
def __init__(self, A, **kwagrs):
|
||||
|
||||
self.ctx = DMumpsContext(sym=0, par=1)
|
||||
|
||||
if self.ctx.myid ==0:
|
||||
self.A = A
|
||||
self.ctx.set_icntl(14, 60)
|
||||
self.ctx.set_centralized_sparse(A)
|
||||
|
||||
self.ctx.set_silent()
|
||||
self.ctx.run(job=4) # Factorization
|
||||
|
||||
def solve(self,b):
|
||||
self.x = b.copy()
|
||||
self.ctx.set_rhs(self.x)
|
||||
self.ctx.run(job=3) # Solve
|
||||
|
||||
return self.x
|
||||
|
||||
def clean(self):
|
||||
self.ctx.destroy()
|
||||
@@ -1,4 +0,0 @@
|
||||
try:
|
||||
from Mumps import MumpsSolver
|
||||
except Exception, e:
|
||||
pass
|
||||
+31
-126
@@ -1,5 +1,5 @@
|
||||
import unittest
|
||||
from SimPEG import Solver
|
||||
from SimPEG import *
|
||||
from SimPEG.Mesh import TensorMesh
|
||||
from SimPEG.Utils import sdiag
|
||||
import numpy as np
|
||||
@@ -8,136 +8,41 @@ import scipy.sparse as sparse
|
||||
TOL = 1e-10
|
||||
numRHS = 5
|
||||
|
||||
def dotest(solver, multi=False, **solverOpts):
|
||||
h1 = np.ones(10)*100.
|
||||
h2 = np.ones(10)*100.
|
||||
h3 = np.ones(10)*100.
|
||||
|
||||
h = [h1,h2,h3]
|
||||
|
||||
M = TensorMesh(h)
|
||||
|
||||
D = M.faceDiv
|
||||
G = -M.faceDiv.T
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
|
||||
Ainv = Solver(A, **solverOpts)
|
||||
if multi:
|
||||
e = np.ones(M.nC)
|
||||
else:
|
||||
e = np.ones((M.nC, numRHS))
|
||||
rhs = A * e
|
||||
x = Ainv * rhs
|
||||
return np.linalg.norm(e-x,np.inf) < TOL
|
||||
|
||||
class TestSolver(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
h1 = np.ones(10)*100.
|
||||
h2 = np.ones(10)*100.
|
||||
h3 = np.ones(10)*100.
|
||||
def test_direct_spsolve_1(self): self.assertTrue(dotest(Solver, False))
|
||||
def test_direct_spsolve_M(self): self.assertTrue(dotest(Solver, True))
|
||||
|
||||
h = [h1,h2,h3]
|
||||
def test_direct_splu_1(self): self.assertTrue(dotest(SolverLU, False))
|
||||
def test_direct_splu_M(self): self.assertTrue(dotest(SolverLU, True))
|
||||
|
||||
M = TensorMesh(h)
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
self.A = A
|
||||
self.M = M
|
||||
|
||||
def test_directFactored_1(self):
|
||||
solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
|
||||
e = np.ones(self.M.nC)
|
||||
rhs = self.A.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
|
||||
|
||||
def test_directFactored_M(self):
|
||||
solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'})
|
||||
e = np.ones((self.M.nC,numRHS))
|
||||
rhs = self.A.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
|
||||
def test_directSpsolve_1(self):
|
||||
solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'})
|
||||
e = np.ones(self.M.nC)
|
||||
rhs = self.A.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
|
||||
def test_directSpsolve_M(self):
|
||||
solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'})
|
||||
e = np.ones((self.M.nC, numRHS))
|
||||
rhs = self.A.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'})
|
||||
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)
|
||||
|
||||
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)
|
||||
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)
|
||||
rhs = AU.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
|
||||
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))
|
||||
rhs = AU.dot(e)
|
||||
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={})
|
||||
e = np.ones(self.M.nC)
|
||||
rhs = AD.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
|
||||
def test_directDiagonal_M(self):
|
||||
AD = sdiag(self.A.diagonal())
|
||||
solve = Solver(AD, doDirect=True, flag='D', options={})
|
||||
e = np.ones((self.M.nC,numRHS))
|
||||
rhs = AD.dot(e)
|
||||
x = solve.solve(rhs)
|
||||
self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
|
||||
def test_iterative_cg_1(self): self.assertTrue(dotest(SolverLU, False))
|
||||
def test_iterative_cg_M(self): self.assertTrue(dotest(SolverLU, True))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -3,7 +3,7 @@ import scipy.sparse as sp
|
||||
import unittest
|
||||
from TestUtils import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import Utils, Solver
|
||||
from SimPEG import *
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh']
|
||||
|
||||
@@ -54,7 +54,7 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
|
||||
err = np.linalg.norm((q-V*q_anal), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
solver = Solver(A, doDirect=False, options={'M':'J','iterSolver':'CG','backend':'scipy','maxIter':1000})
|
||||
solver = SolverCG(A, maxiter=1000)
|
||||
xc = solver.solve(rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
@@ -220,7 +220,7 @@ class Test1D_InhomogeneousNeumann(OrderTest):
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
@@ -228,7 +228,7 @@ class Test1D_InhomogeneousNeumann(OrderTest):
|
||||
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
@@ -268,12 +268,12 @@ class Test2D_InhomogeneousNeumann(OrderTest):
|
||||
j_anal = np.r_[jX_anal,jY_anal]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
|
||||
|
||||
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
|
||||
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
|
||||
gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
|
||||
gBCx = self.M.gridCC[(cxm|cxp),:]
|
||||
gBCy = self.M.gridCC[(cym|cyp),:]
|
||||
@@ -307,7 +307,7 @@ class Test2D_InhomogeneousNeumann(OrderTest):
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
@@ -315,7 +315,7 @@ class Test2D_InhomogeneousNeumann(OrderTest):
|
||||
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
@@ -384,7 +384,7 @@ class Test1D_InhomogeneousMixed(OrderTest):
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
@@ -392,7 +392,7 @@ class Test1D_InhomogeneousMixed(OrderTest):
|
||||
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
@@ -432,12 +432,12 @@ class Test2D_InhomogeneousMixed(OrderTest):
|
||||
j_anal = np.r_[jX_anal,jY_anal]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
|
||||
|
||||
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
|
||||
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
|
||||
gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
|
||||
gBCx = self.M.gridCC[(cxm|cxp),:]
|
||||
gBCy = self.M.gridCC[(cym|cyp),:]
|
||||
@@ -471,7 +471,7 @@ class Test2D_InhomogeneousMixed(OrderTest):
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
@@ -479,7 +479,7 @@ class Test2D_InhomogeneousMixed(OrderTest):
|
||||
err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
|
||||
@@ -6,7 +6,7 @@ import numpy as np
|
||||
class TestLinear(unittest.TestCase):
|
||||
|
||||
def test_running(self):
|
||||
Linear.run(100)
|
||||
Linear.run(100, plotIt=False)
|
||||
self.assertTrue(True)
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
+29
-18
@@ -2,10 +2,20 @@ import numpy as np
|
||||
from matutils import mkvc
|
||||
import warnings
|
||||
|
||||
def DSolverWrap(fun, factorize=True, destroy = False, checkAccuracy=True, accuracyTol=1e-6):
|
||||
def _checkAccuracy(A, b, X, accuracyTol):
|
||||
nrm = np.linalg.norm(mkvc(A*X - b), np.inf)
|
||||
nrm_b = np.linalg.norm(mkvc(b), np.inf)
|
||||
if nrm_b > 0:
|
||||
nrm /= nrm_b
|
||||
if nrm > accuracyTol:
|
||||
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
|
||||
print msg
|
||||
warnings.warn(msg, RuntimeWarning)
|
||||
|
||||
|
||||
def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
|
||||
|
||||
def __init__(self, A, **kwargs):
|
||||
|
||||
self.A = A.tocsc()
|
||||
self.kwargs = kwargs
|
||||
if factorize:
|
||||
@@ -28,27 +38,27 @@ def DSolverWrap(fun, factorize=True, destroy = False, checkAccuracy=True, accura
|
||||
X[:,i] = fun(self.A, b[:,i], **self.kwargs)
|
||||
|
||||
if checkAccuracy:
|
||||
nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b))
|
||||
if nrm > accuracyTol:
|
||||
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
|
||||
print msg
|
||||
warnings.warn(msg, RuntimeWarning)
|
||||
_checkAccuracy(self.A, b, X, accuracyTol)
|
||||
return X
|
||||
|
||||
def clean(self):
|
||||
if destroy == True:
|
||||
if hasattr(self.solver, 'clean'):
|
||||
return self.solver.clean()
|
||||
else:
|
||||
return True
|
||||
|
||||
return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean})
|
||||
|
||||
def __mul__(self, val):
|
||||
if type(val) is np.ndarray:
|
||||
return self.solve(val)
|
||||
raise TypeError('Can only multiply by a numpy array.')
|
||||
|
||||
return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "clean": clean, "__mul__": __mul__})
|
||||
|
||||
|
||||
|
||||
def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5):
|
||||
|
||||
def __init__(self, A, **kwargs):
|
||||
self.A = A.tocsc()
|
||||
self.A = A
|
||||
self.kwargs = kwargs
|
||||
|
||||
def solve(self, b):
|
||||
@@ -74,11 +84,12 @@ def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5):
|
||||
X[:,i] = out
|
||||
|
||||
if checkAccuracy:
|
||||
nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b))
|
||||
if nrm > accuracyTol:
|
||||
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
|
||||
print msg
|
||||
warnings.warn(msg, RuntimeWarning)
|
||||
_checkAccuracy(self.A, b, X, accuracyTol)
|
||||
return X
|
||||
|
||||
return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve})
|
||||
def __mul__(self, val):
|
||||
if type(val) is np.ndarray:
|
||||
return self.solve(val)
|
||||
raise TypeError('Can only multiply by a numpy array.')
|
||||
|
||||
return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve, "__mul__": __mul__})
|
||||
|
||||
@@ -1,64 +0,0 @@
|
||||
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
|
||||
+1
-1
@@ -1,7 +1,7 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import Utils
|
||||
from Solver import Solver
|
||||
from Solver import *
|
||||
import Mesh
|
||||
import Maps
|
||||
import Problem
|
||||
|
||||
Reference in New Issue
Block a user