Complete SOR Code

implementSOR.py

@@ -0,0 +1,118 @@
+# -*- coding: utf-8 -*-
+"""
+
+
+
+"""
+import numpy as np
+
+def get_euclid_norm(matrix):
+    """
+    this function accepts one vector
+    the calling function should compute this vector as the difference of two vectors
+    the standard Euclidean distance formula sqrt(x1^2+x2^2+....+xn^2) is applied
+    """
+    sum = 0
+    for each_val in matrix:
+        sum+=each_val**2
+    return np.sqrt(sum)
+
+def get_residual(A, x, b):
+    """
+    return the residual error Ax-b for an approximation of x
+    input parameters should be numpy arrays
+    this is not as simple as using numpy.dot because A is in CSR format 
+    """
+    adotx = []    
+    for i in range(0,len(b)):
+        adotx.append(0.0)    
+    #i should really do a DOT function instead of coding it explicitly here and in SOR also
+    for j in range (0, len(b)): 
+        first_nz = A[1][(A[2][j]-1)] - 1 #pos of first non zero on row             
+        for k in range(A[2][j]-1, A[2][j+1]-1):            
+            adotx[j] = adotx[j] + A[0][k] * x[k - (A[2][j]-1) + first_nz]
+    return get_euclid_norm(np.subtract(adotx, b))    
+
+def get_x_seq(xold, xnew):
+    """
+    this function computes Euclidean distance between successive iterations of x
+    input parameters should be numpy arrays
+    """
+    return get_euclid_norm(np.subtract(xnew, xold))
+
+def chk_diverge(xold, xnew, A, b):
+    """
+    check if previous approx of x was closer than new approx of x
+    """
+    dist_old = get_residual(A, xold, b)
+    dist_new = get_residual(A, xnew, b)
+    if dist_old < dist_new:
+        return True
+    else:
+        return False   
+
+def chk_converge(A, xnew, b, xold, x_seq_tol, res_tol, flag):
+    #checks both residual and x_seq for convergence 
+    if flag == True:
+        return -1 #required to enter sparse_sor loop
+    elif get_residual(A, xnew, b) < res_tol:
+        return 2 #dict value for this stopping reason
+    elif get_x_seq(xold, xnew) < x_seq_tol:
+        return 1
+    elif chk_diverge(xold, xnew, A, b) == True:
+        return 4
+    return -1
+
+
+
+def sparse_sor(matrix_a, matrix_b, matrix_x, dimension_n, max_it, \
+        x_tol, res_tol, w ):
+
+    num_it = 1
+    stop_reason = 6 #something has gone wrong if this does not get overwritten later
+    matrix_x_last_it = np.array([0.0])
+    matrix_x_new = np.array(matrix_x)
+    matrix_a_np = np.array(matrix_a)
+    matrix_b_np = np.array(matrix_b)
+
+    flag = True #required to enter while loop first time only
+
+    while num_it <= max_it and \
+        chk_converge(matrix_a_np, matrix_x_new, matrix_b_np, matrix_x_last_it, 
+                     x_tol, res_tol, flag) == -1:
+
+        flag = False        
+        matrix_x_last_it = np.array(matrix_x_new[:])
+        for i in range(0,len(matrix_b)):
+            sum = 0
+            first_nz = matrix_a_np[1][(matrix_a_np[2][i]-1)] - 1 #pos of first non zero on row
+            for j in range(matrix_a_np[2][i]-1, matrix_a_np[2][i+1]-1):
+                sum = sum + matrix_a_np[0][j] * matrix_x_new[j - \
+                            (matrix_a_np[2][i]-1) + first_nz]
+                if matrix_a_np[1][j] == i+1:
+                    d = matrix_a[0][j]
+            matrix_x_new[i] = matrix_x_new[i] + w * (matrix_b_np[i] - sum) / d   
+        num_it+=1
+    conv = chk_converge(matrix_a_np, matrix_x_new, matrix_b_np, matrix_x_last_it, \
+                    x_tol, res_tol, False)
+    if num_it-1 == max_it:
+        stop_reason = 3
+    elif conv != -1:
+        stop_reason = conv
+#    processIO.output(err.args[1], max_it, num_it, x_tol, res_tol, matrix_x, output_filename,err.args[0])           
+#    set_output(stop_reason, max_it, num_it-1, x_tol, res_tol,matrix_x_new, )
+    return (stop_reason, max_it, num_it-1, matrix_x_new )
+
+def dense_SOR(matrix_a,matrix_b,dimension_n,max_iteration,w,matrix_x):
+    print("Dense SOR Calculation")
+    it_counter=0
+    while(it_counter<=max_iteration):
+        for row_counter in range(dimension_n):
+            sum=0.0
+            for col_counter in range(dimension_n):
+                sum = sum + matrix_a[row_counter,col_counter]*matrix_x[col_counter]
+            matrix_x[row_counter]=matrix_x[row_counter]+w*(matrix_b[row_counter]-sum) \
+                                                /matrix_a[row_counter,row_counter]
+#        print("Iteration: ",it_counter,"\n","X:",matrix_x)    
+        it_counter+=1
+    return

