From fff2375acd1f2b58a334a77416b54ce434b9949e Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 21 Nov 2013 14:34:32 -0800 Subject: [PATCH 1/3] Moved scipy direct code to own function. --- SimPEG/utils/Solver.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 16db872f..85d5f640 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -120,6 +120,7 @@ class Solver(object): """ 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 @@ -129,6 +130,20 @@ class Solver(object): assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' + if backend == 'scipy': + X = self.solveDirect_scipy(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) From 623ff674ea1da5f97cf62453bc3a32411e2bfa8d Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 21 Nov 2013 15:00:51 -0800 Subject: [PATCH 2/3] Added pymumps direct solver backend. --- SimPEG/utils/Solver.py | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 85d5f640..c451a7ca 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -19,6 +19,11 @@ except Exception, e: DEFAULTS['forward'] = 'python' DEFAULTS['backward'] = 'python' +try: + import mumps +except Exception, e: + print 'Warning: mumps solver not available.' + class Solver(object): """ Solver is a light wrapper on the various types of @@ -113,6 +118,9 @@ class Solver(object): 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 @@ -132,6 +140,8 @@ class Solver(object): if backend == 'scipy': X = self.solveDirect_scipy(b, factorize) + elif backend == 'mumps': + X = self.solveDirect_mumps(b, factorize) return X @@ -165,6 +175,40 @@ class Solver(object): 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: + return self.dsolve(b) + else: + return mumps.spsolve(self.A, b) + + else: + raise NotImplementedError('Multiple RHS not yet implemented with mumps solver.') + def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50): if backend is None: backend = DEFAULTS['iter'] From 1439052172dc2048494feadf5a63d4c8cb032a3a Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 22 Nov 2013 15:39:06 -0800 Subject: [PATCH 3/3] Multiple RHS using mumps backend. --- SimPEG/utils/Solver.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index c451a7ca..5861c146 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -187,7 +187,7 @@ class Solver(object): if factorize and self.dsolve is None: self.mctx = mumps.DMumpsContext() self.mctx.set_icntl(14, 60) - self.mctx.set_silent() + # self.mctx.set_silent() self.mctx.set_centralized_sparse(self.A) self.mctx.run(job=4) @@ -202,12 +202,20 @@ class Solver(object): if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS if factorize: - return self.dsolve(b) + X = self.dsolve(b) else: - return mumps.spsolve(self.A, b) + X = mumps.spsolve(self.A, b) else: - raise NotImplementedError('Multiple RHS not yet implemented with mumps solver.') + # 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']