solver-eoc.py

#!/usr/bin/python

import sys
import getopt
import struct

from math import sqrt
from numpy import array, matrix
from scipy import sparse
from scipy.sparse import coo_matrix

####################################################
#  Print program usage
####################################################
def usage():
    print 'This program reads a matrix and an iteration history from disk'
    print 'and compares the iterates with the last one, the presumed exact solution.'
    print 'Errors in the energy norm and convergence rates are plotted.'
    print 'Options are:'
    print ' --help, -h:  show help'
    print ' --matrix, -m <file>         Matrix filename'
    print ' --iterates, -i <path>       Base name of the iterates.  A number will be appended'
    print ' --from, -f <arg>            Number of the first iterate (default: 0)'
    print ' --to, -t <arg>              Number of the last iterate'
    print ' --solution, -s <file>       The solution file (default: the last iterate)'
    
####################################################
#  Read a binary vector
####################################################
def read_binary_vector(filename, expected_size):

    vector = [None]*expected_size

    f = file(filename,'r')

    for i in range(0,expected_size):
        c = f.read(struct.calcsize('d'))
        if c=="":
            print 'Vector is too short!'
            exit(1)
        d = struct.unpack('d', c)[0]
        vector[i] = d
        #print '%0.20f' % d
        
    f.close()
    return matrix(vector[:],float)
    
####################################################
#  Main programm
####################################################

try:                                
    opts, args = getopt.getopt(sys.argv[1:], "hm:s:i:f:t:",
                               ["help", "matrix=", "solution:", "iterates=", "from=", "to="])
except getopt.GetoptError:
    usage()
    sys.exit(2)

# Default variables
fro = 0
solution_filename = None
iterates_filename = None

for opt, arg in opts:
    if opt in ("-h", "--help"):
        usage()                     
        sys.exit()                  
    elif opt in ("-m", "--matrix"):
        matrix_filename = arg    
    elif opt in ("-s", "--solution"):
        solution_filename = arg    
    elif opt in ("-i", "--iterates"):
        iterates_filename = arg    
    elif opt in ("-f", "--from"):
        fro = int(arg)    
    elif opt in ("-t", "--to"):
        to = int(arg) + 1

# Error if no filename for the iterates is given
if iterates_filename = None:
    print "Please provide a name for the iterates (option '-i')"
    usage()
    sys.exit(2)

# If no explicit solution file has been given we use the last one
# in the list of iterates
if solution_filename == None:
    solution_filename = iterates_filename + '%04d' %(to-1)

print 'Matrix filename:', matrix_filename
print 'Solution filename:', solution_filename

####################################################
#  Read a matrix.
#  This is not in a separate function because I didn't get
#  call-by-reference for the matrix to work yet.
####################################################

f = open(matrix_filename)

counter = 0;
for line in f:
    if (counter % 100000)==0:
        print counter
    counter += 1

f.close()
print 'matrix file consists of', counter, 'lines'

rows = array([0]*counter,int)
cols = array([0]*counter,int)
values = array([0]*counter,float)
counter = 0;
f = open(matrix_filename)
for line in f:
    line_entries = line.split()
    rows[counter] = int(line_entries[0]) - 1
    cols[counter] = int(line_entries[1]) - 1
    values[counter] = float(line_entries[2])
    if (counter % 100000)==0:
        print counter
    counter += 1
        
print 'Reading a', max(rows)+1, 'x', max(cols)+1, 'matrix'

matrix_size = max(rows)+1
A = coo_matrix((values,(rows,cols)), dims=(max(rows)+1,max(cols)+1))

# The vector taken as the 'exact' solution
solution = read_binary_vector(solution_filename,matrix_size)

# Main loop
conv_rate         = None
conv_rate_product = 1.
total_conv_rate   = None

for i in range(fro,to):

    # Compute the error vector
    leading_zeros =  ''.join(['0']*(4-len(str(i))))
    
    current_iterate = read_binary_vector(iterates_filename + '%04d' %i, matrix_size)
    current_iterate -= solution

    # Compute the energy norm of the error
    error = sqrt(current_iterate*A*current_iterate.transpose())

    # Compute the convergence rate
    if (i!=fro):
        conv_rate = error / old_error
        conv_rate_product *= conv_rate
        
    old_error = error
    if (i!=fro):
        total_conv_rate = pow(conv_rate_product, 1/float(i))

    # Output
    print 'Iteration:', i, 'error:', error, '\t rate:', conv_rate, '\t total rate:', total_conv_rate