# -*- coding: utf-8 -*-
"""
Created on Thu Dec 18 17:37:35 2014

@author: grivet
Listing 6.1 : Factorisation LU sans permutation de lignes
"""
import numpy as np

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

N = 4
A = np.random.rand(N,N)+ np.eye(N,N)
print( 'A = ')
affiche_mat(A)
print()
A0 = np.copy(A)
#pas A0 = A, qui copie dans AO l'adresse de A:
#on aurait alors deux noms pour la même matrice

for k in range(N-1):
#    lignes = range(k+1,N)
    A[k+1:,k] = - A[k+1:,k]/A[k,k]
    A[k+1:,k+1:] = A[k+1:,k+1:] + np.outer(A[k+1:,k],A[k,k+1:])
#Comme son nom l'indique, outer(a,b) forme le produit extérieur
# des deux vecteurs: colonnne X ligne = matrice
#    print 'A = '
#    affiche_mat(A)
U = np.zeros_like(A); L = np.eye(N,N)
#L est la matrice identité d'ordre N
for k in range(N):
   U[k,k:] = A[k,k:]
for k in range(1,N):
   L[k,0:k] = - A[k,0:k]
print ('L = ')
affiche_mat(L)
print ('U = ')
affiche_mat(U)
print ('L*U = ')
affiche_mat(np.dot(L,U)); print()
print("A0 = ")
affiche_mat(A0)