<reponame>doctorsrn/hw_competition
"""Paquete que muestra la resolucion y factorizacion de un sistema de ecuaciones lineales
	 aplicando diferentes metodos.

   Nota: 
	- La mayoria de las funciones y las variables estan en ingles, y los comentarios en espanol.	
	- Aun falta implementar varios metodos, y mejorar los existentes.

        - El paquete no soporta matrices singulares.

   Definiciones:
   A, M : Matriz original.
   a : Matriz aumentada (A|b)
   b : Vector independiente. 
   L : Matriz triangular inferior unitaria.
   M : Matriz triangular inferior unitaria.
   U : Matriz triangular superior unitaria.
   D : Matriz diagonal.
   P : Matriz de permutacion
   Pr : Matriz de permutacion de filas
   Pc : Matriz de permutacion de columnas
   G : Matriz superior hallada con el metodo de Cholesky.
   Gt : Transpuesta de G.
   A = LU : significa factorizacion sin pivoteo
   PA = LU : significa factorizacion con pivoteo parcial
   PAQ = LU : significa factorizacion con pivoteo total

"""
from math import fabs
from math import sqrt
from sys import exit


def inputMatrix():
	"""Ingreso de datos de la matriz"""

	n = int(input("Ingrese el orden de la matriz: "))
	print "Ingrese los elementos de la matriz A fila por fila con un espacio luego enter"
        print "Por ejemplo, para una matriz de orden 3:","\n", "3 2 1"
        print "1 7 5","\n\n"
	A = [[0.0]*n for i in range(n)]
	aux = [[0.0]*n for i in range(n)]
	for i in range(n):
		temp = raw_input()
		A[i] = temp.split()
		for j in range(n):
			A[i][j] = float(A[i][j])
	if A == aux:
		exit('Matriz nula, vuelva a escribir la matriz')
	return A		

def inputVector(M):
	"""Ingreso de datos del vector independiente"""

	n = len(M)
	print 'Ingrese los elementos del vector b'
	b = [0]*n
	#for i in range(n):
	#	b[i] = float(input())

	temp = raw_input();
	b = temp.split()
	for i in range(n):
		b[i] = float(b[i])

	return b	
	
def printMatrixD(M):
	"""Imprime la matriz, version debugger."""

	for i in range(len(M)):
		# print '|',
		for j in range(len(M[i])):
			if(i == 0 and j == 0):
				print " ",
				for x in range(len(M[i])):
					print '{0:8}'.format(x),
				print
				print
			if(j == 0):
				print i,

			if(j == len(M)):
				print '|',
				print '{0:8.4f}'.format(M[i][j]),
			else:
				print '{0:8.4f}'.format(M[i][j]),	
		print '|'
	print

def printMatrix(M):
	"""Imprime la matriz de una forma legible."""

	for i in range(len(M)):
		print '|',
		for j in range(len(M[i])):
			if(j == len(M)):
				print '|',
				print '{0:8.4f}'.format(M[i][j]),
			else:
				print '{0:8.4f}'.format(M[i][j]),	
		print '|'
	print  

def norm1(A):
	"""Calcula la norma 1"""

	summation = 0
	for i in range(len(A)):
		summation += abs(A[i][0])

	for j in range(1, len(A)):
		temp = 0
		for i in range(len(A[i])):
			temp += abs(A[i][j])
		summation = temp if (temp > summation) else summation

	return summation 

def infinityNorm(A):
	"""Calcula la norma infinita"""

	summation = 0
	for j in range(len(A)):
		summation += abs(A[0][j])

	for i in range(1, len(A)):
		temp = 0
		for j in range(len(A[i])):
			temp += abs(A[i][j])
		summation = temp if (temp > summation) else summation
			
	return summation  

def symmetricMatrix(A):
	"""Calcula la matriz simetrica de A"""

	for i in range(len(A)):
		for j in range(i+1,len(A)):
			if A[i][j] != A[j][i]:
				return False
	return True 

def trans(M):
	"""Calcula la matriz transpuesta de M"""

	n = len(M)
	return [[ M[i][j] for i in range(n)] for j in range(n)]

def reverseSub(a):
	"""Se usa para resolver la matriz superior de gauss/gauss-Jordan"""

	n = len(a)
	x = [0]*n
	for j in range(n-1, -1, -1):
            if a[j][j] != 0:#Una manera de evitar una excepcion en matriz singular
                    x[j] = (a[j][n] - sum(a[j][k]*x[k] for k in range(j+1, n)))/float(a[j][j])
            else:
                print 'Matriz singular, la ecuacion no tiene solucion o tiene infinitas soluciones .'
                break
	return x	

def solMatrixSup(M, b):
	"""Resuelve la matriz superior con el metodo inverso"""	

	x = []
	for i in range(len(M)-1,-1,-1):
		x.append((1.0/(M[i][i]))*(b[i]-sum(M[i][len(M)-j-1]*x[j] for j in range(len(x)))))
	x.reverse()
	
	return x

def solMatrixInf(M, b):
	"""Resuelve la matriz inferior con el metodo inverso"""

	x = []
	for i in range(len(M)):
		x.append((1.0/(M[i][i]))*(b[i]-sum(M[i][j]*x[j] for j in range(len(x)))))
	return x
	
def augmentedMatrix(M, b):
	"""Retorna la matriz aumentada"""
	
	a = M
	for i in range(len(M)):
		a[i].append(b[i])

	return a	

def maxColum(M, c):
	"""Devuelve la fila del mayor valor debajo de la diagonal en la columna c"""

	r = c #fila
	maximum = M[c][c]
	for i in range(c+1,len(M)):
		if(fabs(maximum) < fabs(M[i][c])):
			maximum = M[i][c]
			r = i
	return r

def maxSubMatrix(M, c):
	"""Retorna el mayor elemento de la submatriz A[c]"""

	row = c
	colum = c
	n = len(M)
	maximum = M[c][c]
		
	for j in range(c, n):
		for k in range(c, n):
			maxTemp = M[k][j]
			if(fabs(maximum) < fabs(maxTemp)):
				maximum = maxTemp
				colum = j
				row = k

	return row, colum			

def exchangeRows(M, r1, r2):
	"""Intercambia las filas r1 y r2 de M"""

	M[r1], M[r2] = M[r2], M[r1]
	return M

def exchangeCols(M, c1, c2):
	"""Intercambia las columnas c1 y c2 de M"""

	for k in range(len(M)):
		M[k][c1] , M[k][c2] = M[k][c2], M[k][c1]
	return M	

def pivot(a, P, Q, colum, piv=0):
	"""Se encarga del pivoteo en cualquier metodo."""

	if piv > 2 or piv < 0:
		exit('Valores invalidos para el parametro pivoteo, valores validos: 0, 1, 2.')
	n = len(a)
		
	temp = a[colum][colum]
	if(temp == 0.0 and piv == 0):
		row_maxColumn = maxColum(a, colum)
		a = pivotP(a, row_maxColumn, colum)
		P = exchangeRows(P, row_maxColumn, colum)
		print 'P(%d,%d)' % (row_maxColumn, colum)
		printMatrix(Pr(n, row_maxColumn, colum))
		
	elif piv == 1:
		row_maxColumn = maxColum(a, colum)
		if row_maxColumn != colum:
			a = pivotP(a, row_maxColumn, colum)
			P = exchangeRows(P, row_maxColumn, colum)
			print 'P(%d,%d)' % (row_maxColumn, colum)
			printMatrix(Pr(n, row_maxColumn, colum))
			
	elif piv == 2:
		row, c = maxSubMatrix(a, colum)	
		if (row != colum) or (c != colum) :
			a = pivotT(a, colum)
			P = exchangeCols(P, row, colum)
			Q = exchangeCols(Q, colum, c)
			print 'P(%d,%d):' % (colum, row)
			printMatrix(Pr(n, row, colum))
			print 'Q(%d,%d):' % (colum, c)
			printMatrix(Pr(n, c, colum))
	
	#retorna la matriz, la matriz de permutacion de filas y la de columnas.	
	return 	a, P, Q		

