diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py
index af384971..05ae09df 100644
--- a/SimPEG/tests/HTMLTestRunner.py
+++ b/SimPEG/tests/HTMLTestRunner.py
@@ -241,7 +241,8 @@ function showClassDetail(cid, count) {
for (var i = 0; i < count; i++) {
tid = id_list[i];
if (toHide) {
- document.getElementById('div_'+tid).style.display = 'none'
+ var divTid = document.getElementById('div_'+tid);
+ if(divTid !== null){divTid.style.display = 'none';}
document.getElementById(tid).className = 'hiddenRow';
}
else {
diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py
index 1cc2bbae..bfb150c3 100644
--- a/SimPEG/tests/TestUtils.py
+++ b/SimPEG/tests/TestUtils.py
@@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
-def getQuadratic(A, b):
+def getQuadratic(A, b, c=0):
"""
- Given A and b, this returns a quadratic, Q
+ Given A, b and c, this returns a quadratic, Q
.. math::
- \mathbf{Q( x ) = 0.5 x A x + b x}
+ \mathbf{Q( x ) = 0.5 x A x + b x} + c
"""
def Quadratic(x, return_g=True, return_H=True):
- f = 0.5 * x.dot( A.dot(x)) + b.dot( x )
+ f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c
out = (f,)
if return_g:
g = A.dot(x) + b
diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py
index a1194487..82c4b709 100644
--- a/SimPEG/utils/Solver.py
+++ b/SimPEG/utils/Solver.py
@@ -1,9 +1,10 @@
import numpy as np
-import scipy.sparse as sparse
+import scipy.sparse as sp
import scipy.sparse.linalg as linalg
-from SimPEG.utils import mkvc
+from SimPEG.utils import mkvc, sdiag
+import warnings
-DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
+DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
try:
import TriSolve
@@ -45,13 +46,41 @@ class Solver(object):
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 type(M) is sp.linalg.LinearOperator:
+ return
+ elif type(M) is tuple:
+ PreconditionerList = ['J','GS']
+ 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':
+ LL = sp.tril(M[1])
+ UU = sp.triu(M[1])
+ DD = sdiag(M[1].diagonal())
+ Uinv = Solver(UU, flag='U')
+ Linv = Solver(LL, 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):
"""
@@ -118,8 +147,20 @@ class Solver(object):
return X
- def solveIter(self, b, M=None, iterSolver='CG'):
- pass
+ 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}
+ assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
+ alg = algorithms[iterSolver]
+
+ if len(b.shape) == 1 or b.shape[1] == 1:
+ x, self.info = alg(self.A, b, M=M, 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):
"""
@@ -132,9 +173,8 @@ class Solver(object):
: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)
+ 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
@@ -164,7 +204,7 @@ class Solver(object):
:return: x
"""
if backend is None: backend = DEFAULTS['forward']
- if type(self.A) is not sparse.csr.csr_matrix:
+ 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
@@ -240,13 +280,13 @@ if __name__ == '__main__':
print np.linalg.norm(e-x,np.inf)
- n = 6000
+ 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 = sparse.csr_matrix(L)
+ A = sp.csr_matrix(L)
pSolve = Solver(A,flag='L',options={'backend':'python'});
fSolve = Solver(A,flag='L',options={'backend':'fortran'})
tic = time()
@@ -257,3 +297,17 @@ if __name__ == '__main__':
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