#
# matrix.py
#
# Module for handling and calculating with Matrices
# 
# Author : Christian Hopp
#

#
# - Matrizen werden representiert durch Listen von Spaltenvektoren
#
# - Koordinaten werden gegeben als Paar aus Spalte und Zeile
#

# Achtung : bei solve_by_iterate2 ist keine SICHERE Überprüfung der
# Konvergenz gegeben! 


import types
import math
import cmath

def isMatrix(obj):
	return hasattr(obj,'matrix') 

def isMatrix_and_except(obj, col=None, lin=None):
	if not hasattr(obj,'matrix'):
		raise Matrix.Error,  "'Not a Matrix' error!"
	if (( col != None ) and ( lin != None)):
		if (( col != obj.col ) and ( lin != obj.lin)) :
			raise Matrix.Error,  "The matrices are of the wrong size!"
	return 1

def check_and_except(func):
	value = func()
	if (value):
		return value
	else:
		exception = "Error while checking : " + repr(func)
		raise Matrix.Error, exception
	
		
class Matrix:
#{{{ Konstruktor
	Error="Matrix Error"

	def __init__(self , col=0 , lin=0):
		if (type(col)==types.ListType):
			lin_ok=1
			for item in col:
				if len(col[0])!=len(item):
					raise Matrix.Error, "The listentries must have"\
						  "equal sizes!"
			
			self.matrix=col
			self.lin=len(col[0])
			self.col=len(col)
		elif isMatrix(col):
			self.col = col.col
			self.lin = col.lin
			self.matrix = []
			item = []
			for counter in range(col.lin):
				item.append(0)
			for counter in range(col.col):
				self.matrix.append(item[:])
		else:
			self.col = col
			self.lin = lin
			self.matrix = []
			item = []
			for counter in range(lin):
				item.append(0)
			for counter in range(col):
				self.matrix.append(item[:])
		return None

#}}}
#{{{ Verschiedene Operatoren

#{{{ __getitem__         --- Elementen abfragen

	def __getitem__(self, key):
		return self.matrix[key[0]][key[1]]

#}}}
#{{{ __setitem__         --- Elementen setzen

	def __setitem__(self, key, value):
		self.matrix[key[0]][key[1]]=value
		return None

#}}}
#{{{ __str__             --- Matrix in String umwandeln

	def __str__(self):
		return str(self.matrix)

#}}}
#{{{ __len__             --- Größe der Matrix

	def __len__(self):
		return (self.col*self.lin)-1

#}}}
#{{{ __getslice__        --- Ein Slice zurückliefern

	def __getslice__(self, leftkey, rightkey):
		if (type(leftkey) == types.TupleType):
			_leftkey = leftkey
		else:
			_leftkey = [int(leftkey)/int(self.col),int(leftkey)%int(self.col)] 
		if (type(rightkey) == types.TupleType):
			_rightkey = rightkey
		else:
			_rightkey = [int(rightkey)/int(self.col),int(rightkey)%\
						 int(self.col)]
			
		print _rightkey,_leftkey
		return self.partialmatrix(_leftkey,_rightkey)

#}}}
#{{{ __nonzero__         --- Matrix nicht gleich der Nullmatrix
	
	def __nonzero__(self):
		return (self!=Matrix(self.col,self.len))
	
#}}}


#}}}
#{{{ Rechenoperatoren

#{{{ __cmp__             --- Matrix-Vergleich

	def __cmp__(self, item):
		isMatrix_and_except(self)
		isMatrix_and_except(item,self.col,self.lin)
		for i in range(self.col):
			for j in range(self.lin):
				if self[i,j] != item[i,j] :
					return 1
		return 0

#}}}
#{{{ __neg__             --- Matrix-Negation

	def __neg__(self):
		return self.func_on_entry(lambda x : -x)