def pivotP(M, r1, r2):
	"""Permuta la fila r1 con la fila r2 de la matriz M"""

	return exchangeRows(M, r1, r2)

def pivotT(M, i):
	"""Busca el mayor elemento de la submatriz A[i] y permuta filas y columnas"""

	r, c = maxSubMatrix(M, i)
	M = pivotP(M, i, r)		

	return exchangeCols(M, c, i)

def L(n, ri, c):
	"""Calcula las matrices Li.

		ri: a[i][c] / a[c][c]
		c : columna c-esima 
	"""
	#Matriz identidad
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	for k in range(c+1, n):
		L[k][c] = -ri[k]
			
	return L

def matrixMulti(A, B):
	"""Multiplica dos matrices, C = A*B """

	rowsA, colsA = len(A), len(A[0])
	rowsB, colsB = len(B), len(B[0])

	if colsA != rowsB:
		exit('Dimensiones incorrectas')

	C = [[0 for row in range(colsB)] for col in range(rowsA)]

	for i in range(rowsA):
		for j in range(colsB):
			for k in range(colsA):
				C[i][j] += A[i][k]*B[k][j]
	return C

def T(n, ri, c):
	"""Calcula las matrices Ti (Para el metodo Gauss-Jordan)

		ri: a[i][c] / a[c][c]
		c : columna c-esima 
	"""
	T = [[float(i == j) for j in range(n)] for i in range(n)]
	for k in range(n):
		T[k][c] = -ri[k]

	return T

def Pr(n, r1, r2):
	"""Calcula la matriz de permutacion de fila"""

	#Matriz identidad
	I = [[float(i == j) for j in range(n)] for i in range(n)]
	return exchangeRows(I, r1, r2)

def Pc(n, c1, c2):
	"""Calcula la matriz de permutacion de columna"""

	#Matriz identidad
	I = [[float(i == j) for j in range(n)] for i in range(n)]
	return exchangeCols(I, c1, c2)	

def multi(A,x):
	"""Multiplica una matriz con un vector"""

	c=[]
	for i in range(len(A)):
		c.append(sum(A[i][j]*x[j] for j in range(len(A))))
	return c

def invLU(A):
	P, L, U = croutLU1(A)
	for x in range(1,10):
		pass
		
def gauss(A, b, piv=0):
	"""Metodo de Gauss

		A: Matriz
		b: vector independiente
		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial
			 piv=2 con pivoteo total

	"""
	a = augmentedMatrix(A, b)
	n = len(A)
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]
	for j in range(len(a[0])):
		r = [0]*n
		# pivoteo
		if j < n:
			a, P, Q = pivot(a, P, Q, j, piv)

		for i in range(j+1, n):		
			
			#Eliminacion de gauss
			r[i] = a[i][j] / float(a[j][j])

			for k in range(j, len(a[0])):
				a[i][k] = a[i][k] - r[i]*a[j][k]
					
		if j < n:
			temp = [0]*n
			if (temp != r):
				print ''#'L%d' % (j+1)
#				printMatrix(L(n, r, j))			

#	print 'A|b:'			
#	printMatrix(a)

	if piv == 2:
		return multi(Q, reverseSub(a))	

	return reverseSub(a)			

def gaussJordan(A, b, piv=0):
	"""Metodo de Gauss-Jordan
		A: Matriz
		b: vector independiente
		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial
			 piv=2 con pivoteo total  	
	"""
	a = augmentedMatrix(A, b)
	n_rows = len(a)
	n_cols = len(a[0])
	P = [[float(i == j) for j in range(n_rows)] for i in range(n_rows)]
	Q = [[float(i == j) for j in range(n_rows)] for i in range(n_rows)]
	for j in range(n_cols): #0, 1, 2, 3, 4
		r = [0]*n_rows
		for i in range(j+1, n_rows+1):		
			#pivoteo
			if i < n_rows:
				a, P, Q = pivot(a, P, Q, j, piv)

			if i == j+1:			
			#Eliminacion de gauss			
				for m in range(j):
					r[m] = a[m][j] / float(a[j][j])

				for l in range(j):
					for k in range(j, n_cols):
						if l < n_rows:
							a[l][k] = a[l][k] - r[l]*a[j][k]

			if i < n_rows:		
				r[i] = a[i][j] / float(a[j][j])				
				for k in range(j, n_cols):
					a[i][k] = a[i][k] - r[i]*a[j][k]				

		if j < n_rows:
			temp = [0]*n_rows
			if (temp != r):
				print 'T%d' % (j+1)
				printMatrix(T(n_rows, r, j))			

	#para convertir la diagonal en 1's.
	for i in range(n_rows):
		temp = float(a[i][i])
		for j in range(n_cols):
			a[i][j] = a[i][j] / temp

	print 'A|b:'			
	printMatrix(a)
	if piv == 2:
		return multi(Q, reverseSub(a))	

	return reverseSub(a)	

def croutLU1(A, piv=0):
	"""Calcula el metodo LU1"""

	n = len(A)
	L = [[0.0]*n for j in range(n)];
	U = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(k, n):
			L[i][k] = A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))
		for i in range(k+1, n):
			if L[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			U[k][i] = (A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))) / float(L[k][k])	

	if piv == 2:
		return P, Q, L, U	
	return P, L, U

def croutL1U(A, piv=0):
	"""Calcula el metodo L1U"""

	n = len(A)
	U = [[0.0]*n for j in range(n)];
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(k, n):
			U[k][i] = A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))
		for i in range(k+1, n):
			if U[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			L[i][k] = (A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))) / float(U[k][k])	
	
	if piv == 2:
		return P, Q, L, U

	return P, L, U

def Doolittle(A, piv=0):

	n = len(A)
	U = [[0.0]*n for j in range(n)]
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(n):
			U[k][i] = A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))
		for i in range(k, n):
			if U[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			L[i][k] = (A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))) / float(U[k][k])	

	if piv == 2:
		return P, Q, L, U
	return P, L, U

def LDMt(A, piv=0):
	""""""
	if piv == 2:
		P, I, L, U = croutLU1(A, piv)
	else:
		P, L, U = croutLU1(A, piv)

	D, M = diagonalDMt(U)
	return L, D, M

def diagonalDMt(U):

	n = len(U)
	D = [[float(i == j) for j in range(n)] for i in range(n)]
	for i in range(n):
		D[i][i] = float(U[i][i])
	for i in range(n):
		for j in range(n):
			U[i][j] /= (D[i][i])
	return D, U

def cholesky(A):
	if not symmetricMatrix(A):
		exit('La matriz no es simetrica')
	n = len(A)
	G = [[0.0]*n for j in range(n)]
	for i in range(n):
		suma = A[i][i]
		for k in range(i):
			suma -= A[k][i]**2
		if suma < 0:
			exit('No es definida positiva')	
		A[i][i] = sqrt(suma)
		for j in range(i+1, n):
			suma = A[i][j]
			for k in range(i):
				suma -= A[k][i]*A[k][j]
			A[i][j] = suma / A[i][i]

	for j in range(n):
		for i in range(n):
			if(i > j):
				A[i][j] = 0.0
	return A	

def norm(x):
    """Norma euclideana del vector x, ie, ||x||2"""

    return sqrt(sum([x_i**2 for x_i in x]))

def Q_i(Q_min, i, j, k):
    """Rellenando elementos de la matrix Q_t"""

    if i < k or j < k:
        return float(i == j)
    else:
        return Q_min[i-k][j-k]

