python gauss.py < steuer

# simple gauss elimination
# Written by Rainer Gerhards 2009-05-26 * <rgerhards@adiscon.com>
# WARNING: no syntax/semantic checks done so far on input stream!
# Needs the numpy package for matrix operations
#

import sys, operator
import numpy
#import numpy.linalg


def getNewWords():
   """get list of words from a single line of stdin"""
   while True:
      line = sys.stdin.readline()
      if line == "":
         return
      if line[0] == '#' or line == "\n":
         continue
      # a bit ugly, but works...
      mylist = []
      for item in line.strip().split():
         if item != '':
            mylist.append(item)
      return mylist

def initMat(m, rows, cols):
   m['rows'] = int(rows)
   m['cols'] = int(cols)
   m['mat']  = numpy.zeros((rows, cols), numpy.float_, 'C')


def fillDataMat(m):
   words = getNewWords()
   for i in range(m['rows']):
      if words[0] == "enddata":
         return
      for j in range(m['cols']):
         m['mat'][i, j] = float(words[j])
      words = getNewWords()


def matCreate(name):
   """create a new matrix object"""
   m = {}
   m['name'] = name
   words = getNewWords()
   while words:
      if words[0] == "dim":
         initMat(m, int(words[1]), int(words[2]))
      elif words[0] == "data":
         fillDataMat(m)
      elif words[0] == "endmatrix":
         return m
      else:
         print "unexpected token in input: '" + words[0] + "'"
      words = getNewWords()


def matPrint(M, name):
   """print matrix """
   print "Matrix " + name + "" + str(M.shape) + ":"
   print numpy.array_str(M)
   print


def matPrintAll():
   for name in allMat:
      matPrint(allMat[name], name)


def matDup(name, newname):
   """ duplicate matrix """
   new = {}
   new['name'] = newname
   allMat[newname] = new
   mat = allMat[name]
   M = mat['mat'].copy()
   new['mat'] = M
   new['rows'] = mat['rows']
   new['cols'] = mat['cols']


# helpers mentioned, but not provided, in 1142 script
def swaprows(M, i, k):
   if verbose:
      print "swap rows " + str(i) + "," + str(k)
   for j in range(M.shape[1]):
      x = M[i, j]
      M[i, j] = M[k, j]
      M[k, j] = x

def swapcolumns(M, index, i, k):
   """ swap columns in a matrix """
   """ Note that we must also swap variable indexes, as we rename them """
   if verbose:
      print "swap columns " + str(i) + "," + str(k)
   x = index[k]
   index[k] = index[i]
   index[i] = x
   for j in range(M.shape[0]):
      x = M[j, i]
      M[j, i] = M[j, k]
      M[j, k] = x


def findpivot(M, index, k):
   """ find pivot element """
   # k is a_kk which would be a pivot element, but is zero
   for i in range(k, M.shape[0]): # all remaining rows
      for j in range(k, M.shape[1]): # all remaining cols
         if(M[i, j] != 0):
            return (1, i, j)
   # if we reach this point, all elements were zero and there
   # is no pivot element
   return (0, k, k)


def isSolvable(M):
   """ test if linear system of equations is solvable """
   m = M.shape[0]
   n = M.shape[1] - 1 # result vector b (Ax=b) inside matrix
   print m, n
   if m > n:
      return False;   # irregular!
   # check if rank A = rank A|b
   i = m - 1;
   while i >= 0 and M[i, i] == 0:
      if M[i, n] != 0:
         if verbose:
            print "rank different in " + str(i+1) + ": " + \
               str(M[i,i]) + "/" + str(M[i,n])
         return False;   # rank different!
      i = i - 1
   return True;
# end helpers mentioned in 1142 script


### functions imported from FeU Hagen Script 1142, Chapter 5.4
# lightly adopted for this framework

def gausselim(A, i, j):
   """ gausselimination on the given matrix for a given rowset """
   # i: start line, j: column in question
   if verbose:
      print "doing gausselim at " + str(i) + "," + str(j)
   for k in range(i + 1, A.shape[0]):
      A_kj=A[k,j]
      A_ij=A[i,j]
      Q_kj=-A_kj/A_ij
      A[k,:]=A[k,:]+Q_kj*A[i,:]
   if verbose:
      matPrint(A, 'gausselim result')

# not yet done!
def compute_x(C,ind,d):
   # C=[A|b], wobei A eine obere Dreiecksmatrix mit m <= n ist.
   # Wir loesen das System Ax=b unter Beruecksichtigung der
   # Permutation index fuer x.

   ## added rgerhards ##
   solvable = isSolvable(C)
   n = C.shape[0]
   m = C.shape[1]
   x = numpy.zeros((n, 1), numpy.float_, 'C')
   ## end rgerhards ##

   if solvable:
      x[ind[d]]=C[d,n]/C[d,d]
      for i in range(1,d+1):
         x[ind[d-i]]=C[d-i,n]
         for j in range(0,i):
            x[ind[d-i]]=x[ind[d-i]]-C[d-i,d-j]*x[ind[d-j]]
            x[ind[d-i]]=x[ind[d-i]]/C[d-i,d-i]
   return solvable, x