mainSOR.py

@@ -0,0 +1,75 @@
+# -*- coding: utf-8 -*-
+"""
+Coding Practices:
+All variables and functions are named in lower-case letters seperated by 
+underscores between words. e.g. `max_n_users`. 
+All Classes are usually named with CamelCase. e.g. `DatabaseConnection`.
+
+"""
+
+
+
+#import modules
+import sys
+import processIO
+import implementSOR
+import validateMatrix
+
+
+def main():
+
+    try: 
+        #Input Processing from input file
+        if(len(sys.argv) >= 2):
+            [dimension_n,matrix_a,matrix_b] = processIO.input(sys.argv[1])
+        else:
+            [dimension_n,matrix_a,matrix_b]=processIO.input()        
+
+        #Constants
+        max_it = 50
+        num_it = 0
+        w=1.25
+        x_tol=1e-13
+        res_tol=1e-13
+        matrix_x = processIO.np.zeros(dimension_n)
+
+
+        #Output processing
+        if(len(sys.argv) >=3):
+            output_filename = sys.argv[2]
+        else:
+            output_filename = "nas_Sor.out" #="nas_Sor.out"
+
+        #checks on A and B Matrix
+        if(not validateMatrix.non_zero_diagonal_check(matrix_a,dimension_n)):
+            raise Exception("Determinant check failed for Matrix A.",5)
+
+        if(not validateMatrix.det_check(matrix_a)):
+            raise Exception("Determinant check failed for Matrix A.",6)
+
+        if((not validateMatrix.col_diagonally_dominant(matrix_a,dimension_n)) and \
+            (not validateMatrix.row_diagonally_dominant(matrix_a,dimension_n))):
+            raise Exception("Matrix A wont converge as it is not row / column diagonally dominant.",6)
+
+        #Dense SOR Calculation
+#       implementSOR.dense_SOR(matrix_a,matrix_b,dimension_n,max_it,w,matrix_x)
+
+        #Sparse SOR Calculation
+        csr_matrix_a = processIO.to_csr_format(matrix_a)
+
+        [stop_reason, max_it, num_it, matrix_x ] = implementSOR.sparse_sor(csr_matrix_a, \
+        matrix_b, matrix_x, dimension_n, max_it, x_tol, res_tol, w )
+
+        processIO.output(stop_reason, max_it, num_it, x_tol, res_tol, matrix_x,\
+                output_filename)           
+
+    except Exception as err:
+#        print("Error occured: ",err.args[0],err.args[1])
+        processIO.output(err.args[1], max_it, num_it, x_tol, res_tol, matrix_x, \
+            output_filename,err.args[0])
+
+    return
+
+
+if __name__ == "__main__":
+    main()    

nas_Sor.in

@@ -0,0 +1,5 @@
+3
+3 -1 1
+-1 3 -1
+1 -1 3
+-1 7 -7

nas_Sor.in.err

@@ -0,0 +1,7 @@
+5
+4 1 2 -3 5
+-3 3 -1 4 -2
+-1 2 5 1 3
+5 4 3 -1 2
+1 -2 3 -4 5
+-16 20 -4 -10 3

nas_Sor_incorrect_format.in

@@ -0,0 +1,6 @@
+3
+3 -1 1
+-1 0 -1
+1 -1 3
+1 -1 3  
+-1 7 -7

processIO.py