def Householder(A):
	""" Transformacion por Householder, donde Hn...H2*H1*A = R, Q_t*A = R
		Q = H1*H2 ... Hn, Q es ortogonal (Q * Qt = I)
		Retorna Q y R, donde A = Q*R 
	"""

	n = len(A)

	R = A
	Q = [[0.0] * n for i in range(n)]

	for k in range(n-1):                                                                    
		I = [[float(i == j) for i in range(n)] for j in range(n)]

		# Se crea los vectores x, e y un escalar alpha
		x = [row[k] for row in R[k:]]
		e = [row[k] for row in I[k:]]
		#La funcion cmp(a, b) retorna -1 si a<b, 1 si a>b, 0 si a==b
		alpha = -cmp(x[0],0) * norm(x) 

		#Se crea los vectores u, v
		u = map(lambda p,q: p + alpha * q, x, e)
		norm_u = norm(u)
		v = map(lambda p: p/norm_u, u)

		#Se crea la matriz menor Q_t
		Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in range(n-k)] for j in range(n-k) ]

		#Se rellena la matriz menor Q (Q_min)
		Q_t = [[ Q_i(Q_min,i,j,k) for i in range(n)] for j in range(n)]

		#Si esta en la primera ejecutada, entonces se calcula Q_t*A = R
		#Sino, Q_tn .. Q_t1*A = R
		if k == 0:
			Q = Q_t
			R = matrixMulti(Q_t,A)
		else:
			Q = matrixMulti(Q_t,Q)
			R = matrixMulti(Q_t,R)

	#Se retorna la transpuesta de los Q_t, ie, trans(Q_tn* ... *Q_t1) = Q
	return trans(Q), R

def givens(A):
	"""	Gn* ... G2*G1*A = R
		Q_t = Gn* ... G2*G1
		A = Q*R, de la propiedad Q_t * Q = I
	"""
	n = len(A)

	An = A	
	Gn = [[float(i == j) for j in range(n)] for i in range(n)]
	Q_t = [[float(i == j) for j in range(n)] for i in range(n)]

	a = An[0][n-2]
	b = An[0][n-1]
	index = 1
	for i in range(n):
		for j in range(n-1, i, -1):
			a = An[j-1][i]
			b = An[j][i]
			if a*a + b*b == 0:
				continue

			cosX = a / (sqrt(a*a + b*b)) 
			sinX = -b / (sqrt(a*a + b*b))

			Gn[j][j] = cosX
			Gn[j][j-1] = sinX
			Gn[j-1][j] = -sinX
			Gn[j-1][j-1] = cosX

			print 'G' +str(index) + ':'
			printMatrix(Gn)

			An = matrixMulti(Gn, An)

			print 'A' +str(index) + ':'
			printMatrix(An)

			Q_t = matrixMulti(Gn, Q_t)
			#Volviendo la matriz Gn a la identidad
			Gn = [[float(k == l) for l in range(n)] for k in range(n)]
			index += 1
	return trans(Q_t), An

def normaInfVector(L):
	""" Calcula la norma infinita de un vector:
		||x|| = max {|xi|}, i = 0, 1, ... n.
	"""

	maximum = fabs(L[0])
	for i in range(1, len(L)):
		maximum = max(maximum, fabs(L[i]))
	return maximum	

def jacobi(A, b, prec=1e-7):
	"""Metodo que calcula la solucion Ax = b, usando tecnicas iterativas"""
	
	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')

	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	for j in range(n):
		dominancia = 0.0
		for i in range(j+1, n):
			if j != i:
				dominancia += fabs(A[i][j])
		if A[i][i] < dominancia:
			exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation = sum(A[i][j]*Xk1[j] if i!=j else 0 for j in range(n))
			Xk1[i] = (1.0/A[i][i])*(b[i] - sumation)

		print Xk1
	return Xk1			

def gaussSeidel(A, b, prec=1e-7):
	"""Metodo que calcula la solucion Ax = b, usando tecnicas iterativas"""

	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')
	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	for j in range(n):
		dominancia = 0.0
		for i in range(j+1, n):
			if j != i:
				dominancia += fabs(A[i][j])
		if A[i][i] < dominancia:
			exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation1 = sum(A[i][j]*Xk1[j] for j in range(i))
			sumation2 = sum(A[i][j]*Xk1[j] for j in range(i+1, n))

			Xk1[i] = (1.0/A[i][i])*(b[i] - sumation1 - sumation2)

		print Xk1
	return Xk1

def sor(A, b, prec=1e-7, w=1.5):
	""" Metodo SOR, calcula la solucion de un sistema de ecuaciones
		usando el parametro de relajacion w, que necesariamente varia
		de <0, 2> solo en el caso que la matriz converja.
	"""
	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	#Garantizar que los elementos de la diagonal principal sean distintos de cero
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')

	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	# for j in range(n):
	# 	dominancia = 0.0
	# 	for i in range(n):
	# 		if j != i:
	# 			dominancia += fabs(A[i][j])
	# 	if A[i][i] < dominancia:
	# 		exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation1 = sum(A[i][j]*Xk1[j] for j in range(i))
			sumation2 = sum(A[i][j]*Xk1[j] for j in range(i+1, n))

			Xk1[i] = (float(w)/A[i][i])*(b[i] - sumation1 - sumation2) + (1-w)*Xk[i]

#		print Xk1
	return Xk1