#}}}
#{{{ __add__             --- Matrix-Addition
	
	def __add__(self, item):
		isMatrix_and_except(self)
		isMatrix_and_except(item,self.col,self.lin)
		return_value = Matrix ( self.col , self.lin );
		for i in range(self.col):
			for j in range(self.lin):
				return_value[i,j] = self[i,j] + item[i,j]
		return return_value
	
#}}}
#{{{ __radd__            --- Rechtsseitige Matrix-Addition

	__radd__ = __add__

#}}}
#{{{ __sub__             --- Matrix-Subtraktion

	__sub__ = lambda a,b : a+(-b)
 
#}}}
#{{{ __rsub__            --- Rechtsseitige Matrix-Subtraktion

	__rsub__ = lambda a,b : (-a)+b

#}}}
#{{{ __mul__             --- Matrix-Multipliktion

	def __mul__(self, item):
		if (type(item) == type (1)) or (type(item) == type(1.0)) \
		   or (type(item) == type(1.0j)):
			return_value = Matrix ( self.col, self.lin );
			for i in range(self.lin):
				for j in range(self.col):
					return_value[j,i] = item * self[j,i]
		else:
			isMatrix_and_except(item)
			if ( self.col != item.lin ) :
				raise Matrix.Error,  "The matrices are of the wrong size!"
			return_value = Matrix ( item.col, self.lin );
			for i in range(self.lin):
				for j in range(item.col):
					value = 0
					for k in range(self.col) :
						value = value + self[k,i] * item[j,k]
						return_value[j,i] = value
		return return_value

#}}}
#{{{ __rmul__            --- Rechtsseitige Matrix-Multipliktion (!Types)

	def __rmul__(self, item):
		if (type(item) == type (1)) or (type(item) == type(1.0)) \
		   or (type(item) == type(1.0j)):
			return self*item

#}}}
#{{{ __div__             --- Matrix-Division

	def __div__(self, item):
		if ( type(item) == type (1) ) or ( type(item) == type(1.0)) \
		   or (type(item) == type(1.0j)):
			return self*(1.0/item)
		else:
			return self*item.inverse()

#}}}
#{{{ __rdiv__            --- Rechtsseitige Matrix-Division  (!Types)

	def __rdiv__(self, item):
		if ( type(item) == type (1) ) or ( type(item) == type(1.0)) \
		   or (type(item) == type(1.0j)):
			return item * self.inverse()

#}}}
#{{{ __pow__             --- Potentieren von Matrizen

	def __pow__(self, item):
		if (type(item) == type (1)) or (type(item) == type(1.0))\
		   and (item == int(item)):
			if (item == 0):
				return self.E_matrix()
			if (item == 1):
				return self
			if (self.col != self.lin):
				raise Matrix.Error,  "Matrix does not have equal side"\
					  "lengthes!"
			if (item < 0 ):
				return self.inverse()**abs(item)
			else:
				return self*(self**(item-1))
		else:
			raise Matrix.Error,  "The Exponent must be an positive or"\
				  " negativ, natural Number !"

#}}}

#}}}
#{{{ einige Methoden

#{{{ E_matrix            --- Einheitsmatrix zurück geben

	def E_matrix(self, col=None,lin=None):
		if (lin == None) and (col != None):
			lin = col
		if (col == None) and (lin != None):
			col = lin
		if (col == None):
			col=self.col
		if (lin == None):
			lin=self.col
		return_value = Matrix(col,lin)
		for i in range(col):
			for j in range(lin):
				if i==j:
					return_value[i,j] = 1.0
				else:
					return_value[i,j] = 0.0
		return return_value

#}}}

#{{{ func_on_entry       --- Die Elemente der Matrix jeweils mit
#                            einer Funktion bearbeiten

	def func_on_entry(self,func):
		if (type(func) != types.FunctionType) and \
		   (type(func) != types.BuiltinFunctionType) and \
		   (type(func) != types.MethodType):
			raise Matrix.Error, "Parametrer must be a function"
		
		return_value = Matrix(self.col,self.lin)
		for i in range(self.col):
			for j in range(self.lin):
					return_value[i,j] = func(self[i,j])
		return return_value