def gaussalg(matA, matB):
   A = matA['mat']
   m = matA['rows']
   n = matA['cols']
   d=min(m, n)
   index=range(0,n) # Permutation der Variablen
   C=numpy.concatenate((A,matB['mat']),1) # Verschmelze A,b nebeneinander
   for k in range(d):
      if verbose:
         print "pivot element at " + str(k) + " is: " + str(C[k,k])
      if C[k,k]==0:
         pivotfound,i,j = findpivot(C,index,k)
         if pivotfound == 0:
            break
         else:
            if i != k: # vertausche Zeilen i and k
               swaprows(C,i,k)
            if j != k: # vertausche Spalten
               swapcolumns(C,index,j,k) # und Variablen j and k
      gausselim(C,k,k)
   matA['mat'] = C
   matA['index'] = index
   solvable, x = compute_x(C,index,k)
   return solvable, x



### end functions imported from FeU Hagen Script 1142, Chapter 5.4

def solve(matA, matB):
   """ (try to) solve a system of linear equations """
   solvable, x = gaussalg(matA, matB)
   print "system is solvable: " + str(solvable)
   if solvable:
      matPrint(x, 'result matrix')
   

# ---------------------------------------------------------------------- #
# main program
print "Gauss 0.1, written by Rainer Gerhards (20090526)"
verbose = False
allMat = {}
words = getNewWords()
while words:
   cmd = words[0].lower()
   if cmd == "verbose":
         if words[1] == "on":
            verbose = True
         else:
            verbose = False
   elif cmd == "matrix":
      allMat[words[1]] = matCreate(words[1])
   elif cmd == "printmatrix":
      matPrint(allMat[words[1]]['mat'], words[1])
   elif cmd == "printallmatrix":
      matPrintAll()
   elif cmd == "duplicatematrix":
      matDup(words[1], words[2])
   elif cmd == "issolvable":
      print "matrix " + words[1] + " is solvable: " + \
            str(isSolvable(allMat[words[1]]['mat']))
   elif cmd == "gausselim":
      gausselim(allMat[words[1]]['mat'], int(words[2])-1, int(words[3])-1)
   elif cmd == "gaussalg":
      gaussalg(allMat[words[1]], allMat[words[2]])
      print "associated index: " + str(allMat[words[1]]['index'])
   elif cmd == "solve":
      solve(allMat[words[1]], allMat[words[2]])
   elif cmd == "swaprows":
      swaprows(allMat[words[1]]['mat'], int(words[2])-1, int(words[3])-1)
   elif cmd == "swapcols":
      index = range(0, allMat[words[1]]['cols'])
      swapcolumns(allMat[words[1]]['mat'], index,
             int(words[2])-1, int(words[3])-1)
      print "associated index: " + str(index)
   else:
      print "unexpected token in input: '" + cmd + "'"
   words = getNewWords()

# Steuerdatei für Beispiel 5.4.1
matrix A
dim 4 4
data
 1  1  1  1
-1  0  0  0
 0 -1  0  0
 0  0 -1  0
enddata
endmatrix

# result vector:
matrix b
dim 4 1
data
 10
-4
-3
-2
enddata
endmatrix

#now come commands
verbose off
printMatrix A
printMatrix b
solve A b
printMatrix A

Gauss 0.1, written by Rainer Gerhards (20090526)
Matrix C(4, 4):
[[ 1.  1.  1.  1.]
 [-1.  0.  0.  0.]
 [ 0. -1.  0.  0.]
 [ 0.  0. -1.  0.]]

Matrix b(4, 1):
[[ 10.]
 [ -4.]
 [ -3.]
 [ -2.]]

4 4
system is solvable: True
Matrix result matrix(4, 1):
[[ 4.]
 [ 3.]
 [ 2.]
 [ 1.]]

Matrix C(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.   0.   1.   1.   3.]
 [  0.   0.   0.   1.   1.]]

Gauss 0.1, written by Rainer Gerhards (20090526)
Matrix A(4, 4):
[[ 1.  1.  1.  1.]
 [-1.  0.  0.  0.]
 [ 0. -1.  0.  0.]
 [ 0.  0. -1.  0.]]

Matrix b(4, 1):
[[ 10.]
 [ -4.]
 [ -3.]
 [ -2.]]

pivot element at 0 is: 1.0
doing gausselim at 0,0
Matrix gausselim result(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.  -1.   0.   0.  -3.]
 [  0.   0.  -1.   0.  -2.]]

pivot element at 1 is: 1.0
doing gausselim at 1,1
Matrix gausselim result(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.   0.   1.   1.   3.]
 [  0.   0.  -1.   0.  -2.]]

pivot element at 2 is: 1.0
doing gausselim at 2,2
Matrix gausselim result(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.   0.   1.   1.   3.]
 [  0.   0.   0.   1.   1.]]

pivot element at 3 is: 1.0
doing gausselim at 3,3
Matrix gausselim result(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.   0.   1.   1.   3.]
 [  0.   0.   0.   1.   1.]]

4 4
system is solvable: True
Matrix result matrix(4, 1):
[[ 4.]
 [ 3.]
 [ 2.]
 [ 1.]]

Matrix A(4, 5):
[[  1.   1.   1.   1.  10.]
 [  0.   1.   1.   1.   6.]
 [  0.   0.   1.   1.   3.]
 [  0.   0.   0.   1.   1.]]