# -*- coding: utf-8 -*-
"""
Created on Tue Dec 23 16:17:33 2014

@author: grivet
Listing 6.2 : Factorisation LU, avec permutation de lignes

"""
from pylab import *

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

N = 4;
A = array([ [2, 1, 0, 4],[-4, -2, 3, -7],[4, 1, -2, 8],[0,-3,-12, -1]],double)
affiche_mat('A = ',A)
print()
A0 = copy(A)
#p = zeros(N)
for k in range (N-1):
    m = amax(abs(A[k:,k]))
    km = argmax(abs(A[k:,k]))
    km = km + k 
    temp = copy(A[k,:])
    A[k,:] = A[km,:]; A[km,:] = temp
    p[k] = km;
#    print m,km,k
#    affiche_mat('apres permut: ',A)
#    print
#   pause;
    if A[k,k] != 0:
        A[k+1:,k] = - A[k+1:,k]/A[k,k]
        A[k+1:,k+1:] = A[k+1:,k+1:] + outer(A[k+1:,k],A[k,k+1:])
#        affiche_mat('apres elim: ', A)
#        print
#   pause;
affiche_mat('A = ',A)
U = zeros((N,N))
L = eye(N,N)
PLU = zeros((N,N))
for k in range(N):
    U[k,k:] = A[k,k:]
    L[k,:k] = - A[k,:k];
affiche_mat('U = ',U)
affiche_mat('L = ',L)
LU = dot(L,U)
affiche_mat('L*U = ',LU)
for k in range(N):
    PLU[k,:] = LU[p[k],:]
affiche_mat('matrice de depart = ', A0)
affiche_mat('matrice L*U permutee = ',PLU)