#}}}
#{{{ entry_to_int        --- Elemente in int konvertieren

	def entry_to_int(self):
		return self.func_on_entry(int)

#}}}
#{{{ entry_to_float      --- Elemente in float konvertieren

	def entry_to_float(self):
		return self.func_on_entry(float)

#}}}
#{{{ entry_to_long       --- Elemente in long konvertieren

	def entry_to_long(self):
		return self.func_on_entry(long)

#}}}

#{{{ partialmatrix       --- Bilden Teilmatrix nach Koordinaten

	def partialmatrix(self,i,j):
		if i[0]>j[0]:
			# i[0],j[0] = j[0],i[0]
			return Matrix(0,0)
		if i[1]>j[1]:
			#i[1],j[1] = j[1],i[1]
			return Matrix(0,0)
		if (i[0]<0) or (i[1]<0) or (j[0]>self.col) or (j[1]>self.lin):
			raise Matrix.Error,  "Index is out of range!"
		return_value = Matrix(j[0]-i[0]+1, j[1]-i[1]+1)
		for k in range(j[0]-i[0]+1):
			for l in range(j[1]-i[1]+1):
				return_value[k,l] = self[i[0]+k,i[1]+l]
		return return_value

#}}}
#{{{ concatmatrix_by_col --- Zusammenfügen zweier Matrizen gegen Spalten

	def concatmatrix_by_col(a,b):
		if a.col == 0 and a.lin == 0 :
			if b.col == 0 and b.lin == 0 :
				return Matrix(0,0)
			else:
				return b
		else:
			if b.col == 0 and b.lin == 0 :
				return a
		if (a.lin != b.lin):
			raise Matrix.Error,  "Matrices do not have equal linelength!"
		return_value = Matrix(a.col+b.col, a.lin)
		for j in range(a.lin):
			for i in range(a.col):
				return_value[i,j] = a[i,j]
			for i in range(b.col):
				return_value[a.col+i,j] = b[i,j]
		return return_value

#}}}
#{{{ concatmatrix_by_lin --- Zusammenfügen zweier Matrizen gegen Zeilen

	def concatmatrix_by_lin(a,b):
		if a.col == 0 and a.lin == 0 :
			if b.col == 0 and b.lin == 0 :
				return Matrix(0,0)
			else:
				return b
		else:
			if b.col == 0 and b.lin == 0 :
				return a
		if (a.col != b.col):
			raise Matrix.Error,  "Matrices do not have equal columnlengthes!"
		return_value = Matrix(a.col, b.lin + a.lin)
		for i in range(a.col):
			for j in range(a.lin):
				return_value[i,j] = a[i,j]
			for j in range(b.lin):
				return_value[i,a.lin+j] = b[i,j]
		return return_value

#}}}