def menu():
	print "\t\t\t  Bienvenidos al Programa de"
	print "\t\t\tResolucion de Ecuaciones Lineales"
	opc = 1
	while(opc):
		print "Elija el metodo a usar:\n"
		print """
			0.  Definiciones
			1.  Gauss
			2.  Gauss-Jordan
			3.  Crout LU1
			4.  Crout L1U
			5.  Doolittle
			6.  LDMt
			7.  Cholesky
			8.  <NAME> (falta)
			9.  Aasen (falta)
			10. Householder
			11. Givens
			12. Jacobi
			13. Gauss-Seidel
			14. SOR
			15. Condicionamiento (incompleto)
			16. Salir
			  """
		opc = input('Ingrese opcion\n')
		if type(opc) != type(1):
			exit('Debe ingresar un valor entero, vuela a intentarlo')
		if opc == 16:
			exit('Gracias por usar el programa')
		if opc < 0 or opc > 16:
			print 'Opcion incorrecta'
			continue
		#Ver la documentacion	
		if opc == 0:
			print __doc__
			opc = 1
			continue

		Matrix = inputMatrix()		 		
		#Matrix = [
		#		[1, -4, 2, 1],
		 #		[2, -6, 1, 4],
		 #		[-1, 2, 3, 4],
		 #		[0, -1, 1, 1]
		 #		]		 		
		#El condicionamiento no necesita el vector b
		if opc != 15:
			
			b = inputVector(Matrix)
	 	 
			#b = [-4, 1, 12, 0]	 	 
			restric = [7, 10, 11, 12, 13, 14]
			if opc not in restric:
				print "Elija el pivoteo a usar:\n"
				print """
					0.  Sin pivoteo
					1.  Pivoteo parcial
					2.  Pivoteo total
						
					  """
				piv = int(raw_input('Ingrese opcion\n'))
		#Gauss
		if opc == 1:
			x = gauss(Matrix, b, piv)
			print 'x:'
			print x

		#Gauss-Jordan	
		elif opc == 2:
			x = gaussJordan(Matrix, b, piv)
			print 'x:'
			print x

		#croutLU1, croutL1U, Doolittle
		elif opc == 3 or opc == 4 or opc == 5:
			#croutLU1
			if opc == 3:
				if piv == 2:
					P, Q, L, U = croutLU1(Matrix, piv)
				else:		
					P, L, U = croutLU1(Matrix, piv)
			#croutL1U		
			elif opc == 4:
				if piv == 2:
					P, Q, L, U = croutL1U(Matrix, piv)
				else:		
					P, L, U = croutL1U(Matrix, piv)
			#Doolittle		
			elif opc == 5:
				if piv == 2:
					P, Q, L, U = Doolittle(Matrix, piv)
				else:		
					P, L, U = Doolittle(Matrix, piv)			

			#Imprimir las matrices	
			print 'P:'
			printMatrix(P)
			if piv == 2:
				print 'Q:'
				printMatrix(Q)		
			print 'L:'
			printMatrix(L)
			print 'U:'
			printMatrix(U)
		
			#calcula las soluciones de Lz=b y Ux=z	
			if piv == 2:
				z = solMatrixInf(L, b)
				y = solMatrixSup(U,z)
				x = multi(P, y)
			else:
				z = solMatrixInf(L, b)
				x = solMatrixSup(U,z) 	

			#imprime las soluciones
			print "z:"
			print z
			print "x:"
			print x

		#LDMt
		elif opc == 6:
			L,D,M = LDMt(Matrix, piv)
			print 'L:'
			printMatrix(L)		
			print "D:"
			printMatrix(D)
			print "M:"
			printMatrix(M)

			#calcula las soluciones
			y = solMatrixInf(L, b)
			z = solMatrixInf(D, y)
			x = solMatrixSup(M, z)
			
			#imprime las soluciones
			print "y:"
			print y
			print "z:"
			print z
			print "x:"
			print x

		#Cholesky
		elif opc == 7:
			G = cholesky(Matrix)
			Gt = trans(G) 
			print "G:"
			printMatrix(G)
			print "Gt:"
			printMatrix(Gt)

			#calcula las soluciones
			z = solMatrixInf(Gt, b)
			x = solMatrixSup(G, z)
			
			#imprime las soluciones
			print "z:"
			print z
			print "x:"
			print x

		#Parlett y Reid	
		elif opc == 8:
			print 'Falta implementar'	
		
		#Aasen	
		elif opc == 9:
			print 'Falta implementar'

		#Householder y Givens	
		elif opc == 10 or opc == 11:
			#Householder
			if opc == 10:	
				Q, R = Householder(Matrix)
			#Givens
			elif opc == 11:
				Q, R = givens(Matrix)

			print 'Q:'
			printMatrix(Q)			  

			print 'R:'
			printMatrix(R)			  
			#Comprobar la igualdad
			print 'A = Q*R :'
			printMatrix(matrixMulti(Q, R))

			b_prima = multi(trans(Q), b)
			print 'b\':'
			print b_prima

			x = solMatrixSup(R, b_prima)
			print 'x:'
			print x

			print 'Ax = b :'
			print multi(Matrix, x)
		
		#Jacobi
		elif opc == 12:
			prec = float(raw_input('Ingrese la precision, ejemplo : 0.001  '))
			x = jacobi(Matrix, b, prec)
			print 'x:'
			print x

			print 'Ax = b :'
			print multi(Matrix, x)

		#Gauss-Seidel	
		elif opc == 13:
			prec = float(raw_input('Ingrese la precision, ejemplo : 0.001  '))
			x = gaussSeidel(Matrix, b, prec)
			print 'x:'
			print x
		
		#SOR (Successive Over relaxation)
		elif opc == 14:
			prec = float(raw_input('Ingrese la precision, ejemplo: 0.001 : '))
			w = float(raw_input('Ingrese el parametro de relajacion <0, 2>, ejemplo 0.5 : '))
			x = sor(Matrix, b, prec, w)
			print 'x:'
			print x
		
		#Condicionamiento	
		elif opc == 15:
			
			while 1:
				print "Elija el tipo de condicionamiento:\n"
				print """
					1.  Norma 1
					2.  Norma 2
					3.  Norma infinita
					4.  Condicionamiento con norma 1
					5.  Condicionamiento con norma 2
					6.  Condicionamiento con norma infinita
					7.  Salir
						
					  """
				nor = int(raw_input('Ingrese Opcion\n'))
				#Norma 1
				if nor == 1:
					print 'norm1(M):'
					print norm1(Matrix)
				
				#Norma 2	
				elif nor == 2:
					print 'Falta implementar'
					# print norma2(Matrix)	
				
				#Norma infinita	
				elif nor == 3:
					print 'infinityNorm(M):'
					print infinityNorm(Matrix)	

				#Condicionamiento de la matriz usando norma 1	
				elif nor == 4:
					print 'Falta implementar'	

				#Condicionamiento de la matriz usando norma 2	
				elif nor == 5:
					print 'Falta implementar'	

				#Condicionamiento de la matriz usando norma infinita	
				elif nor == 6:
					print 'Falta implementar'	

				#Salir del menu	
				elif nor == 7:
					exit('Gracias por usar el programa')
				else:
					print 'Opcion incorrecta'
					continue
				print 'Resultado obtenido con exito'		 		

	print 'Resultado obtenido con exito'	

#menu()

"""Paquete que muestra la resolucion y factorizacion de un sistema de ecuaciones lineales
	 aplicando diferentes metodos.

   Nota: 
	- La mayoria de las funciones y las variables estan en ingles, y los comentarios en espanol.	
	- Aun falta implementar varios metodos, y mejorar los existentes.

        - El paquete no soporta matrices singulares.

   Definiciones:
   A, M : Matriz original.
   a : Matriz aumentada (A|b)
   b : Vector independiente. 
   L : Matriz triangular inferior unitaria.
   M : Matriz triangular inferior unitaria.
   U : Matriz triangular superior unitaria.
   D : Matriz diagonal.
   P : Matriz de permutacion
   Pr : Matriz de permutacion de filas
   Pc : Matriz de permutacion de columnas
   G : Matriz superior hallada con el metodo de Cholesky.
   Gt : Transpuesta de G.
   A = LU : significa factorizacion sin pivoteo
   PA = LU : significa factorizacion con pivoteo parcial
   PAQ = LU : significa factorizacion con pivoteo total

"""
from math import fabs
from math import sqrt
from sys import exit


def inputMatrix():
	"""Ingreso de datos de la matriz"""

	n = int(input("Ingrese el orden de la matriz: "))
	print "Ingrese los elementos de la matriz A fila por fila con un espacio luego enter"
        print "Por ejemplo, para una matriz de orden 3:","\n", "3 2 1"
        print "1 7 5","\n\n"
	A = [[0.0]*n for i in range(n)]
	aux = [[0.0]*n for i in range(n)]
	for i in range(n):
		temp = raw_input()
		A[i] = temp.split()
		for j in range(n):
			A[i][j] = float(A[i][j])
	if A == aux:
		exit('Matriz nula, vuelva a escribir la matriz')
	return A		

def inputVector(M):
	"""Ingreso de datos del vector independiente"""

	n = len(M)
	print 'Ingrese los elementos del vector b'
	b = [0]*n
	#for i in range(n):
	#	b[i] = float(input())

	temp = raw_input();
	b = temp.split()
	for i in range(n):
		b[i] = float(b[i])

	return b	
	
def printMatrixD(M):
	"""Imprime la matriz, version debugger."""

	for i in range(len(M)):
		# print '|',
		for j in range(len(M[i])):
			if(i == 0 and j == 0):
				print " ",
				for x in range(len(M[i])):
					print '{0:8}'.format(x),
				print
				print
			if(j == 0):
				print i,

			if(j == len(M)):
				print '|',
				print '{0:8.4f}'.format(M[i][j]),
			else:
				print '{0:8.4f}'.format(M[i][j]),	
		print '|'
	print

def printMatrix(M):
	"""Imprime la matriz de una forma legible."""

	for i in range(len(M)):
		print '|',
		for j in range(len(M[i])):
			if(j == len(M)):
				print '|',
				print '{0:8.4f}'.format(M[i][j]),
			else:
				print '{0:8.4f}'.format(M[i][j]),	
		print '|'
	print  

def norm1(A):
	"""Calcula la norma 1"""

	summation = 0
	for i in range(len(A)):
		summation += abs(A[i][0])

	for j in range(1, len(A)):
		temp = 0
		for i in range(len(A[i])):
			temp += abs(A[i][j])
		summation = temp if (temp > summation) else summation

	return summation 

def infinityNorm(A):
	"""Calcula la norma infinita"""

	summation = 0
	for j in range(len(A)):
		summation += abs(A[0][j])

	for i in range(1, len(A)):
		temp = 0
		for j in range(len(A[i])):
			temp += abs(A[i][j])
		summation = temp if (temp > summation) else summation
			
	return summation  

def symmetricMatrix(A):
	"""Calcula la matriz simetrica de A"""

	for i in range(len(A)):
		for j in range(i+1,len(A)):
			if A[i][j] != A[j][i]:
				return False
	return True 

