C Gaussian Elimination Implementation

sombriks, I also wrote a parallel version of gaussian elimination using MPI. My idea now is to benchmark them.

Cool! I’ve pasted here:

		import sys

		#############

		# Buggy Gauss Elimination

		# (Zero division is not checked)

		#

		def det(rows):

			 v = None

			 if len(rows) == 2:

				  r1 = rows[0]

				  r2 = rows[1]

				  v = r1[0] * r2[1] - r1[1] * r2[0]

			 else:

				  firstRow = rows[0]

				  aboveRows = rows[1:]        

				  subDets = []

				  # At time I din't know the existence of enumerate

				  for c in range(0, len(firstRow)):

				      subMatrix = []

				      for ar in aboveRows:

				          subRow = []

				          for c2 in range(0, len(ar)):

				              if c != c2:

				                  subRow.append(ar[c2])

				          subMatrix.append(subRow)

				      subDets.append(det(subMatrix) * firstRow[c])

				  evens = [subDets[e] for e in range(0, len(subDets), 2)]

				  odds = [subDets[e] for e in range(1, len(subDets), 2)]

				  v = reduce(lambda x, y: x+y, evens) - reduce(lambda x, y: x+y, odds)        

			 return v

		#non-recursive                

		def solveSystem(rows):

			 if det(rows) == 0:

				  return None

			 zerorColumnNth = 0

			 solvedSystem = rows

			 for workRowNth in range(0, len(rows)-1):

				  for prodRowNth in range(workRowNth+1, len(rows)):

				      workRow = solvedSystem[workRowNth]

				      prodRow = solvedSystem[prodRowNth]

				      mul = -prodRow[zerorColumnNth] / workRow[zerorColumnNth] 

				      newProdRow = map(lambda a, b: a + b, map(lambda x: x*mul, workRow), prodRow)

				      solvedSystem[prodRowNth] = newProdRow

				  zerorColumnNth += 1

			 return solvedSystem

		#recursive   

		def solveSystem2(rows):    

			 def solveLoop(rows, workNth, prodNth, zeroColumn, nrows):            

				  if prodNth < nrows:            

				      workRow = rows[workNth]

				      prodRow = rows[prodNth]

				      mul = -prodRow[zeroColumn] / workRow[zeroColumn] 

				      rows[prodNth] = map(lambda a, b: a + b, map(lambda x: x*mul, workRow), prodRow)

				      return solveLoop(rows, workNth, prodNth+1, zeroColumn, nrows)

				  else:

				      if workNth = 0:

				      row = rs[nth] 

				      vt = [vals[v]*row[v] for v in range(nth+1, len(row)-1)]                

				      if len(vt) > 0:

				          coff = -reduce(lambda a, b: a+b,  vt)

				      else:

				          coff = 0

				      vals[nth] = (coff + row[len(row)-1])/row[nth]

				      return retroLoop(rs, nth-1, vals)

				  return vals

			 nrows = len(rows)

			 return retroLoop(rows, nrows-1, [0 for i in range(0, nrows)])

		if __name__ == '__main__':

			 try:        

				  f = open(sys.argv[1], 'r')

			 except IOError, msg:

				  print(msg)

				  sys.exit(1)

			 rows = []

			 for line in f:

				  if line != '\n':

				      atoms = line.split(' ')

				      r = []

				      for a in atoms:

				          r.append(float(a))

				      rows.append(r)

			 f.close();

			 maxColumn = 0

			 for r in rows:

				  cs = len(r)

				  if cs > maxColumn:

				      maxColumn = cs

			 for r in rows:

				  cs = len(r)

				  if cs != maxColumn:

				      print("Matrix nao quad")

				      sys.exit(1)

			 ss = solveSystem2(rows)

			 r = retroSub(ss)

			 print(r)

		##############

		# Sample input file:

		#

		# 10 2 3 4 5

		# 6 17 8 9 10

		# 11 12 23 14 15

		# 16 17 18 29 20

		# Should output:

		# [0.28248587570621464, 0.24858757062146891, 0.24293785310734467, 0.23728813559322035]