#{{{ determinant_recurse --- Determinante (nach der Zeile/Spalte mit
#                            möglichst vielen Nullen)

	def determinant_recurse(self, lin=-1, col=-1):
		check_and_except(self.quadratic)
		
		if (lin >= 0) and (col >= 0):
			raise Matrix.Error,  "Specify only a Line or a Column!"
		if (lin > self.lin-1) or (col > self.col-1):
			raise Matrix.Error,  "Index out of Range!"
		if (self.col == 0): # det einer Matrix(0,0) = 0 ?
			return 0
		if (self.col == 1):
			return self[0,0]
		
		#Beste Zeile oder Spalte finden
		max_lin=0
		max_col=0
		by_col=0

		if (lin >= 0):
			by_col = 0
		elif (col >= 0):
			by_col = 1
		else:  # Suche nach der Zeile/Spalte nach der entwickelt wird
			if (self.col > 4):
				for i in range(self.col):
					sum_lin=0
					sum_col=0
					for j in range(self.col):
						if isMatrix(self[i,j]):
							by_col = 0
							break
						if self[i,j]==0:
							sum_lin=sum_lin+1
						if self[j,i]==0:
							sum_col=sum_col+1
					if sum_lin > max_lin:
						max_lin=sum_lin
						lin=i
					if sum_col > max_col:
						max_col=sum_col
						col=i
					
				if max_col > max_lin:
					by_col=1
			else:
				lin=0
				col=0
		if (lin < 0):
			lin = 0
		if (col < 0):
			col = 0

		det_func = lambda i,j,m : (- 2 * (((i + j) % 2 ) - 0.5)) * m[i,j] * \
				   m.submatrix(i,j).determinant()
		
		if by_col==1: # Entwicklung nach Spalten			
			return_value = det_func(0,col,self)
			for i in range(1,self.col):
				if isMatrix(self[lin,i]) or self[i,col]!=0:
					return_value = return_value + det_func(i,col,self)
				else:
					return_value = return_value + self[i,col]
		else:         # Entwicklung nach Zeilen
			return_value = det_func(lin,0,self)
			for i in range(1,self.col):
				if isMatrix(self[lin,i]) or self[lin,i]!=0:
					return_value = return_value + det_func(lin,i,self)
				else:
					return_value = return_value + self[lin,i]
			
		return return_value

#}}}
#{{{ determinant_gauss   --- Determinante durch LR-Zerlegung
	def determinant_gauss(self, lin=-1, col=-1):
		check_and_except(self.quadratic)
		try:
			LR           = self.gauss_LR_apart()
		except ZeroDivisionError:
			return 0
			
		R            = LR[1]
		return_value = 1.0
		
		for i in range(self.col):
			return_value = return_value * R[i,i]
		return return_value
		

#}}}
#{{{ determinant         --- ---> determinant_recurse

	determinant = determinant_recurse

#}}}

#{{{ symetric            --- Überprüfen auf Symetrie

	def symetric (self):
		if (self == self.transpose()):
			return 1==1
		else:
			return 1==0
		
#}}}	
#{{{ quadratic           --- Überprüfen ob quadratisch

	def quadratic (self):
		if (self.col == self.lin):
			return 1==1
		else:
			return 1==0
		
#}}}	

#{{{ submatrix           --- Bilden der i,j 'ten Untermatrix

	def submatrix (self, i,j):
		check_and_except(self.quadratic)

		if (self.col == 1):
			return Matrix(0,0)
		
		a = self.partialmatrix([0,0],[i-1,j-1])
		b = self.partialmatrix([i+1,0],[self.col-1,j-1])
		c = self.partialmatrix([0,j+1],[i-1,self.lin-1])
		d = self.partialmatrix([i+1,j+1],[self.col-1,self.lin-1])
		
		return Matrix.concatmatrix_by_lin ( Matrix.concatmatrix_by_col(a,b), \
											Matrix.concatmatrix_by_col(c,d))

#}}}
#{{{ adjunct             --- Bilden der adjunktierten Matrix

	def adjunct (self):
		check_and_except(self.quadratic)

		if (self.col == 1): # Adjunkte einer Matrix(1,1) ?
			return Matrix().E_matrix(1,1)
		return_value = Matrix(self.col,self.lin)
		for i in range(self.col):
			for j in range(self.lin):
				return_value[i,j] = (- 2 * (((i+j) % 2 ) - 0.5 )) *\
									self.submatrix(i,j).determinant()
		return return_value.transpose()

#}}}
#{{{ scalar_prod         --- Bilden des Produkt mit einem Skalar (obsolete!)

	def scalar_prod (self, scalar):
		return_value = Matrix(self.col,self.lin)
		for i in range(self.col):
			for j in range(self.lin):
				return_value[i,j] = self[i,j]*scalar
		return return_value

#}}}

