diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py
index 9b5cc0e7..bfa57761 100644
--- a/SimPEG/tests/test_Solver.py
+++ b/SimPEG/tests/test_Solver.py
@@ -58,17 +58,33 @@ class TestSolver(unittest.TestCase):
         x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
-    def test_directLower_1(self):
+    def test_directLower_1_fortran(self):
         AL = sparse.tril(self.A)
-        solve = Solver(AL, doDirect=True, flag='L', options={})
+        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(self):
+    # 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)
+
+    def test_directLower_1_python(self):
         AL = sparse.tril(self.A)
-        solve = Solver(AL, doDirect=True, flag='L', options={})
+        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)
@@ -90,6 +106,23 @@ class TestSolver(unittest.TestCase):
         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={})
diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py
index d9ac4a1d..1cc2c743 100644
--- a/SimPEG/utils/Solver.py
+++ b/SimPEG/utils/Solver.py
@@ -2,6 +2,13 @@ import numpy as np
 import scipy.sparse as sparse
 import scipy.sparse.linalg as linalg
 
+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
 
 class Solver(object):
     """
@@ -55,11 +62,11 @@ class Solver(object):
         elif self.flag is None and not self.doDirect:
             return self.solveIter(b, **self.options)
         elif self.flag == 'U':
-            return self.solveBackward(b)
+            return self.solveBackward(b, **self.options)
         elif self.flag == 'L':
-            return self.solveForward(b)
+            return self.solveForward(b, **self.options)
         elif self.flag == 'D':
-            return self.solveDiagonal(b)
+            return self.solveDiagonal(b, **self.options)
         else:
             raise Exception('Unknown flag.')
         pass
@@ -120,12 +127,15 @@ class Solver(object):
         vals = self.A.data
         rowptr = self.A.indptr
         colind = self.A.indices
-        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]
+        if backend == 'fortran':
+            x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size)
+        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='python'):
@@ -144,12 +154,15 @@ class Solver(object):
         vals = self.A.data
         rowptr = self.A.indptr
         colind = self.A.indices
-        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]
+        if backend == 'fortran':
+            x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size)
+        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='python'):
@@ -205,3 +218,20 @@ if __name__ == '__main__':
     print np.linalg.norm(e-x,np.inf)
 
 
+    n = 6000
+    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)
+    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
diff --git a/SimPEG/utils/TriSolve.f b/SimPEG/utils/TriSolve.f
new file mode 100644
index 00000000..6d15d013
--- /dev/null
+++ b/SimPEG/utils/TriSolve.f
@@ -0,0 +1,54 @@
+c File TriSolve.f
+       subroutine forward(al, ial, jal, b, nv, n,  x)
+       double precision al(nv)
+       integer ial(n+1)
+       integer jal(nv)
+       double precision b(n)
+       double precision x(n)
+       integer nv
+       integer n
+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(out) :: x
+       real ( kind = 8 ) t
+
+       do k = 1, n
+          t = b(k)
+          do j = ial(k)+1, ial(k+1)
+             t = t - al(j) * x(jal(j)+1)
+          end do
+          x(k) = t/al(ial(k+1))
+       end do
+       end subroutine forward
+
+
+       subroutine backward(au,iau, jau, b, nv, n, x)
+       double precision au(nv)
+       integer iau(n+1)
+       integer jau(nv)
+       double precision b(n)
+       double precision x(n)
+       integer nv
+       integer n
+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(out) :: x
+       real ( kind = 8 ) t
+
+       do k = n, 1, -1
+          t = b(k)
+          do j = iau(k)+1, iau(k+1)
+             t = t - au(j) * x(jau(j)+1)
+          end do
+          x(k) = t/au(iau(k)+1)
+       end do
+
+       end subroutine backward