function [] = comparar_metodos_sol(nombrearchivo,tol,maxit)
% funcion para comparar la descomposicion de Cholesky y los metodos del gradiente
% conjugado y el gradiente conjugado precondicionado aplicados a la
% solucion de un sistema de ecuaciones lineales de muchas ecuaciones y
% con matriz sparse.

% Entrada: nombrearchivo puede ser data-1.mat, data-2.mat o data-3.mat.
% tol y maxit se usan para llamar a los metodos iterativos.

% cargando datos desde el archivo data-1, data-2 o data-3.mat
data = load(nombrearchivo);
% Los datos almacenados son: A, b matriz y parte derecha del sistema de
% ecuaciones a resolver para determinar desviacion de membrana
A = data.A;

% A es sparse, el comando para estimar su numero de condicion es condest
disp('condicion de A')
condest(A)

b = data.b;
% Los datos siguientes son necesarios para mostrar la desviacion de la
% membrana, comprenderlos no es parte de este laboratorio
FreeNodes = data.FreeNodes;
Elements = data.Elements3;
Coordinates = data.Coordinates;

figure(1);
spy(A);
title('Estructura de matriz A');

disp('cantidad de elementos distintos de cero en A')
nnz(A)

% resolviendo sistema de ecuaciones
% con descomposicion LU
tic;
R = chol(A);
solchol = R\(R'\b);
tiempochol = toc;
disp('tiempo solucion Cholesky')
tiempochol
disp('cantidad de elementos distintos de cero en R')
nnz(R)
figure(2)
spy(R)
title('Estructura de matriz R');
disp('norma 2 de residuo en solucion LU')
norm(b - A*solchol)/norm(b)
disp('almacenamiento Cholesky')
nnz(A) + nnz(R) + length(b) + length(solchol) 

% resolviendo con metodos iterativos
% Gradiente conjugado
tic;
[solgc,flaggc,relrgc,itgc] = pcg(A,b,tol,maxit);
tiempogc = toc;
disp('tiempo gradiente conjugado')
tiempogc
if flaggc ~= 0
    disp('gradiente conjugado no alcanzo tolerancia requerida')
end
if flaggc == 0
    figure(4)
    mostrarmembrana(solgc,Coordinates,FreeNodes,Elements,'Membrana, GC');
    disp('norma 2 de residuo en solucion gradiente conjugado')
    relrgc
end
disp('almacenamiento GC')
% A, b, rk, rkp1, pk, xk, alphak, betak,
% para calcular betak necesito 2 residuales
% xkp1 puede escribirse encima de xk, pkp1 puede escribirse encima de pk
nnz(A) + length(b) + 4*length(solgc) + 2

% Gradiente conjugado precondicionada
tic;
R = cholinc(A,1e-2);
P = R'*R;
figure(5)
spy(P)
title('Estructura de precondicionador')
[solgcp,flaggcp,relrgcp,itgcp] = pcg(A,b,tol,maxit,P);
tiempogcp = toc;
disp('tiempo gradiente conjugado precondicionado')
tiempogcp
if flaggcp ~= 0
    disp('gradiente conjugado precondicionado no alcanzo tolerancia requerida')
end
if flaggcp == 0
    figure(6)
    mostrarmembrana(solgcp,Coordinates,FreeNodes,Elements,'Membrana, GCP');
    disp('norma 2 de residuo en solucion gradiente conjugado')
    relrgcp
end
disp('almacenamiento GCP')
% A, b, rk, rkp1, pk, xk, zk, zkp1, P, alphak, betak
nnz(A) + length(b) + 6*length(solgcp) + nnz(P) + 2