#{{{ inverse_adjunct     --- Bilden der Inversen über die adjunktierte Matrix

	def inverse_adjunct (self):
		check_and_except(self.quadratic)
		if (self.col == 1):
			a = Matrix(1,1)
			a[0,0] = -1 * self[0,0]
			return a
		det = check_and_except(self.determinant)
		return self.adjunct()*(1/det)

#}}}
#{{{ inverse_swap        --- Bilden der Inversen mit dem Austauschverfahren 

	def inverse_swap (self):
		# Prüfen ob möglich!
		check_and_except(self.quadratic)
		if (self.col == 1):
			a = Matrix(1,1)
			a[0,0] = -1 * self[0,0]
			return a
		# Anlegen der Daten
		a_old = self
		a_new = Matrix(self.col, self.lin)
		x_new = []
		x_old = []
		y_new = []
		y_old = []
		for i in range(self.col):
			x_new[i:] = [0]
			y_new[i:] = [0]
			x_old[i:] = [i]
			y_old[i:] = [i]
		# print x_old,y_old
		for i in range(self.col):
			# Finden des Pivotelement
			k  = 0
			l  = 0
			
			for k in range(self.col):
				ok = 0
				for l in range(self.col):
					if (x_old[k] >= 0) and (y_old[l] >= 0) and\
					   (a_old[k,l] != 0):
						x_old[k]= -5
						x_new[k]= l
						y_old[l]= -5
						y_new[l]= k
						ok = 1
						# print "|",k,l,"|", x_old,y_old
						break
					elif (k==5) and (l==5) and (a_old[k,l] == 0):
						raise Matrix.Error,  "Inverting Problem!"
				if ok==1:
					break
			# Austauschen des Pivotelement
			for m in range(self.col):
				for n in range(self.lin):
					try:
						if (m==k) and (n==l):
							a_new[m,n]=1.0/a_old[m,n]
						elif (n==l) and (m!=k):
							a_new[m,n]=1.0*a_old[m,n]/a_old[k,l]
						elif (m==k) and (n!=l):
							a_new[m,n]=-1.0*a_old[m,n]/a_old[k,l]
						else:
							a_new[m,n]=a_old[m,n]-1.0*a_old[m,l]*\
										a_old[k,n]/a_old[k,l]
					except ZeroDivisionError:
						raise Matrix.Error,  "Matrix has no inverse!"
			a_old=a_new
			# print "ende:",l,k, a_old, x_old, y_old, x_new, y_new
			a_new=Matrix(self.col,self.lin)
		# Zurück sortieren
		for i in range(self.col-1):
			min_x = i
			min_y = i
			for j in range(i+1,self.col):
				if (x_new[j] < x_new[min_x]) :
					min_x = j
				if (y_new[j] < y_new[min_y]) :
					min_y = j
			temp = x_new[min_x]
			x_new[min_x] = x_new[i]
			x_new[i] = temp
			temp = y_new[min_y]
			y_new[min_y] = y_new[i]
			y_new[i] = temp
			for j in range(self.col):
				temp=a_old[min_x,j]
				a_old[min_x,j]=a_old[i,j]
				a_old[i,j]=temp
			for j in range(self.col):
				temp=a_old[j,min_y]
				a_old[j,min_y]=a_old[j,i]
				a_old[j,i]=temp
		return a_old  

#}}}
#{{{ inverse_gauss       --- Bilden der Inversen mit dem Gauss-Algorhythmus

	def inverse_gauss (self):
		check_and_except(self.quadratic)
		if (self.col == 1):
			a = Matrix(1,1)
			a[0,0] = -1 * self[0,0]
			return a
		return_value = Matrix(0,self.lin)
		e = Matrix(1,self.lin)
		
		LR = self.gauss_LR_apart()
		L,R = LR[0],LR[1]

		for i in range(self.lin):
			e[0,i] = 1
			try:
				v = R.backwardinsert(e.gauss_get_r(L))
			except ZeroDivisionError:
				raise Matrix.Error,  "Matrix has no inverse!"
			e[0,i] = 0
			return_value = Matrix.concatmatrix_by_col(return_value,v)
		return return_value

