% main.m
%
%
% May 2013
%
% This script allows to perform some numerical tests by giving the results
% obtained when solving a linear system A*x=b by using the CMRH algorithm
%
% You have to choose :
%    - the size of the system, 
%    - the coefficient matrix, 
%    - the exact solution (so that we can compute the error), 
%    - the tolerance 
%


fprintf('Hello world \n')
fprintf('Enter the matrix choice :\n')
fprintf('1- A(i,j)= (alpha*j-beta) / (n-i+j)  if j <= i    \n')
fprintf('           (alpha*i-beta) / (n-i+j)  if j > i   \n') 
fprintf('2- A(i,j)= i if j <= i    \n')
fprintf('           j if j > i  \n')
fprintf('3- A(i,j)= abs(i-j)+1/(i-j) if j < i  and j > i  \n')
fprintf('           0  if i = j   \n')
fprintf('4- A(i,j)= i+j    \n')
fprintf('5- helsing matrix \n')
fprintf('6- Matrix Market matrix \n')
matrixchoice=input('Enter your choice : ');
if (matrixchoice <6)
    n=input('Enter the size of the matrix : ');
    if (matrixchoice <5)
         A = matrice(n,matrixchoice);
    end
    if (matrixchoice==5)
        A = helsing(n);
    end
else
    filename = input('Enter the matrix file name :','s');
    [A,rows,cols,entries,rep,field,symm] = mmread(filename);
    n = rows;
end

fprintf('Enter the solution choice :\n')
fprintf('1-sol=random vector \n')
fprintf('2-sol=[1,1,...,1]^T \n')
fprintf('3-sol=[1,2,...,n]^T \n')
solchoice = input('Enter your choice : ');
[b,sol] = solution(A,solchoice);

tol = input('Enter the tolerance : ');
iprint = input('Print results : 1-Yes, 2-No : ');


V = A; bc = b;
[x,pseudores,flag] = CMRH(A,b,tol,n,zeros(n,1),iprint);

fid = fopen('pseudores.dat','w');
for i=1:length(pseudores)
    fprintf(fid,'%i \n',pseudores(i));
end
fclose(fid); 

b = x;
if (iprint==1)
    fprintf('\n approx. and exact solutions :\n');
    if (n < 10)
       for i = 1:n
         fprintf(' %6.4e %6.4e\n',b(i),sol(i));
       end
    else
       for i = 1:10
         fprintf(' %6.4e %6.4e\n',b(i),sol(i));
       end
    end
    % we compute the norm of the residu
    bc = -V*b+bc;
    re_cmrh = norm(bc,2);
    % we compute the norm of the error
    sol = -b+sol;
    er_cmrh = norm(sol,2);
    fprintf('residual norm: %6.4e  error norm: %6.4e \n',re_cmrh,er_cmrh);
end