From f3170a8a7f4db3b308920dba01becfcb5ba36131 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 16 May 2014 18:33:57 -0700 Subject: [PATCH] #69 Removed most solvers, replaced by wrappers. --- SimPEG/Examples/Linear.py | 19 +- SimPEG/Optimization.py | 4 +- SimPEG/Solver.py | 407 +-------------------------- SimPEG/Solvers/Mumps.py | 28 -- SimPEG/Solvers/__init__.py | 4 - SimPEG/Tests/test_Solver.py | 157 ++--------- SimPEG/Tests/test_boundaryPoisson.py | 28 +- SimPEG/Tests/test_example_linear.py | 2 +- SimPEG/Utils/SolverUtils.py | 47 ++-- SimPEG/Utils/TriSolve.f | 64 ----- SimPEG/__init__.py | 2 +- 11 files changed, 92 insertions(+), 670 deletions(-) delete mode 100644 SimPEG/Solvers/Mumps.py delete mode 100644 SimPEG/Solvers/__init__.py delete mode 100644 SimPEG/Utils/TriSolve.f diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Linear.py index d56cc3e7..931af126 100644 --- a/SimPEG/Examples/Linear.py +++ b/SimPEG/Examples/Linear.py @@ -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 diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 7e9c4e37..12dc385e 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -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 diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index e91affb9..8f96923c 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -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) diff --git a/SimPEG/Solvers/Mumps.py b/SimPEG/Solvers/Mumps.py deleted file mode 100644 index 251ccbdb..00000000 --- a/SimPEG/Solvers/Mumps.py +++ /dev/null @@ -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() diff --git a/SimPEG/Solvers/__init__.py b/SimPEG/Solvers/__init__.py deleted file mode 100644 index f49b54d5..00000000 --- a/SimPEG/Solvers/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -try: - from Mumps import MumpsSolver -except Exception, e: - pass diff --git a/SimPEG/Tests/test_Solver.py b/SimPEG/Tests/test_Solver.py index 76177990..9135e851 100644 --- a/SimPEG/Tests/test_Solver.py +++ b/SimPEG/Tests/test_Solver.py @@ -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__': diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index 5e028dee..d4ac504b 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -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): diff --git a/SimPEG/Tests/test_example_linear.py b/SimPEG/Tests/test_example_linear.py index 6b271939..b42d8a70 100644 --- a/SimPEG/Tests/test_example_linear.py +++ b/SimPEG/Tests/test_example_linear.py @@ -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__': diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 282ea802..af721d82 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -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__}) diff --git a/SimPEG/Utils/TriSolve.f b/SimPEG/Utils/TriSolve.f deleted file mode 100644 index ab9048d6..00000000 --- a/SimPEG/Utils/TriSolve.f +++ /dev/null @@ -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 diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index ea317e3c..299bd01e 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -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