Multiple RHSs on solvers in Fortran. ~2x speed up on matlab implementation for a single RHS. for multiple RHS there are still some problems.

Someone with some knowledge of how fortran works should look at this code.

Added a setup.py script that complies things. f2py should work on most computers, because it is included in the numpy distribution.

SimPEG/utils/Solver.py

@@ -1,14 +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:
-    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
+    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):
     """
@@ -76,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.
 
@@ -85,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:
@@ -111,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.
 
@@ -121,6 +131,7 @@ 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)
@@ -128,7 +139,11 @@ class Solver(object):
         rowptr = self.A.indptr
         colind = self.A.indices
         if backend == 'fortran':
-            x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size)
+            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])):
@@ -138,7 +153,7 @@ class Solver(object):
                 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.
 
@@ -148,6 +163,7 @@ 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)
@@ -155,7 +171,11 @@ class Solver(object):
         rowptr = self.A.indptr
         colind = self.A.indices
         if backend == 'fortran':
-            x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size)
+            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]):
@@ -165,7 +185,7 @@ class Solver(object):
                 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.
 
@@ -175,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