% ============
% Apartado 1.1
% ============
%
%
% ##################### Inicio del fichero ####################
function R=opera(A,operacion,B)
%
%function function R=opera(A,operacion,B)
%Realiza operaciones con una y otra matriz dependiendo de la operacion
%Ejemplo de uso de la instrucción case
%
%Entrada
%	A,B         matrices
%	operacion   suma('+'),resta('-'),multiplicación('*') o inversión(#')
%Salida
%	R           resultado
%
% Verificacion de los argumentos de entrada
switch nargin
    case 2 % La unica operacion admitida es la inversion
        if operacion~='#';error('Dos argumentos y la operacion no es inversion');end
        [fA,cA]=size(A);
        if fA~=cA;error('La matriz a invertir no es cuadrada');end
        if fA~=rank(A);error('La matriz a invertir no es regular');end
        R=inv(A);
    case 3
        [fA,cA]=size(A);[fB,cB]=size(B);
        switch operacion
            case {'+','-'}
                if fA~=fB|cA~=cB;error('Las matrices no son congruentes');end
                if operacion='-';signo=-1;else;signo=1;end
                R=A+signo*B;
            otherwise
                error('Operación no valida');
        end
    otherwise
        error('Número de argumentos no válido');
end
% ##################### fin del fichero ####################
%
% * * * EJERCICIO 1 * * *
% Añadir las siguientes instrucciones case
            case '*'
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A*B;
            case '\';
                if fA~=cA;error('La matriz a invertir no es cuadrada');end
                if fA~=rank(A);error('La matriz a invertir no es regular');end
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A\B;
            case '/';
                if fB~=cB;error('La matriz a invertir no es cuadrada');end
                if fB~=rank(B);error('La matriz a invertir no es regular');end
                if cA~=fB;error('Las matrices no son congruentes');end
                R=A\B;
% ============
% Apartado 1.2
% ============
A=[2,2,-1,-2;4,0,2,1;1,3,-2,0;3,2,0,5];b=[10;7;3;2];x=A\b
inv(A)\b
max(abs(eig(A))) % Radio espectral de la matriz A
%
% * * * EJERCICIO 1 * * *
A=[112,84,48,28;84,81,77,-7;48,77,146,-83;28,-7,-83,85];
b=[-204;-186;-193;38];
x=A\b
% Fin del ejercicio
orden=round(logspace(1,5,40));N=length(orden);t=zeros(N,1);
for n = 1:N
	A = rand(n,n);b = rand(n,1);
	tic; x = A\b;t(n) = toc;
end
plot(t)
%
% * * * EJERCICIO 3 * * *
orden=round(logspace(1,5,40));N=length(orden);t=zeros(N,1);
for n = 1:N
	A = tril(rand(n,n));b = rand(n,1);
	tic; x = A\b;t(n) = toc;
end
plot(t)
% Fin del ejercicio
% ============
% Apartado 1.3
% ============
Ab=[A,b] % Matriz ampliada
i=1;j=2;m=Ab(j,i)/Ab(i,i); % calculo del pivote
Ab(j,:)=Ab(j,:)-m*Ab(i,:); % anulamos la fila j-esima
% * * * EJERCICIO 4 * * *
A=[112,84,48,28;84,81,77,-7;48,77,146,-83;28,-7,-83,85];
b=[-204;-186;-193;38];
[xg,mg,AB,Lg]= T31gaussps(A,b,1)
% * * * EJERCICIO 5 * * *
% Sustituir el cálculo del pivote por las siguientes líneas
        escala  =transpose(max(abs( A(i:n,i:n)'))); %Columna con el maximo de cada fila
        if ~all(escala); error('sistema indeterminado o incompatible');end
        pivot=abs(A(i:n,i)./escala);
        [m,ifil]=max(pivot); %Maximo de la columna escalada y posición

% * * * EJERCICIO 6 * * *
% Sustituir el cálculo del pivote y el intercambio de filas por las siguientes
% líneas.
        [mval,ifil]=max(abs( A(i:n,i:n)) );  %Maximo de cada columna
        [mval,icol]=max(mval);               %Maximo global
        nfil=ifil(icol)+i-1;ncol=icol+i-1;   %Posicion del maximo
        if mval==0; error('sistema indeterminado o incompatible');end
        if icol>1;
            j=icol+i-1;
            p=A(:,i);A(:,i)=A(:,j);A(:,j)=p;  %Intercambio de columnas
            p=M(i,2);M(i,2)=M(j,2);M(j,2)=p;%Historico del intercambio columnas
        end
        if ifil(icol)>1; 
            j=ifil(icol)+i-1;
            p=A(i,:);A(i,:)=A(j,:);A(j,:)=p; %Intercambio de filas
            p=M(i,1);M(i,1)=M(j,1);M(j,1)=p; %Historico del intercambio de filas 
        end

% ============
% Apartado 1.4
% ============
% * * * EJERCICIO 7 * * *
A=[112,84,48,28;84,81,77,-7;48,77,146,-83;28,-7,-83,85];
b=[-204;-186;-193;38]
[L,U,M] = lu(A)
y=L\(M*b)
x=U\y
L= chol(A,'lower')
y=L\b
x=L'\y
% ============
% Apartado 1.5
% ============
orden=10;A=hilb(orden);b=A(:,orden);x=A\b
c1=cond(A,1),c2=cond(A,2),ci=cond(A,inf)
%
c=round(b*100000)/100000;y=A\c
%
e1=100*norm(x-y,1)/norm(x,1) % error porcentual
co1=100*c1*norm(b-c,1)/norm(b,1) % cota del error