def trans(M):
	"""Calcula la matriz transpuesta de M"""

	n = len(M)
	return [[ M[i][j] for i in range(n)] for j in range(n)]

def reverseSub(a):
	"""Se usa para resolver la matriz superior de gauss/gauss-Jordan"""

	n = len(a)
	x = [0]*n
	for j in range(n-1, -1, -1):
            if a[j][j] != 0:#Una manera de evitar una excepcion en matriz singular
                    x[j] = (a[j][n] - sum(a[j][k]*x[k] for k in range(j+1, n)))/float(a[j][j])
            else:
                print 'Matriz singular, la ecuacion no tiene solucion o tiene infinitas soluciones .'
                break
	return x	

def solMatrixSup(M, b):
	"""Resuelve la matriz superior con el metodo inverso"""	

	x = []
	for i in range(len(M)-1,-1,-1):
		x.append((1.0/(M[i][i]))*(b[i]-sum(M[i][len(M)-j-1]*x[j] for j in range(len(x)))))
	x.reverse()
	
	return x

def solMatrixInf(M, b):
	"""Resuelve la matriz inferior con el metodo inverso"""

	x = []
	for i in range(len(M)):
		x.append((1.0/(M[i][i]))*(b[i]-sum(M[i][j]*x[j] for j in range(len(x)))))
	return x
	
def augmentedMatrix(M, b):
	"""Retorna la matriz aumentada"""
	
	a = M
	for i in range(len(M)):
		a[i].append(b[i])

	return a	

def maxColum(M, c):
	"""Devuelve la fila del mayor valor debajo de la diagonal en la columna c"""

	r = c #fila
	maximum = M[c][c]
	for i in range(c+1,len(M)):
		if(fabs(maximum) < fabs(M[i][c])):
			maximum = M[i][c]
			r = i
	return r

def maxSubMatrix(M, c):
	"""Retorna el mayor elemento de la submatriz A[c]"""

	row = c
	colum = c
	n = len(M)
	maximum = M[c][c]
		
	for j in range(c, n):
		for k in range(c, n):
			maxTemp = M[k][j]
			if(fabs(maximum) < fabs(maxTemp)):
				maximum = maxTemp
				colum = j
				row = k

	return row, colum			

def exchangeRows(M, r1, r2):
	"""Intercambia las filas r1 y r2 de M"""

	M[r1], M[r2] = M[r2], M[r1]
	return M

def exchangeCols(M, c1, c2):
	"""Intercambia las columnas c1 y c2 de M"""

	for k in range(len(M)):
		M[k][c1] , M[k][c2] = M[k][c2], M[k][c1]
	return M	

def pivot(a, P, Q, colum, piv=0):
	"""Se encarga del pivoteo en cualquier metodo."""

	if piv > 2 or piv < 0:
		exit('Valores invalidos para el parametro pivoteo, valores validos: 0, 1, 2.')
	n = len(a)
		
	temp = a[colum][colum]
	if(temp == 0.0 and piv == 0):
		row_maxColumn = maxColum(a, colum)
		a = pivotP(a, row_maxColumn, colum)
		P = exchangeRows(P, row_maxColumn, colum)
		print 'P(%d,%d)' % (row_maxColumn, colum)
		printMatrix(Pr(n, row_maxColumn, colum))
		
	elif piv == 1:
		row_maxColumn = maxColum(a, colum)
		if row_maxColumn != colum:
			a = pivotP(a, row_maxColumn, colum)
			P = exchangeRows(P, row_maxColumn, colum)
			print 'P(%d,%d)' % (row_maxColumn, colum)
			printMatrix(Pr(n, row_maxColumn, colum))
			
	elif piv == 2:
		row, c = maxSubMatrix(a, colum)	
		if (row != colum) or (c != colum) :
			a = pivotT(a, colum)
			P = exchangeCols(P, row, colum)
			Q = exchangeCols(Q, colum, c)
			print 'P(%d,%d):' % (colum, row)
			printMatrix(Pr(n, row, colum))
			print 'Q(%d,%d):' % (colum, c)
			printMatrix(Pr(n, c, colum))
	
	#retorna la matriz, la matriz de permutacion de filas y la de columnas.	
	return 	a, P, Q		

def pivotP(M, r1, r2):
	"""Permuta la fila r1 con la fila r2 de la matriz M"""

	return exchangeRows(M, r1, r2)

def pivotT(M, i):
	"""Busca el mayor elemento de la submatriz A[i] y permuta filas y columnas"""

	r, c = maxSubMatrix(M, i)
	M = pivotP(M, i, r)		

	return exchangeCols(M, c, i)

def L(n, ri, c):
	"""Calcula las matrices Li.

		ri: a[i][c] / a[c][c]
		c : columna c-esima 
	"""
	#Matriz identidad
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	for k in range(c+1, n):
		L[k][c] = -ri[k]
			
	return L

def matrixMulti(A, B):
	"""Multiplica dos matrices, C = A*B """

	rowsA, colsA = len(A), len(A[0])
	rowsB, colsB = len(B), len(B[0])

	if colsA != rowsB:
		exit('Dimensiones incorrectas')

	C = [[0 for row in range(colsB)] for col in range(rowsA)]

	for i in range(rowsA):
		for j in range(colsB):
			for k in range(colsA):
				C[i][j] += A[i][k]*B[k][j]
	return C

def T(n, ri, c):
	"""Calcula las matrices Ti (Para el metodo Gauss-Jordan)

		ri: a[i][c] / a[c][c]
		c : columna c-esima 
	"""
	T = [[float(i == j) for j in range(n)] for i in range(n)]
	for k in range(n):
		T[k][c] = -ri[k]

	return T

def Pr(n, r1, r2):
	"""Calcula la matriz de permutacion de fila"""

	#Matriz identidad
	I = [[float(i == j) for j in range(n)] for i in range(n)]
	return exchangeRows(I, r1, r2)

def Pc(n, c1, c2):
	"""Calcula la matriz de permutacion de columna"""

	#Matriz identidad
	I = [[float(i == j) for j in range(n)] for i in range(n)]
	return exchangeCols(I, c1, c2)	

def multi(A,x):
	"""Multiplica una matriz con un vector"""

	c=[]
	for i in range(len(A)):
		c.append(sum(A[i][j]*x[j] for j in range(len(A))))
	return c

def invLU(A):
	P, L, U = croutLU1(A)
	for x in range(1,10):
		pass
		
def gauss(A, b, piv=0):
	"""Metodo de Gauss

		A: Matriz
		b: vector independiente
		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial
			 piv=2 con pivoteo total

	"""
	a = augmentedMatrix(A, b)
	n = len(A)
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]
	for j in range(len(a[0])):
		r = [0]*n
		# pivoteo
		if j < n:
			a, P, Q = pivot(a, P, Q, j, piv)

		for i in range(j+1, n):		
			
			#Eliminacion de gauss
			r[i] = a[i][j] / float(a[j][j])

			for k in range(j, len(a[0])):
				a[i][k] = a[i][k] - r[i]*a[j][k]
					
		if j < n:
			temp = [0]*n
			if (temp != r):
				print ''#'L%d' % (j+1)
#				printMatrix(L(n, r, j))			

#	print 'A|b:'			
#	printMatrix(a)

	if piv == 2:
		return multi(Q, reverseSub(a))	

	return reverseSub(a)			

