from numpy import *

def  householder(A): # computes a QR factorization upper triangular matrix R 
    """matrices W, R pour la factorisation QR'
    
    Usage: 
    W, R = householder(A)

    Argument:
    A: matrice mxn avec m>=n
    
    Retour:
    W: matrice de vecteurs de reflexions de householder
    R: matrice R de la factorisation QR

    Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia 
    """
    m,n = A.shape
    ident = eye(m)
    W = zeros((m,n))
    R = copy(A)
 
    for k in range(0,n):
        e1=ident[k:m,k]
        x = R[k:m,k]
        W[k:m,k]= ( int(x[0]>=0) - int(x[0]<0) )*linalg.norm(x,2)*e1+x
        W[k:m,k]=W[k:m,k]/linalg.norm(W[k:m,k],2)
        R[k:m][:,k:n]=R[k:m][:,k:n]-dot(outer(2*W[k:m,k],W[k:m,k]),R[k:m][:,k:n])
 
    return W, R
 
def formQ(W): 
    """Matrice Q de la factorisation QR
    
    Usage: 
    Q = formQ(W)
    
    Argument:
    W: matrice mxn de reflexions obtenue de la décomposition Householder
    
    Retour:
    Q: Matrice Q de la factorisation QR
    
    Note:
    La matrice W se calcule avec W, R = householder(A)
    """
    
    m,n = W.shape
    Q=zeros((m,m))
 
    # compute the product QI = Q, column by column
    for i in range(0,m):
        Q[i,i]=1
        for k in range(n-1,-1,-1): 
            Q[k:m,i]=Q[k:m,i]-2*dot(outer(W[k:m,k],W[k:m,k]),Q[k:m,i])
 
    return Q

n = 5
A0 = random.random((n,n))-0.5
A0 = dot(A0.T,A0) # on construit une matrice réelle symétrique
print A0

A = A0
# Algorithme QR 'pur'
for k in range(0,20):
    W, R = householder(A)
    Q = formQ(W)
    A = dot(R,Q)

# affiche A avec 4 décimales
for i in range(0,n):
    for j in range(0,n):
        print '{:+.4f}'.format(A[i,j]),
    print '\n'