Rotina python p/ Inversão de Matrizes pelo método de fatoração QR

Aqui está postado o código em python 2.6.5 que faz a inversão de uma matriz inversível.

# -*- coding: cp1252 -*-

# Inverção da matriz [K] pelo método de fatoração QR

from numpy import *

# Matriz [K]

K = array([[3, 1, 0],[1, 3, 1],[0, 1, 3]],float)

print 'Matriz [K] = '

print K

# Ordem da matriz

n=len(K) # conta o número de colunas da matriz [K]

# Inicia Pg

Pg=zeros([n,n],float)

for i in range(n):

____Pg[i][i]=1

# Fatoração QR (Givens)

for i in range(n-1):

____for h in range(i+1,n):

________if (K[i][i])**2+(K[h][i])**2==0:

____________cs=1

____________sn=0

________else:

____________# Determina cossenos e senos

____________cs=(K[i][i])/((K[i][i])**2+(K[h][i])**2)**0.5

____________sn=(K[h][i])/((K[i][i])**2+(K[h][i])**2)**0.5

________for j in range(n):

____________k1,k2 = 0, 0

____________k1,k2 = K[i][j], K[h][j]

____________# Cálculo da matriz ~K ( vai transformando K em Sg)

____________K[i][j] = k1*cs + k2*sn

____________K[h][j] = -k1*sn + k2*cs

____________# Cálcula Pg transposto

____________p1,p2 = 0, 0

____________p1,p2 = Pg[i][j], Pg[h][j]

____________Pg[i][j] = p1*cs + p2*sn

____________Pg[h][j] = -p1*sn + p2*cs

# Obtém a matriz inversa de [K] armazenando cada vetor coluna ci resultante da solução

#do sistema Sg.ci=li onde li é o vetor linha de Pg.

# No código abaixo K foi tranformado em Sg.

R=zeros([n,1],float)

for k in range(n):

____for i in range(n):

________R[i][0]=Pg[i][k]

____U=zeros([n,1],float)

____for i in range(n):

________SS=0

________for j in range(n):

____________SS=SS+(K[n-1-i][n-1-j]*U[n-1-j])

________U[n-1-i]=(R[n-1-i]-SS)/K[n-1-i][n-1-i]

____for i in range(n):

________Pg[i][k]=U[i][0]# Matriz inversa [K]-1

print Pg

# Obs: Ao copiar o arquivo, antes de rodar a rotina é necessário editar o código substituindo cada underline por um espaço.

Para que esta rotina funcione é necessário que vc tenha instalado em seu computador o compilador python e o modulo de ferramentas numpy. Ambos podem ser baixados do site

http://www.python.org.

Posteriormente, vou editar esta publicação adicionando a teoria que fundamenta o código desta rotina.