# -*- coding: utf-8 -*-
"""
Created on Tue Jan 20 13:37:37 2015

Listing_10.3 : Diagonalisation par QR, ˆ l'aide de la fonction 
 QRdec Žcrite pour l'occasion.
@author: grivet

"""
import numpy as np
import scipy.linalg as lg

def affiche_mat(titre,M):
    print (titre)
    for row in M:
        for val in row:
            print (' %10.5f' %val,end = '')
        print()

def QRdec(N,A):
    V = np.zeros((N,N)); Q = np.eye(N,N)
    for k in range(N-1):
        x = A[k:N,k]
        e1 = np.zeros_like(x)
        e1[0] = 1
        v = x + np.sign(x[0])*lg.norm(x)*e1
        v /= lg.norm(v)
        V[k:N,k] = v
        A[k:N,k:N] = A[k:N,k:N] -2*np.outer(v,np.dot(v.transpose(),A[k:N,k:N]))
    for k in range(N):
        Q = Q - 2*np.outer(np.dot(Q,V[:,k]),V[:,k].transpose())
    return Q,A

N = 3
F0 = np.array([[1,1,0.5],[1,1,0.25],[0.5,0.25,2]])
F = np.copy(F0); S = np.eye(N)
nit = 20;
for it in range(nit):
    Q,R = QRdec(N,F)
    F = np.dot(R,Q)
    S = np.dot(S,Q)
affiche_mat('matrice de depart: ',F0)
affiche_mat('matrice diagonalisee: ',F)
affiche_mat('matrice des vecteurs propres: ',S)
affiche_mat('dernier R: ',R)
affiche_mat('dernier Q: ',Q)
print()
print('trace(F0) = %12.6f\t trace(F) = %12.6f' % (np.trace(F0),np.trace(F)))
print('det(F0) = %12.6f\t det(F) = %12.6f ' % (lg.det(F0),lg.det(F)))