% pregunta 2, laboratorio 4, 521230, S2-2012

N = 50;

% crear matriz tridiagonal directamente como sparse
A = matriz_preg1_lab4_sparse(N,4,-1,-1,1);
% parte derecha de sistema a resolver
b = ones(N,1);
for tol = [1e-2, 1e-4, 1e-6]
    % gradiente conjugado, gradiente conjugado precondicionado
    [xgc,fgc,relres,itscg] = pcg(A,b,tol,100);
    % si fgc es distinto de cero, no se pudo calcular aproximacion
    if ~fgc
        disp('número de iteraciones de gradiente conjugado')
        itscg
        disp('residuo relativo')
        relres
    end
    % gradiente conjugado precondicionado, buscando precondicionador con ichol
    % (cholinc es obsoleta en última versión de matlab)
    tolp = 0.2;
    % L = ichol(A,struct('type','ict','droptol',tolp));
    % P = L*L';
    R = cholinc(A,tolp);
    % llamado a pcg con valor para P (gradiente conjugado precondicionado)
    [xgcp,fgcp,relres,itspcg] = pcg(A,b,tol,100,R'*R);
    if ~fgcp
        disp('número de iteraciones de gradiente conjugado precondicionado')
        itspcg
        disp('residuo relativo')
        relres
    end
end

% gradiente conjugado con matriz mal condicionada
A = gallery('poisson',10);
condest(A)
b = ones(100,1);
xlu = A\b;

[xgc,fgc,relres,itscg] = pcg(A,b,1e-6);
% si fgc es distinto de cero, no se pudo calcular aproximacion
if ~fgc
    disp('número de iteraciones de gradiente conjugado')
    itscg
end
for tolp = [0.2 0.1 0.01 0.001]
    % gradiente conjugado precondicionado, buscando precondicionador con ichol
    % (cholinc es obsoleta en última versión de matlab)
    % L = ichol(A,struct('type','ict','droptol',tolp));
    % P = L*L';
    R = cholinc(A,tolp);
    norm(R'*R-A,inf)
    % llamado a pcg con valor para P (gradiente conjugado precondicionado)
    [xgcp,fgcp,relres,itspcg] = pcg(A,b,1e-6,100,R'*R);
    if ~fgcp
        disp('número de iteraciones de gradiente conjugado precondicionado')
        itspcg
    end
end