function [A,f,u,n]=LoadMatrixSymphony(nx)
%
% Function LoadMatrixSymphony
%
% [A,f,u,n]=LoadMatrixSymphony(nx);
%
% Load a non-symmetric benchmark matrix detailed [1] or used in 
% the references below.
%
% Input: nx, the level of discretization of the problem (Number
%           of elements in the x-direction).  There are three
%           supported values: nx=16, 32, or 64.  
%        files. The data files (matrix16.dat, matrix32.dat, 
%           and matrix64.dat) should be available.
%
% Output: A, the sparse benchmark matrix. Note that A is banded 
%            with only 9 non-zero diagonals.
%         n, the dimension of A.
%         f, the right-hand side vector.
%         u, the solution vector, A * u = f.
%
% Christophe Lecomte
% Boston University
% Department of Aerospace and Mechanical Engineering
%
% 11th of June 2004 (Day of the funeral of former US president Ronald Reagan)
% 
% References:
% 
% [1] C. Vuik. Non-symmetric benchmark problem. Available on the Internet at
%     http://ta.twi.tudelft.nl/nw/users/vuik/laplace/problem.pdf.
% [2] E.F.F. Botta, K. Dekker, Y. Notay, A. van der Ploeg, C. Vuik, F.W. Wubs, 
%     and P.M. de Zeeuw. How fast the Laplace equation was solved in 1995. 
%     Appl. Num. Meth., 24:439R-455, 1997. 
% [3] H.A. van der Vorst and C. Vuik. GMRESR: a family of nested GMRES methods. 
%     Num. Lin. Alg. Appl., 1:369-386, 1994. 
% [4] C. Vuik. Solution of the discretized incompressible Navier-Stokes equations 
%     with the GMRES method. Int. J. for Num. Meth. Fluids, 16:507-523, 1993. 
% [5] C. Vuik. Fast iterative solvers for the discretized incompressible 
%     Navier-Stokes equations. Int. J. for Num. Meth. Fluids, 22:195-210, 1996.
%

% Initialization
% --------------
% Check input 
%
if (nx ~= 16 & nx ~= 32 & nx ~= 64)
    error('The supported values for nx are 16, 32, or 64 - LoadMatrixSymphony.')
end
%
ny=4*nx;
n=(nx+1)*(ny+1);
%
% Read File
% ---------
% Open file
%
file=strcat('matrix',num2str(nx),'.dat');
fid=fopen(file,'rt');
if (fid == -1)
    error(strcat('The file ',file,' could not be opened - LoadMatrixSymphony.'))
end
%
% Read matrix
%
A=zeros(n,9);
for j=1:9
    for i=1:n
        A(i,j)=str2num(fgetl(fid));
    end
end
%
% Read solution vector
%
for i=1:n
    u(i)=str2num(fgetl(fid));
end
u=u(:);
%
% Read right-hand side vector
%
for i=1:n
    f(i)=str2num(fgetl(fid));
end
f=f(:);
%
% Close file
%
fclose(fid);
if (fid == -1)
    error(strcat('The file ',file,' could not be closed - LoadMatrixSymphony.'))
end
%
% Put matrix in sparse format
% ---------------------------
% Shift zeros to bottom or top of array 
% The first row of spdiags' input matrix contains the elements of
% the first column of the square matrix while the matrix read from the
% file contains the first elements of the first row of the square matrix. 
n1=nx+1;
A(:,2)=[A(1+n1+1:end,2);zeros(n1+1,1)];
A(:,3)=[A(1+n1:end,3);zeros(n1,1)];
A(:,4)=[A(1+n1-1:end,4);zeros(n1-1,1)];
A(:,5)=[A(1+1:end,5);zeros(1,1)];
A(:,6)=[zeros(1,1);A(1:end-1,6)];
A(:,7)=[zeros(n1-1,1);A(1:end-n1+1,7)];
A(:,8)=[zeros(n1,1);A(1:end-n1,8)];
A(:,9)=[zeros(n1+1,1);A(1:end-n1-1,9)];
%
A=spdiags(A,[0,-n1-1,-n1,-n1+1,-1,1,n1-1,n1,n1+1],n,n);
%