def gaussJordan(A, b, piv=0):
	"""Metodo de Gauss-Jordan
		A: Matriz
		b: vector independiente
		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial
			 piv=2 con pivoteo total  	
	"""
	a = augmentedMatrix(A, b)
	n_rows = len(a)
	n_cols = len(a[0])
	P = [[float(i == j) for j in range(n_rows)] for i in range(n_rows)]
	Q = [[float(i == j) for j in range(n_rows)] for i in range(n_rows)]
	for j in range(n_cols): #0, 1, 2, 3, 4
		r = [0]*n_rows
		for i in range(j+1, n_rows+1):		
			#pivoteo
			if i < n_rows:
				a, P, Q = pivot(a, P, Q, j, piv)

			if i == j+1:			
			#Eliminacion de gauss			
				for m in range(j):
					r[m] = a[m][j] / float(a[j][j])

				for l in range(j):
					for k in range(j, n_cols):
						if l < n_rows:
							a[l][k] = a[l][k] - r[l]*a[j][k]

			if i < n_rows:		
				r[i] = a[i][j] / float(a[j][j])				
				for k in range(j, n_cols):
					a[i][k] = a[i][k] - r[i]*a[j][k]				

		if j < n_rows:
			temp = [0]*n_rows
			if (temp != r):
				print 'T%d' % (j+1)
				printMatrix(T(n_rows, r, j))			

	#para convertir la diagonal en 1's.
	for i in range(n_rows):
		temp = float(a[i][i])
		for j in range(n_cols):
			a[i][j] = a[i][j] / temp

	print 'A|b:'			
	printMatrix(a)
	if piv == 2:
		return multi(Q, reverseSub(a))	

	return reverseSub(a)	

def croutLU1(A, piv=0):
	"""Calcula el metodo LU1"""

	n = len(A)
	L = [[0.0]*n for j in range(n)];
	U = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(k, n):
			L[i][k] = A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))
		for i in range(k+1, n):
			if L[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			U[k][i] = (A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))) / float(L[k][k])	

	if piv == 2:
		return P, Q, L, U	
	return P, L, U

def croutL1U(A, piv=0):
	"""Calcula el metodo L1U"""

	n = len(A)
	U = [[0.0]*n for j in range(n)];
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(k, n):
			U[k][i] = A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))
		for i in range(k+1, n):
			if U[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			L[i][k] = (A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))) / float(U[k][k])	
	
	if piv == 2:
		return P, Q, L, U

	return P, L, U

def Doolittle(A, piv=0):

	n = len(A)
	U = [[0.0]*n for j in range(n)]
	L = [[float(i == j) for j in range(n)] for i in range(n)]
	P = [[float(i == j) for j in range(n)] for i in range(n)]
	Q = [[float(i == j) for j in range(n)] for i in range(n)]

	#pivot
	if piv:
		for j in range(len(A)):
			A, P, Q = pivot(A, P, Q, j, piv)

	for k in range(n):
		for i in range(n):
			U[k][i] = A[k][i] - sum(L[k][p]*U[p][i] for p in range(k))
		for i in range(k, n):
			if U[k][k] == 0:
				exit('Debe usarse pivoteo parcial')
			L[i][k] = (A[i][k] - sum(L[i][p]*U[p][k] for p in range(k))) / float(U[k][k])	

	if piv == 2:
		return P, Q, L, U
	return P, L, U

def LDMt(A, piv=0):
	""""""
	if piv == 2:
		P, I, L, U = croutLU1(A, piv)
	else:
		P, L, U = croutLU1(A, piv)

	D, M = diagonalDMt(U)
	return L, D, M

def diagonalDMt(U):

	n = len(U)
	D = [[float(i == j) for j in range(n)] for i in range(n)]
	for i in range(n):
		D[i][i] = float(U[i][i])
	for i in range(n):
		for j in range(n):
			U[i][j] /= (D[i][i])
	return D, U

def cholesky(A):
	if not symmetricMatrix(A):
		exit('La matriz no es simetrica')
	n = len(A)
	G = [[0.0]*n for j in range(n)]
	for i in range(n):
		suma = A[i][i]
		for k in range(i):
			suma -= A[k][i]**2
		if suma < 0:
			exit('No es definida positiva')	
		A[i][i] = sqrt(suma)
		for j in range(i+1, n):
			suma = A[i][j]
			for k in range(i):
				suma -= A[k][i]*A[k][j]
			A[i][j] = suma / A[i][i]

	for j in range(n):
		for i in range(n):
			if(i > j):
				A[i][j] = 0.0
	return A	

def norm(x):
    """Norma euclideana del vector x, ie, ||x||2"""

    return sqrt(sum([x_i**2 for x_i in x]))

def Q_i(Q_min, i, j, k):
    """Rellenando elementos de la matrix Q_t"""

    if i < k or j < k:
        return float(i == j)
    else:
        return Q_min[i-k][j-k]

def Householder(A):
	""" Transformacion por Householder, donde Hn...H2*H1*A = R, Q_t*A = R
		Q = H1*H2 ... Hn, Q es ortogonal (Q * Qt = I)
		Retorna Q y R, donde A = Q*R 
	"""

	n = len(A)

	R = A
	Q = [[0.0] * n for i in range(n)]

	for k in range(n-1):                                                                    
		I = [[float(i == j) for i in range(n)] for j in range(n)]

		# Se crea los vectores x, e y un escalar alpha
		x = [row[k] for row in R[k:]]
		e = [row[k] for row in I[k:]]
		#La funcion cmp(a, b) retorna -1 si a<b, 1 si a>b, 0 si a==b
		alpha = -cmp(x[0],0) * norm(x) 

		#Se crea los vectores u, v
		u = map(lambda p,q: p + alpha * q, x, e)
		norm_u = norm(u)
		v = map(lambda p: p/norm_u, u)

		#Se crea la matriz menor Q_t
		Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in range(n-k)] for j in range(n-k) ]

		#Se rellena la matriz menor Q (Q_min)
		Q_t = [[ Q_i(Q_min,i,j,k) for i in range(n)] for j in range(n)]

		#Si esta en la primera ejecutada, entonces se calcula Q_t*A = R
		#Sino, Q_tn .. Q_t1*A = R
		if k == 0:
			Q = Q_t
			R = matrixMulti(Q_t,A)
		else:
			Q = matrixMulti(Q_t,Q)
			R = matrixMulti(Q_t,R)

	#Se retorna la transpuesta de los Q_t, ie, trans(Q_tn* ... *Q_t1) = Q
	return trans(Q), R

def givens(A):
	"""	Gn* ... G2*G1*A = R
		Q_t = Gn* ... G2*G1
		A = Q*R, de la propiedad Q_t * Q = I
	"""
	n = len(A)

	An = A	
	Gn = [[float(i == j) for j in range(n)] for i in range(n)]
	Q_t = [[float(i == j) for j in range(n)] for i in range(n)]

	a = An[0][n-2]
	b = An[0][n-1]
	index = 1
	for i in range(n):
		for j in range(n-1, i, -1):
			a = An[j-1][i]
			b = An[j][i]
			if a*a + b*b == 0:
				continue

			cosX = a / (sqrt(a*a + b*b)) 
			sinX = -b / (sqrt(a*a + b*b))

			Gn[j][j] = cosX
			Gn[j][j-1] = sinX
			Gn[j-1][j] = -sinX
			Gn[j-1][j-1] = cosX

			print 'G' +str(index) + ':'
			printMatrix(Gn)

			An = matrixMulti(Gn, An)

			print 'A' +str(index) + ':'
			printMatrix(An)

			Q_t = matrixMulti(Gn, Q_t)
			#Volviendo la matriz Gn a la identidad
			Gn = [[float(k == l) for l in range(n)] for k in range(n)]
			index += 1
	return trans(Q_t), An

def normaInfVector(L):
	""" Calcula la norma infinita de un vector:
		||x|| = max {|xi|}, i = 0, 1, ... n.
	"""

	maximum = fabs(L[0])
	for i in range(1, len(L)):
		maximum = max(maximum, fabs(L[i]))
	return maximum	

def jacobi(A, b, prec=1e-7):
	"""Metodo que calcula la solucion Ax = b, usando tecnicas iterativas"""
	
	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')

	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	for j in range(n):
		dominancia = 0.0
		for i in range(j+1, n):
			if j != i:
				dominancia += fabs(A[i][j])
		if A[i][i] < dominancia:
			exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation = sum(A[i][j]*Xk1[j] if i!=j else 0 for j in range(n))
			Xk1[i] = (1.0/A[i][i])*(b[i] - sumation)

		print Xk1
	return Xk1			

