% rutero preg3, lab4, 521230, S2-2013
elsurdata = load('elsurdata');

G = elsurdata.Gelsur;
Urls = elsurdata.Uelsur;

A = creatematrix_pagerank(G);

[L,U] = lu(A);
spy(A)
title('spy(A)')

figure()
spy(L)
title('spy(L)')

figure()
spy(U)
title('spy(U)')

disp('numero de elementos distintos de cero en A,L,U');
nnz(A)
nnz(L)
nnz(U)

% si A no es simetrica, no es posible aplicar gradiente conjugado

% descomponer A para crear matrices de iteracion de Jacobi y Gauss-Seidel
% matriz diagonal con diagonal principal igual a diagonal principal de A
D = diag(diag(A));
% parte triangular inferior de A
E = -tril(A,-1);
% parte triangular superior de A
F = -triu(A,1);

% matriz de iteracion de Jacobi
J = D\(E+F);
% matriz de iteracion de Gauss-Seidel
GS = (D-E)\F;
% radios espectrales de ambas matrices
disp('radios espectrales de matrices de iteracion, Jacobi y GS')
max(abs(eigs(J)))
max(abs(eigs(GS)))

norm(J,inf)
norm(GS,inf)
norm(J,1)
norm(GS,1)

% intentanto con Jacobi
[jacsol,convjac,itsjac] = jacobisol(A,ones(size(A,1),1),...
    zeros(size(A,1),1),1e-6,100);
[gssol,convgs,itsgs] = gaussseidelsol(A,ones(size(A,1),1),...
    zeros(size(A,1),1),1e-6,100);
if convjac && convgs
    disp('norma de diferencia entre aproximación de Jacobi y aproximación de Gauss-Seidel')
    norm(jacsol-gssol,inf)
    disp('número de iteraciones de Jacobi')
    itsjac
    disp('número de iteraciones de Gauss-Seidel')
    itsgs
    mostrarrangos(gssol,G,Urls)
elseif ~convjac || ~convgs
    if ~convjac
        disp('el método de Jacobi no encontró aproximación con precisión deseada')
    else
        disp('el metodo de Gauss-Seidel no encontró aproximación con precisión deseada')
    end
end
