function [solgs,tolalcanzada,iteraciones] = gaussseidelsol(A,bv,varargin)
% se resuelve sistema Ax = b con metodo de Gauss-Seidel
% los primeros dos parametros de entrada son matriz A y parte derecha b
% se resolveran tantos sistemas como columnas tenga b
% los sistemas a resolver son Ax=b(:,i)
% cada una de las soluciones se almacena en solgs
% los restantes parametros de entrada son opcionales, ellos son
% tolerancia, aproximacion inicial y maximo numero de iteraciones. 

% Si se dan, deben darse en
% cell-arrays (uno por cada parametro opcional de entrada)
% si se quiere dar un valor a x0, un cell array debe ser
% tal que su primera componente sea la cadena de caracteres 'x0' y su
% segunda, el vector que debe usarse como aproximacion inicial a la
% solucion exacta del sistema Ax = b. Valor por defecto: vector nulo.
% Para, por ejemplo, resolver Ax=b y usar un vector de 1 como aproximacion
% inicial, se debe escribir
%[sol, tol, it] = gaussseidelsol(A,b,{'x0',ones(length(b),1)})

% si se quiere dar valor a la tolerancia con que se desea calcular la
% aproximacion a la solucion exacta de Ax=b, un cell array de
% entrada debe ser de la forma {'tol', valor}. Valor por defecto: 1e-6

% si se desea dar un valor al maximo numero de iteraciones, un cell
% array de entrada debe ser de la forma {'maxiter',valor}. Valor por
% defecto: 100.

% Por ejemplo, si se quiere resolver un sistema Ax = b usando [1,1] 
% como aproximacion inicial a su solucion exacta, 
% con una tolerancia de 1e-4 y permitiendo un maximo de 10 iteraciones, 
% el llamado a esta funcion debe ser
% [gssol,co,it] = gaussseidelsol(A,b,{'tol',1e-4},{'maxiter',10},{'x0',[1 1]})

% Chequeando dimensiones correctas en A y b
[mb,nb] = size(bv);
[mA,nA] = size(A);
% A debe ser cuadrada y el numero de filas de b debe ser igual al tamanno
% de A
if mA ~= nA || mA ~= mb
    error('entrada erronea')
end

% se resuelven tantos sistemas de ecuaciones como columnas tiene b

% valores por defecto para x0, maxiter y tol
% se resolveran nb sistemas de ecuaciones con el metodo de Gauss-Seidel
% usando como tolerancia 1e-6, como vector inicial el vector nulo
% y a lo sumo 100 iteraciones del metodo
x0v = zeros(mb,nb);
tolv = 1e-6*ones(nb,1);
maxiterv = 100*ones(nb,1);

% comprobando parametros opcionales
if ~isempty(varargin)
    for i = 1 : length(varargin)
        input = varargin{i};
        if strcmp(input{1},'x0')
            % x0 es parametro de entrada
            x0values = input{2};
            [m,n] = size(x0values);
            if m == mA && n == nb
                x0v = x0values;
            elseif m == mA && n == 1
                x0v = repmat(x0values,1,nb);
            end
        elseif strcmp(input{1},'maxiter')
            % maxiter es parametro de entrada
            maxitervalues = input{2};
            maxitervalues = maxitervalues(:);
            n = length(maxitervalues);
            if n == 1
                maxiterv = maxitervalues(1)*ones(nb,1);
            elseif n == nb
                maxiterv = maxitervalues;
            else
                error('maxiter: entrada erronea');
            end
        elseif strcmp(input{1},'tol')
            % tolerancia es parametro de entrada
            tolvalues = input{2};
            tolvalues = tolvalues(:);
            n = length(tolvalues);
            if n == 1
                tolv = tolvalues*ones(nb,1);
            elseif n == nb
                tolv = tolvalues;
            else
                error('tol: entrada erronea');
            end
        end           
    end
end

% en sol se guarda la aproximacion a la solucion exacta de los sistemas
% a resolver
solgs = zeros(mb,nb);
% por cada uno de los sistemas, tolalcanzada(i) indica si se alcanzo la
% precision deseada en la solucion numerica
tolalcanzada = zeros(nb,1);
% iteraciones(i) indicara el numero de iteraciones necesarias para obtener la
% precision deseada (si tolalcanzada(i) es 1) o sera igual al numero maximo
% de iteraciones (si tolalcanzada(i) es 0)
iteraciones = zeros(nb,1);

% descomposicion de A
D = diag(diag(A));
E = -(tril(A)-D);
F = -(triu(A)-D);
% matriz de iteracion del metodo de Gauss-Seidel, no es necesario calcularla para
% implementar el metodo, pero queremos comprobar si su radio espectral es
% estrictamente menor que 1
MGS = inv(D-E)*F;
% valores propios de MGS
abs(eig(MGS))
for i = 1 : nb
    % numero maximo de iteraciones para resolver sistema i-esimo
    maxiter = maxiterv(i);
    % tolerancia con que se quiere obtener aproximacion a solucion exacta
    % de sistema i-esimo
    tol = tolv(i);
    % aproximacion inicial para solucion de sistema i-esimo
    x0 = x0v(:,i);
    % necesario para implementar criterio de parada en laminas del curso
    x00 = zeros(length(x0),1);
    % parte derecha de sistema i-esimo
    b = bv(:,i);
    % inicializando numero de iteracion actual
    it = 1;
    % inicializando variable que chequea si se alcanzo la tolerancia 
    % requerida en la solucion aproximada
    conv = tolalcanzada(i);    
    while ~conv && it <= maxiter
        % inicializar vector donde se almacena iterando a calcular en paso
        % actual
        xnew = zeros(mb,1);
        % ciclo del metodo de Gauss-Seidel
        for j = 1 : nA
            suma = 0;
            for k = 1 : j-1
                suma = suma + A(j,k)*xnew(k);
            end
            for k = j+1 : nA
                suma = suma + A(j,k)*x0(k);
            end
            xnew(j) = (b(j)-suma)/A(j,j);
        end
        % se ha calculado el iterando actual, se comprueba si se alcanzo la
        % tolerancia requerida, para ello debe comprobarse si
        % coef*norm(diferencia entre iterando actual e iterando anterior)
        % es menor o igual que la tolerancia requerida, coef se especifica
        % en las transparencias del curso
        if it >= 2
            nM = norm(xnew-x0)/norm(x0-x00);
            coef = nM/(1-nM);
            if coef <= 0
                coef = 1;
            end
        else
            coef = 1;
        end
        conv = (coef*norm(xnew-x0) <= tol);          
        x00 = x0;
        x0 = xnew;
        % aumentando numero de iteraciones
        it = it + 1;
    end
    solgs(:,i) = xnew;
    tolalcanzada(i) = conv;
    iteraciones(i) = it-1;
end