def gaussSeidel(A, b, prec=1e-7):
	"""Metodo que calcula la solucion Ax = b, usando tecnicas iterativas"""

	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')
	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	for j in range(n):
		dominancia = 0.0
		for i in range(j+1, n):
			if j != i:
				dominancia += fabs(A[i][j])
		if A[i][i] < dominancia:
			exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation1 = sum(A[i][j]*Xk1[j] for j in range(i))
			sumation2 = sum(A[i][j]*Xk1[j] for j in range(i+1, n))

			Xk1[i] = (1.0/A[i][i])*(b[i] - sumation1 - sumation2)

		print Xk1
	return Xk1

def sor(A, b, prec=1e-7, w=1.5):
	""" Metodo SOR, calcula la solucion de un sistema de ecuaciones
		usando el parametro de relajacion w, que necesariamente varia
		de <0, 2> solo en el caso que la matriz converja.
	"""
	n = len(A)
	Xk = [0.0]*n
	sumation = 0.0
	#Garantizar que los elementos de la diagonal principal sean distintos de cero
	for i in range(n):
		if A[i][i] == 0:
			exit('Los elementos A[i][i] deben ser diferentes de 0')

	Xk1 = [b[i]/float(A[i][i]) for i in range(n)]
	minus = lambda x, y: [x[i]-y[i] for i in range(n)]

	# for j in range(n):
	# 	dominancia = 0.0
	# 	for i in range(n):
	# 		if j != i:
	# 			dominancia += fabs(A[i][j])
	# 	if A[i][i] < dominancia:
	# 		exit('La matriz no converge')

	while (normaInfVector(minus(Xk1,Xk)) / float(normaInfVector(Xk1))) > prec:
		Xk[:] = Xk1[:]
		for i in range(n):
			sumation1 = sum(A[i][j]*Xk1[j] for j in range(i))
			sumation2 = sum(A[i][j]*Xk1[j] for j in range(i+1, n))

			Xk1[i] = (float(w)/A[i][i])*(b[i] - sumation1 - sumation2) + (1-w)*Xk[i]

#		print Xk1
	return Xk1

def menu():
	print "\t\t\t  Bienvenidos al Programa de"
	print "\t\t\tResolucion de Ecuaciones Lineales"
	opc = 1
	while(opc):
		print "Elija el metodo a usar:\n"
		print """
			0.  Definiciones
			1.  Gauss
			2.  Gauss-Jordan
			3.  Crout LU1
			4.  Crout L1U
			5.  Doolittle
			6.  LDMt
			7.  Cholesky
			8.  <NAME> (falta)
			9.  Aasen (falta)
			10. Householder
			11. Givens
			12. Jacobi
			13. Gauss-Seidel
			14. SOR
			15. Condicionamiento (incompleto)
			16. Salir
			  """
		opc = input('Ingrese opcion\n')
		if type(opc) != type(1):
			exit('Debe ingresar un valor entero, vuela a intentarlo')
		if opc == 16:
			exit('Gracias por usar el programa')
		if opc < 0 or opc > 16:
			print 'Opcion incorrecta'
			continue
		#Ver la documentacion	
		if opc == 0:
			print __doc__
			opc = 1
			continue

		Matrix = inputMatrix()		 		
		#Matrix = [
		#		[1, -4, 2, 1],
		 #		[2, -6, 1, 4],
		 #		[-1, 2, 3, 4],
		 #		[0, -1, 1, 1]
		 #		]		 		
		#El condicionamiento no necesita el vector b
		if opc != 15:
			
			b = inputVector(Matrix)
	 	 
			#b = [-4, 1, 12, 0]	 	 
			restric = [7, 10, 11, 12, 13, 14]
			if opc not in restric:
				print "Elija el pivoteo a usar:\n"
				print """
					0.  Sin pivoteo
					1.  Pivoteo parcial
					2.  Pivoteo total
						
					  """
				piv = int(raw_input('Ingrese opcion\n'))
		#Gauss
		if opc == 1:
			x = gauss(Matrix, b, piv)
			print 'x:'
			print x

		#Gauss-Jordan	
		elif opc == 2:
			x = gaussJordan(Matrix, b, piv)
			print 'x:'
			print x

		#croutLU1, croutL1U, Doolittle
		elif opc == 3 or opc == 4 or opc == 5:
			#croutLU1
			if opc == 3:
				if piv == 2:
					P, Q, L, U = croutLU1(Matrix, piv)
				else:		
					P, L, U = croutLU1(Matrix, piv)
			#croutL1U		
			elif opc == 4:
				if piv == 2:
					P, Q, L, U = croutL1U(Matrix, piv)
				else:		
					P, L, U = croutL1U(Matrix, piv)
			#Doolittle		
			elif opc == 5:
				if piv == 2:
					P, Q, L, U = Doolittle(Matrix, piv)
				else:		
					P, L, U = Doolittle(Matrix, piv)			

			#Imprimir las matrices	
			print 'P:'
			printMatrix(P)
			if piv == 2:
				print 'Q:'
				printMatrix(Q)		
			print 'L:'
			printMatrix(L)
			print 'U:'
			printMatrix(U)
		
			#calcula las soluciones de Lz=b y Ux=z	
			if piv == 2:
				z = solMatrixInf(L, b)
				y = solMatrixSup(U,z)
				x = multi(P, y)
			else:
				z = solMatrixInf(L, b)
				x = solMatrixSup(U,z) 	

			#imprime las soluciones
			print "z:"
			print z
			print "x:"
			print x

		#LDMt
		elif opc == 6:
			L,D,M = LDMt(Matrix, piv)
			print 'L:'
			printMatrix(L)		
			print "D:"
			printMatrix(D)
			print "M:"
			printMatrix(M)

			#calcula las soluciones
			y = solMatrixInf(L, b)
			z = solMatrixInf(D, y)
			x = solMatrixSup(M, z)
			
			#imprime las soluciones
			print "y:"
			print y
			print "z:"
			print z
			print "x:"
			print x

		#Cholesky
		elif opc == 7:
			G = cholesky(Matrix)
			Gt = trans(G) 
			print "G:"
			printMatrix(G)
			print "Gt:"
			printMatrix(Gt)

			#calcula las soluciones
			z = solMatrixInf(Gt, b)
			x = solMatrixSup(G, z)
			
			#imprime las soluciones
			print "z:"
			print z
			print "x:"
			print x

		#Parlett y Reid	
		elif opc == 8:
			print 'Falta implementar'	
		
		#Aasen	
		elif opc == 9:
			print 'Falta implementar'

		#Householder y Givens	
		elif opc == 10 or opc == 11:
			#Householder
			if opc == 10:	
				Q, R = Householder(Matrix)
			#Givens
			elif opc == 11:
				Q, R = givens(Matrix)

			print 'Q:'
			printMatrix(Q)			  

			print 'R:'
			printMatrix(R)			  
			#Comprobar la igualdad
			print 'A = Q*R :'
			printMatrix(matrixMulti(Q, R))

			b_prima = multi(trans(Q), b)
			print 'b\':'
			print b_prima

			x = solMatrixSup(R, b_prima)
			print 'x:'
			print x

			print 'Ax = b :'
			print multi(Matrix, x)
		
		#Jacobi
		elif opc == 12:
			prec = float(raw_input('Ingrese la precision, ejemplo : 0.001  '))
			x = jacobi(Matrix, b, prec)
			print 'x:'
			print x

			print 'Ax = b :'
			print multi(Matrix, x)

		#Gauss-Seidel	
		elif opc == 13:
			prec = float(raw_input('Ingrese la precision, ejemplo : 0.001  '))
			x = gaussSeidel(Matrix, b, prec)
			print 'x:'
			print x
		
		#SOR (Successive Over relaxation)
		elif opc == 14:
			prec = float(raw_input('Ingrese la precision, ejemplo: 0.001 : '))
			w = float(raw_input('Ingrese el parametro de relajacion <0, 2>, ejemplo 0.5 : '))
			x = sor(Matrix, b, prec, w)
			print 'x:'
			print x
		
		#Condicionamiento	
		elif opc == 15:
			
			while 1:
				print "Elija el tipo de condicionamiento:\n"
				print """
					1.  Norma 1
					2.  Norma 2
					3.  Norma infinita
					4.  Condicionamiento con norma 1
					5.  Condicionamiento con norma 2
					6.  Condicionamiento con norma infinita
					7.  Salir
						
					  """
				nor = int(raw_input('Ingrese Opcion\n'))
				#Norma 1
				if nor == 1:
					print 'norm1(M):'
					print norm1(Matrix)
				
				#Norma 2	
				elif nor == 2:
					print 'Falta implementar'
					# print norma2(Matrix)	
				
				#Norma infinita	
				elif nor == 3:
					print 'infinityNorm(M):'
					print infinityNorm(Matrix)	

				#Condicionamiento de la matriz usando norma 1	
				elif nor == 4:
					print 'Falta implementar'	

				#Condicionamiento de la matriz usando norma 2	
				elif nor == 5:
					print 'Falta implementar'	

				#Condicionamiento de la matriz usando norma infinita	
				elif nor == 6:
					print 'Falta implementar'	

				#Salir del menu	
				elif nor == 7:
					exit('Gracias por usar el programa')
				else:
					print 'Opcion incorrecta'
					continue
				print 'Resultado obtenido con exito'		 		

	print 'Resultado obtenido con exito'	