#}}}	
#{{{ inverse             --- ---> inverse_swap 

	inverse = inverse_swap

#}}}
#{{{ transpose           --- Transponieren der Matrix

	def transpose (self):
		return_value = Matrix(self.lin,self.col)
		for i in range(self.col):
			for j in range(self.lin):
				return_value[j,i] = self[i,j]
		return return_value

#}}}

#{{{ cholesky            --- Bilden der Choleskyzerlegung
	
	def cholesky (self):
		check_and_except(self.quadratic)
		check_and_except(self.symetric)
		return_matrix = Matrix(self.col, self.lin)
		
		for k in range(self.col):
			for j in range(k+1):
				if (j == k):         # Diagonalelement
					value = self[k,j] * 1.0
					
					for i in range(k):
						value = value - return_matrix[k,i]\
								* return_matrix[k,i]
						
					if (value > 0):
						value = pow(value,0.5)
					else:
						raise Matrix.Error,  "Matrix is not positiv definit"
				else:                # Nicht Diagonalelement
					value = self[k,j] *1.0
					
					for i in range(j):
						value = value - return_matrix[j,i]\
								* return_matrix[k,i]
						
					value = value / return_matrix[j,j]
					
				return_matrix[k,j] = value
				
		return return_matrix
		
#}}}	
#{{{ backwardinstert     --- Rückwärtseinsetzen von Vektoren in Matrizen

	def backwardinsert (self,x):
		if x.col != 1:
			raise Matrix.Error, "It's not a Vector!"
		
		if self.col != x.lin:
			raise Matrix.Error, "Matrix and Vector are not of the same size!"
		
		if (self.col != self.lin):
			raise Matrix.Error,  "Matrix does not have equal side lengthes!"
		
		return_vector = Matrix(1,self.lin)
		return_vector[0,self.col-1] = x[0,self.col-1] / self[self.lin-1,\
															 self.col-1]
		
		for i in range(self.col-2,-1,-1):
			value = x[0,i] * 1.0
			
			for j in range(i+1,self.col):
				value = value - self[j,i] * return_vector[0,j]
				
			return_vector[0,i] = value / self[i,i]
			
		return return_vector
	
#}}} 
#{{{ forewardinstert     --- Vorwärtseinsetzen von Vektoren in Matrizen

	def forewardinsert (self,x):
		if x.col != 1:
			raise Matrix.Error, "It's not a Vector!"
		
		if self.col != x.lin:
			raise Matrix.Error, "Matrix and Vector are not of the same size!"
		
		if (self.col != self.lin):
			raise Matrix.Error,  "Matrix does not have equal side lengthes!"
		
		return_vector = Matrix(1,self.lin)
		return_vector[0,self.col-1] = x[0,self.col-1] / self[self.lin-1,\
															 self.col-1]
		
		for i in range(self.col):
			value = x[0,i] * 1.0
			
			for j in range(i):
				value = value - self[i,j] * return_vector[0,j]
				
			return_vector[0,i] = value / self[i,i]
			
		return return_vector
	
#}}} 
#{{{ solve_by_gauss      --- Lösen von einem LGS mit dem Gaußalgorhithmus 

	def gauss_LR_apart (self):
		L = Matrix(self.col, self.lin)
		R = Matrix(self.col, self.lin)
		
		for j in range(self.lin):
			for k in range(j,self.lin):
				value = self[k,j] * 1.0
				for p in range(j):
					value = value - L[p,j] * R[k,p]
				R[k,j] = value
			for i in range(j+1,self.lin):
				value = self[j,i] * 1.0
				for p in range(j):
					value = value - L[p,i] * R[j,p]
				L[j,i] = value / R[j,j]

		return_value = [L,R]
		return return_value

	def gauss_get_r (self,L):
		if self.col != 1:
			raise Matrix.Error, "Not a Vector!"

		r = Matrix(1,self.lin)
		
		for j in range(L.lin):
			value = self[0,j] * 1.0
			for p in range(j):
				value = value - L[p,j] * r[0,p]
			r[0,j] = value
		return r	

	def solve_by_gauss (self,x):
		if x.col != 1:
			raise Matrix.Error, "Not a Vector!"
		
		if self.col != x.lin:
			raise Matrix.Error, "Matrix and Vector are not of the same size!"
		
		if (self.col != self.lin):
			raise Matrix.Error,  "Matrix does not have equal side lengthes!"

		LR = self.gauss_LR_apart()
		L = LR[0]
		R = LR[1]
		r = x.gauss_get_r(L)
			
		return R.backwardinsert(r)