@@ -0,0 +1,111 @@
+# -*- coding: utf-8 -*-
+"""
+
+"""
+
+import numpy as np
+import sys
+
+def input(file_name="nas_Sor.in"):
+    try:
+#        print("FIle Name:",file_name)
+        dimension_n = int(np.genfromtxt(file_name,max_rows=1))
+#        print("Dimension n :",dimension_n)
+
+        matrix_a = np.genfromtxt(file_name,skip_header=1,skip_footer=1)
+#        print("Matrix A :",matrix_a, "Shape:",matrix_a.shape)        
+        if(matrix_a.shape[0]!=dimension_n or matrix_a.shape[1]!=dimension_n):
+            raise Exception('Incorrect File Format for A Matrix dimensions')
+
+        matrix_b=np.genfromtxt(file_name,skip_header=dimension_n+1)
+#        print("Matrix B :",matrix_b, "Shape:",matrix_b.shape)
+        if(matrix_b.shape[0]!=dimension_n):
+            raise Exception('Incorrect File Format for B Matrix Dimensions')
+
+    except OSError:
+        print("Input File : ",file_name," not found. ")
+        sys.exit()
+
+    except Exception as err:
+        print(" Input File : ",file_name,"\n","Error: ",err.args[0])
+        print(" Please check format of the file. ","\n","Accepted format :","\n\t" \
+            " First Line should contain dimension of matrix - n","\n\t", \
+            "Next n lines should contain matrix A - n * n ","\n\t",\
+            "Last Line should contain matrix B with n rows")
+        sys.exit()
+
+    return (dimension_n,matrix_a,matrix_b)
+
+def output(stop_reason, max_it, num_it, x_tol, res_tol, matrix_x, file_name, \
+            stop_cause='something has gone wrong'):
+    try:
+        stop_reasons = {
+        1 : "x Sequence convergence"
+        ,2 : "Residual convergence"
+        ,3 : "Max Iterations reached"
+        ,4 : "X Sequence divergence"
+        ,5 : "Zero on diagonal"
+        ,6 : "Cannot proceed"
+        }
+
+        print("Output File - ",file_name)
+
+#        print("Stopping Reason {20} ","Max num of Iterations","\t","Number of Iterations", \
+#        "\t","Machine Epsilon","\t","X Seq Tolerance","\t","Residual Seq Tolerance")
+
+        header_line = "{:<30} {:<30} {:<30} {:<30} {:<30} {:<30}".format("Stopping Reason", \
+                    "Max num of Iterations","Number of Iterations", \
+                    "Machine Epsilon","X Seq Tolerance","Residual Seq Tolerance")
+
+        values_line = "{:<30} {:<30} {:<30} {:<30} {:<30} {:<30}".format(stop_reasons[stop_reason] \
+                , max_it, num_it, np.finfo(float).eps, x_tol, res_tol)
+
+        print(header_line)
+        print(values_line)
+
+        if stop_reason in (1,2,3):
+            print(matrix_x)
+
+        if stop_reason ==6:
+            print(stop_cause)
+
+    except:
+        print("Output File - ",file_name," processing error. ")
+    return
+
+def to_csr_format(matrix_a):
+    """
+    to_csr_format takes a matrix and converts it into dense form
+    row1: values, row2: columns, row3: rowstart
+
+    takes a full line-by-line matrix and compresses into non-zero 'rowstart' CSR format
+    row 1 - values, row 2 - columns, row 3 - rowstart
+    """
+    csr_format_a = []
+    cols = []
+    vals = []
+    rowstart = []
+
+    row_num = 1
+    for each_row in matrix_a:
+        col_num = 1
+        rowstarted = False
+
+        for val in each_row:            
+#            if row_num == col_num and matrix_a[row_num - 1][col_num - 1] == 0:
+#                    return [5] #STOP - zero on diagonal
+            if val != 0:
+                vals.append(val)
+                cols.append(col_num)
+                if rowstarted == False:
+                    rowstart.append(len(vals))
+                    rowstarted = True
+            col_num+=1
+        row_num+=1
+    rowstart.append(len(vals)+1) #final value
+
+    csr_format_a.append(vals)
+    csr_format_a.append(cols)
+    csr_format_a.append(rowstart)  
+
+    return csr_format_a

validateMatrix.py

@@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+"""
+
+"""
+import numpy as np
+
+def non_zero_diagonal_check(matrix_a,dimension_n):
+    for counter in range(dimension_n):
+        if(matrix_a[counter,counter]==0.0):
+            return False
+    return True
+
+
+def det_check(matrix_a):
+    if(np.linalg.det(matrix_a)==0.0):
+        return False
+    return True
+
+def row_diagonally_dominant(matrix_a,dimension_n):
+    for row_counter in range(dimension_n):
+        sum = 0
+        for col_counter in range(dimension_n):
+            if(row_counter != col_counter):
+                sum += abs(matrix_a[row_counter,col_counter])
+#        print("Row Sum: ",sum,"Diag Value:",matrix_a[row_counter,row_counter])
+        if (abs(matrix_a[row_counter,row_counter])> sum):
+            return True
+    return False
+
+def col_diagonally_dominant(matrix_a,dimension_n):
+    for col_counter in range(dimension_n):
+        sum = 0
+        for row_counter in range(dimension_n):
+            if(row_counter != col_counter):
+                sum += abs(matrix_a[row_counter,col_counter])
+#        print("Col Sum: ",sum,"Diag Value:",matrix_a[col_counter,col_counter])
+        if (abs(matrix_a[col_counter,col_counter])> sum):
+            return True
+    return False