#menu()
Paquete que muestra la resolucion y factorizacion de un sistema de ecuaciones lineales 	 aplicando diferentes metodos.     Nota:  	- La mayoria de las funciones y las variables estan en ingles, y los comentarios en espanol.	 	- Aun falta implementar varios metodos, y mejorar los existentes.          - El paquete no soporta matrices singulares.     Definiciones:    A, M : Matriz original.    a : Matriz aumentada (A|b)    b : Vector independiente.     L : Matriz triangular inferior unitaria.    M : Matriz triangular inferior unitaria.    U : Matriz triangular superior unitaria.    D : Matriz diagonal.    P : Matriz de permutacion    Pr : Matriz de permutacion de filas    Pc : Matriz de permutacion de columnas    G : Matriz superior hallada con el metodo de Cholesky.    Gt : Transpuesta de G.    A = LU : significa factorizacion sin pivoteo    PA = LU : significa factorizacion con pivoteo parcial    PAQ = LU : significa factorizacion con pivoteo total Ingreso de datos de la matriz Ingreso de datos del vector independiente #for i in range(n): #	b[i] = float(input()) Imprime la matriz, version debugger. # print '|', Imprime la matriz de una forma legible. Calcula la norma 1 Calcula la norma infinita Calcula la matriz simetrica de A Calcula la matriz transpuesta de M Se usa para resolver la matriz superior de gauss/gauss-Jordan #Una manera de evitar una excepcion en matriz singular Resuelve la matriz superior con el metodo inverso Resuelve la matriz inferior con el metodo inverso Retorna la matriz aumentada Devuelve la fila del mayor valor debajo de la diagonal en la columna c #fila Retorna el mayor elemento de la submatriz A[c] Intercambia las filas r1 y r2 de M Intercambia las columnas c1 y c2 de M Se encarga del pivoteo en cualquier metodo. #retorna la matriz, la matriz de permutacion de filas y la de columnas. Permuta la fila r1 con la fila r2 de la matriz M Busca el mayor elemento de la submatriz A[i] y permuta filas y columnas Calcula las matrices Li.  		ri: a[i][c] / a[c][c] 		c : columna c-esima #Matriz identidad Multiplica dos matrices, C = A*B Calcula las matrices Ti (Para el metodo Gauss-Jordan)  		ri: a[i][c] / a[c][c] 		c : columna c-esima Calcula la matriz de permutacion de fila #Matriz identidad Calcula la matriz de permutacion de columna #Matriz identidad Multiplica una matriz con un vector Metodo de Gauss  		A: Matriz 		b: vector independiente 		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial 			 piv=2 con pivoteo total # pivoteo #Eliminacion de gauss #				printMatrix(L(n, r, j)) #	print 'A|b:' #	printMatrix(a) Metodo de Gauss-Jordan 		A: Matriz 		b: vector independiente 		piv: pivoteo, piv=0 sin pivoteo, piv=1 con pivoteo parcial 			 piv=2 con pivoteo total #0, 1, 2, 3, 4 #pivoteo #Eliminacion de gauss #para convertir la diagonal en 1's. Calcula el metodo LU1 #pivot Calcula el metodo L1U #pivot #pivot Norma euclideana del vector x, ie, ||x||2 Rellenando elementos de la matrix Q_t Transformacion por Householder, donde Hn...H2*H1*A = R, Q_t*A = R 		Q = H1*H2 ... Hn, Q es ortogonal (Q * Qt = I) 		Retorna Q y R, donde A = Q*R # Se crea los vectores x, e y un escalar alpha #La funcion cmp(a, b) retorna -1 si a<b, 1 si a>b, 0 si a==b #Se crea los vectores u, v #Se crea la matriz menor Q_t #Se rellena la matriz menor Q (Q_min) #Si esta en la primera ejecutada, entonces se calcula Q_t*A = R #Sino, Q_tn .. Q_t1*A = R #Se retorna la transpuesta de los Q_t, ie, trans(Q_tn* ... *Q_t1) = Q Gn* ... G2*G1*A = R 		Q_t = Gn* ... G2*G1 		A = Q*R, de la propiedad Q_t * Q = I #Volviendo la matriz Gn a la identidad Calcula la norma infinita de un vector: 		||x|| = max {|xi|}, i = 0, 1, ... n. Metodo que calcula la solucion Ax = b, usando tecnicas iterativas Metodo que calcula la solucion Ax = b, usando tecnicas iterativas Metodo SOR, calcula la solucion de un sistema de ecuaciones 		usando el parametro de relajacion w, que necesariamente varia 		de <0, 2> solo en el caso que la matriz converja. #Garantizar que los elementos de la diagonal principal sean distintos de cero # for j in range(n): # 	dominancia = 0.0 # 	for i in range(n): # 		if j != i: # 			dominancia += fabs(A[i][j]) # 	if A[i][i] < dominancia: # 		exit('La matriz no converge') #		print Xk1 0.  Definiciones 			1.  Gauss 			2.  Gauss-Jordan 			3.  Crout LU1 			4.  Crout L1U 			5.  Doolittle 			6.  LDMt 			7.  Cholesky 			8.  <NAME> (falta) 			9.  Aasen (falta) 			10. Householder 			11. Givens 			12. Jacobi 			13. Gauss-Seidel 			14. SOR 			15. Condicionamiento (incompleto) 			16. Salir #Ver la documentacion #Matrix = [ #		[1, -4, 2, 1], #		[2, -6, 1, 4], #		[-1, 2, 3, 4], #		[0, -1, 1, 1] #		] #El condicionamiento no necesita el vector b #b = [-4, 1, 12, 0] 0.  Sin pivoteo 					1.  Pivoteo parcial 					2.  Pivoteo total #Gauss #Gauss-Jordan #croutLU1, croutL1U, Doolittle #croutLU1 #croutL1U #Doolittle #Imprimir las matrices #calcula las soluciones de Lz=b y Ux=z #imprime las soluciones #LDMt #calcula las soluciones #imprime las soluciones #Cholesky #calcula las soluciones #imprime las soluciones #Parlett y Reid #Aasen #Householder y Givens #Householder #Givens #Comprobar la igualdad #Jacobi #Gauss-Seidel #SOR (Successive Over relaxation) #Condicionamiento 1.  Norma 1 					2.  Norma 2 					3.  Norma infinita 					4.  Condicionamiento con norma 1 					5.  Condicionamiento con norma 2 					6.  Condicionamiento con norma infinita 					7.  Salir #Norma 1 #Norma 2 # print norma2(Matrix) #Norma infinita #Condicionamiento de la matriz usando norma 1 #Condicionamiento de la matriz usando norma 2 #Condicionamiento de la matriz usando norma infinita #Salir del menu #menu()