#}}}
#{{{ solve_by_cholesky   --- Lösen von einem LGS durch Choleskyzerlegung

	def solve_by_cholesky (self,x):
		chol = self.cholesky()
		return chol.backwardinsert(chol.forewardinsert(x))

#}}}
#{{{ solve_by_iterate1   --- Lösen von einem LGS durch Iterieren mit
#                            Gesamtschritt
	def solve_by_iterate1 (self,b,x_start,n_max):
		if b.col != 1:
			raise Matrix.Error, "It's not a Vector!"
		
		if self.col != b.lin:
			raise Matrix.Error, "Matrix and Vector are not of the same size!"
		
		if (self.col != self.lin):
			raise Matrix.Error,  "Matrix does not have equal side lengthes!"
		
		# Überprüfen der Konvergenz
		
		max_1 = 0
		max_2 = 0
		
		for i in range(self.lin):
			value_1 =  0
			value_2 =  0
			for k in range(self.col):
				if not i==k:
					value_1 = value_1 + abs(self[k,i] * 1.0 /self[i,i])
					value_2 = value_2 + abs(self[i,k] * 1.0 /self[k,k])
			if ( value_1 > max_1 ):
				max_1 = value_1
			if ( value_2 > max_2 ):
				max_2 = value_2
		
		if (max_1 >= 1) and (max_2 >= 1):
			raise Matrix.Error, "Not konvergent!"
		
		x_new = Matrix(1,self.lin)
		x_old = Matrix(1,self.lin)
		
		for k in range(self.lin):
			x_old[0,k] = x_start[0,k]
		
		# Iteration
		
		for iter in range(n_max):
			for i in range(self.lin):
				value = b[0,i] * 1.0 / self[i,i]
				for k in range(self.lin):
					if not i==k:
						value = value - 1.0 * x_old[0,k] *  self[k,i] /\
								self[i,i]
				x_new[0,i] = value
			for k in range(self.lin):
				x_old[0,k] = x_new[0,k]
				
		return x_old

#}}}
#{{{ solve_by_iterate2   --- Lösen von einem LGS durch Iterieren mit\
#                            Einzelschritt
	def solve_by_iterate2 (self,b,x_start,n_max):
		if b.col != 1:
			raise Matrix.Error, "It's not a Vector!"
		
		if self.col != b.lin:
			raise Matrix.Error, "Matrix and Vector are not of the same size!"
		
		if (self.col != self.lin):
			raise Matrix.Error,  "Matrix does not have equal side lengthes!"
		
		# Überprüfen der Konvergenz
		
		max_1 = 0
		max_2 = 0
		
		for i in range(self.lin):
			value_1 =  0
			value_2 =  0
			for k in range(self.col):
				if not i==k:
					value_1 = value_1 + abs(self[k,i] * 1.0 /self[i,i])
					value_2 = value_2 + abs(self[i,k] * 1.0 /self[k,k])
			if ( value_1 > max_1 ):
				max_1 = value_1
			if ( value_2 > max_2 ):
				max_2 = value_2
		
		if (max_1 >= 1) and (max_2 >= 1):
			raise Matrix.Error, "Not convergend!"
		
		x_new = Matrix(1,self.lin)
		x_old = Matrix(1,self.lin)
		
		for k in range(self.lin):
			x_old[0,k] = x_start[0,k]
		
		# Iteration
		
		for iter in range(n_max):
			for i in range(self.lin):
				value = b[0,i] * 1.0 / self[i,i]
				for k in range(i):
					value = value - 1.0 * x_new[0,k] *  self[k,i] / self[i,i]
				for k in range(i+1,self.lin):
					value = value - 1.0 * x_old[0,k] *  self[k,i] / self[i,i]
				x_new[0,i] = value
			for k in range(self.lin):
				x_old[0,k] = x_new[0,k]
				
		return x_old

