Matrix solvers

From GeoMod

This file, which I call matrix_solvers.py, contains matrix solvers using the following methods;

• jacobi - the Point Jacobi method
• gauss-seidel - the Gauss-Seidel method (usually converges twice as fast as the Point Jacobi method.
• sor - Successive over-relaxation
  • This method accelerates the Gauss-Seidel method using a relaxation factor (rf > 1.0)

Descriptions and details regarding these methods can be found at

• Lecture notes of James Demmel from UC Berkeley
• Kartik Patel's website on "Relaxation methods for solving Partial Differential Equations"

Contents 1 Discussion 1.1 Efficiency 2 The file: matrix_solvers.py 3 Example of use: test_matrix_solvers.py

Discussion

Efficiency

As programmed the Jacobi solver is actually fastest even though it requires more itterations to solve the matrix. Running a benchmark simulation with a banded 150x150 matrix found that;

Efficiency (150,150) matrix

Solver	Number of iterations	Run Time (seconds)
jacobi	5510	3.8
gauss-seidel	2701	10.0
sor (rf=1.1)	2703	8.2
sor (rf=1.5)	2038	6.1
sor (rf=1.9)	1648	5.0

Efficiency calculations (2000, 2000) matrix

Solver	Number of iterations	Run Time (seconds)
jacobi	8355	909
gauss-seidel
sor (rf=1.1)
sor (rf=1.5)	4137	812
sor (rf=1.9)

The file: matrix_solvers.py

from numpy import *
import time

def jacobi(A, B, tol=0.005, init_val=0.0):
    '''itteratively solving for matrix A with solution vector B
       tol = tolerance for dh/h
       init_val = array of initial values to use in the solver
    '''
    print "jacobi solver"
    if init_val == 0.0:
        h = zeros(shape(B), float)
    else:
        h = init_val
    h = zeros(shape(B), float)
    dmax = 1.0
    n = 0

    t1 = t = time.time()
    while dmax > tol:
        n += 1
        #print "  Iteration = ", n
        hnew = h+ divide(subtract(B, add.reduce(multiply(A,h),1)), diag(A))
##        print_arr(multiply(A,h), "multiply(A,h)")
##        print_arr(add.reduce(multiply(A,h),1), "add.reduce(multiply(A,h))")
##        print_arr(subtract(B, add.reduce(multiply(A,h),1)), "subtract(B, add.reduce(multiply(A,h)))")
##        print_arr(hnew, "hnew")
        dmax = max(abs(divide(subtract(h,hnew), h)))
        #print "   dmax = ", dmax, ": tol = ", tol
        #print "subtract", subtract(h,hnew)
        if (time.time() - t) > 2.0:
            print "Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
            t = time.time()
        h = hnew
##        print n, dmax,  h
##        print " "
    #print_arr(h, "h")
    print "SOLVED at: Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
    print "                           Time = ", time.time() - t1, "seconds"

    return h


def gauss_seidel(A, B, tol=0.005, init_val=0.0):
    '''itteratively solving for matrix A with solution vector B
       tol = tolerance for dh/h
       init_val = array of initial values to use in the solver
    '''
    print "gauss_siedel"
    if init_val == 0.0:
        h = zeros(shape(B), float)
        hnew = zeros(shape(B), float)
    else:
        h = init_val
        hnew[:] = h[:]
    dmax = 1.0
    n = 0

    t1 = t = time.time()
    while dmax > tol:
        n += 1
        #print "  Iteration = ", n
        for i in range(len(B)):
            #print "add.reduce(multiply(A,h)) ", add.reduce(multiply(A,h))
            hnew[i] = h[i] + (B[i] - add.reduce(multiply(A[i,:],hnew)))/ A[i,i]
            #divide(add.reduce(subtract(B, multiply(hnew, A[i,:]))), A[i,i])
            #( B - h * A[i,:]) / A[i,i]
        dmax = max(abs(divide(subtract(h,hnew), h)))
        #print "Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
##        print_arr(h, "h")
##        print_arr(hnew, "hnew")
        if (time.time() - t) > 2.0:
            print "Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
            t = time.time()
        h[:] = hnew[:]
##        print n, ":", dmax,  h   
    print "SOLVED at: Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
    print "                           Time = ", time.time() - t1, "seconds"
    
    return h

def sor(A, B, rf, tol=0.005, init_val=0.0):
    '''itteratively solving for matrix A with solution vector B
       rf = relaxation factor (rf > 1.0; eg rf = 1.1)
       tol = tolerance for dh/h
       init_val = array of initial values to use in the solver
       
    '''
    print "successive over-relaxation (sor)"
    if init_val == 0.0:
        h = zeros(shape(B), float)
        hnew = zeros(shape(B), float)
    else:
        h = init_val
        hnew[:] = h[:]
    dmax = 1.0
    n = 0

    t1 = t = time.time()
    while dmax > tol:
        n += 1
        for i in range(len(B)):
            #print "add.reduce(multiply(A,h)) ", add.reduce(multiply(A,h))
            hnew[i] = h[i] + (B[i] - add.reduce(multiply(A[i,:],hnew)))/ A[i,i]

        #over-relaxation
        dh = subtract(hnew, h)
        hnew = h + multiply(rf,dh)
        dmax = max(abs(divide(subtract(h,hnew), h)))
        if (time.time() - t) > 2.0:
            print "Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
            t = time.time()
        
##        print_arr(h, "h")
##        print_arr(hnew, "hnew")
        h[:] = hnew[:]
    print "SOLVED at: Itteration = ", n, ":   dmax = ", dmax, ": tol = ", tol
    print "                           Time = ", time.time() - t1, "seconds"
    
    return h

Example of use: test_matrix_solvers.py

This bit of code demonstrates how to use the matrix solvers.

from numpy import *
from matrix_solvers import *

def print_arr(arr, title):
    print title, shape(arr)
    print arr

A = array( (4, -1, -1, 0, -1, 4, 0, -1, -1, 0, 4, -1, 0, -1, -1, 4),float)
A = reshape(A, (4,4))
print_arr( A, "A")

B = array( (1,2,0,1), float)
print_arr( B, "B")

print "---"

C = jacobi(A,B)
print_arr(C, "C jacobi")

C = gauss_seidel(A,B)
print_arr(C, "C gauss_seidel")

C = sor(A, B, rf=1.3)
print_arr(C, "C sor")