#}}}

#{{{ lin_sum_norm        --- Zeilensummennorm
	
	def lin_sum_norm (self):
		max = 0
		for i in range(self.col):
			value = 0
			for j in range(self.lin):
				value = value + abs(self[i,j])
			if ( value > max ):
				max = value
		return max
	
#}}} 
#{{{ col_sum_nom         --- Spaltensummennorm
	
	def col_sum_norm (self):
		max = 0
		for i in range(self.lin):
			value = 0
			for j in range(self.col):
				value = value + abs(self[j,i])
			if ( value > max ):
				max = value
		return max
	
#}}} 

#{{{ exp                 --- Exponentialfunktion über eine Matrix nach\
#                            Exponentialreihenentw.

	def exp_ok(self, n=10):
		check_and_except(self.quadratic)
		return_value = self.E_matrix()
		power_matrix = self.E_matrix()
		float_self   = self.entry_to_float()

		for i in range(1,n+1):
			power_matrix = power_matrix * float_self / i
			return_value = return_value + power_matrix 
		return return_value

	def exp_ok2(self, n=10,limit=1e-5):
		check_and_except(self.quadratic)
		return_value = self.E_matrix()
		power_matrix = self.E_matrix()
		float_self   = self.entry_to_float()
		i            = 0
		
		while (i<n) or (power_matrix.col_sum_norm()>limit) :
			i=i+1
			power_matrix = power_matrix * float_self / i
			return_value = return_value + power_matrix
		
		return return_value

	def exp(self, n=3, max_n = 100,limit=1e-5, max_exponent = 6):
		check_and_except(self.quadratic)
		return_value = self.E_matrix()
		power_matrix = self.E_matrix()
		exponent     = int(math.sqrt(self.col_sum_norm()))
		
		if ( exponent < 2 ):
			exponent = 1
		elif ( exponent > max_exponent ):
			exponent = max_exponent

		float_self   = self.entry_to_float()/exponent
		
		i=0
		while ((i<n) or (power_matrix.col_sum_norm()>limit)) and (i<max_n):
			i=i+1
			power_matrix = power_matrix * float_self / i
			return_value = return_value + power_matrix
		return return_value**exponent

#}}}

#{{{ perpendicular       --- Bilden von eines senkrechten Vektor auf einer
#                            von n-1 Vektoren auf gespannten Hyperebene
#                            (für 2  [1x3]-Matrizen => Kreuzprodukt)

	def perpendicular (self,*vectors):
		for item in vectors:
			if vectors[0].lin != item.lin:
				raise Matrix.Error, "The Vectors must have equal sizes!"
			if 1 != item.col:
				raise Matrix.Error, "Some Vectors are not Columnvectors!"
			
		if len(vectors) != (vectors[0].lin - 1):
			raise Matrix.Error, "n-1 Vectors are needed!"
		
		det_matrix = Matrix(vectors[0].lin,vectors[0].lin)
		
		for i in range(vectors[0].lin):
			det_matrix[i,0] = Matrix().E_matrix(vectors[0].lin).\
							  partialmatrix([i,0],[i,vectors[0].lin-1])
			for j in range(vectors[0].lin - 1):
				det_matrix[i,j+1] = vectors[j][0,i]
				
		return det_matrix.determinant(lin=0)

#}}}
#}}}