%
% This is a test driver for the
% testing the LU factorization
% routines.
%
% The actions of this routine are controlled
% by the parameters set at the beginning.
%
% Tests are run for a range of values of the dimension n
% (the outer loop of code)
%
% For each value of n matrices problems are generated
% in a variety of ways.
% (the main inner loop of code)
% This can be done based on a custom 
% specification of A and x and/or b
% (which the user must place in the code)
% or a random number generator based matrix
% with some structural contraints imposed.
% Depending on the problem specified
% Matlab's library and commands may be used to determine
% the "true" answer for comparison to the
% computed answer.
% The matrices and vecotor required are then
% placed in the data structures as required
% by the assignment, the test is run,
% the error information collected.
%
% During this loop, if desired, detailed
% information on the results of each test problem
% can be output. Or, a summary (mean and max) of the
% relative error can be output for each set of tests
% with the same n. A third output option collects the
% information and plots the means and maximums of each set of
% tests with the same n vs the value of n.
% In all cases at the end the maximum value of the relative
% error over all tests is printed.
%

%
% The array M  is not n by n in order
% to have a code similar to that used
% when a range of problems would use the same
% data structure or when other information
% would be stored in M in addition
% to L and  U on output and A on input.
% The same is true for the  vectors b, xcomp, and xtrue.
%
%
% These set the dimensions of the  data structures.
% All problems tested should have n smaller than the smallest of these.
Mrdim=100;
Mcdim=120;
xdim=150;
bdim=110;

%
% set range of matrix dimensions and number
% of tests per dimension
%
% make sure this is made consistent with the dimension
% of the specially structured problems in problemtype == 1
% when they are a single size.
%

nmin=6;
nmax=6;
ndelta=1;
ntests=1;

%
% Indicate the subroutine to be tested by commenting out those not desired.
%
%routinenum = 1; % LU without pivoting
routinenum = 2; % LU with partial pivoting of rows, PA=LU
%routinenum = 3; % LU with complete pivoting PAQ = LU

%
% If you want to have a known solution then
% generate the corresponding right-hand side
% vector here by:
%     1. using b as the known solution xtrue
%        and generating a new b. 
%        (set rhstype = 1)
%     2. using some specified known solution
%        and generating a new b. 
%        (set rhstype = 2 and set the known solution xtrue)
%
%     3. To use b as the right-hand side vector
%        set rhstype to anything other than 1 and 2.
%
    rhstype = 1; 

%
% Set the manner of matrix and vector generation by commenting out those not desired.
%
% In all cases, when generating A via factors, L and U have
% the magnitudes of their elements rescaled to reduce the chance of ill-conditioning.
% Note that if A is produced by generating L and U then setting A=LU. The computed
% factors may be different due to use of pivoting.
% The factors produced by a library routine and the tested routine
% may not be the same due to different pivoting choices (there may be more than one
% element with a maximal magnitude).
%
%___________________________________________________________________________________
problemtype = 1; % User custom manner. Add the code at the appropirate place below.

%___________________________________________________________________________________
%problemtype = 2; % L and U by uniform random numbers on the interval [dmin,dmax].
                 % Note that this interval's endpoints can be set to include both signs
                 % or have a fixed one based on the signs of the endpoints.
                 % this does not affect the signs of the 1's on the diagonal of L.
                 % A is calculated from the product LU
%dmin=-10.0;
%dmax=10.0;
%umin=1.0;
%umax=10.0;
%___________________________________________________________________________________

%problemtype = 3; % L and U by normal random numbers with specified mean and standard deviation.
                 % A is calculated from the product LU
%tempmean=0.0;
%tempstdev=500.0;

%___________________________________________________________________________________

%problemtype = 4; % A filled with uniform random numbers on the interval [dmin,dmax].
                 % Note that this interval's endpoints can be set to include both signs
                 % or have a fixed one based on the signs of the endpoints.
%dmin=-10.0;
%dmax=10.0;
%___________________________________________________________________________________

%problemtype = 5; % A filled with normal random numbers with specified mean and standard deviation.
%tempmean=0.0;
%tempstdev=500.0;



%___________________________________________________________________________________
%
% Set output type by commenting out those not desired.
%

%outputtype = 1; % norm information listed for every test problem
%outputtype = 2; % only mean/max norm information listed for each set with fixed n
%outputtype = 3; % mean/max norm vs n plotted (the arrays are available for 
                % examination after driver runs in maxforthisn and meanforthisn arrays
outputtype = 4  % special display for a specific test problem
                % typically from problemtype = 1 or from a single test problem
                % to verify expected structure of factors or other properties


%___________________________________________________________________________________
%
% There are no parameters for the user to set below this point
% except the problemtype = 1 code to set the the manner in which
% the problem is generated, i.e., the set of numbers from which
% elements of A, xtrue and b are generated.
%


%___________________________________________________________________________________
%
% initialize max relative error and max condition number for all  tests  of all sizes
%
maxrelerror = 0.0;
Amaxrelerror = 0.0;
Akappa = 1.0;

%___________________________________________________________________________________
%
% Output Label of Tests to the Matlab Terminal
%

if routinenum == 1
  fprintf(1,'\n**************************************** \n')
  fprintf(1,'LU without pivoting algorithm tests for sizes n = %g to %g by %g \n',nmin,nmax,ndelta)
  fprintf(1,'**************************************** \n')
elseif routinenum == 2 
  fprintf(1,'\n**************************************** \n')
  fprintf(1,'LU with partial (row) pivoting algorithm tests for sizes n = %g to %g by %g \n',nmin,nmax,ndelta)
  fprintf(1,'**************************************** \n')
elseif routinenum == 3
  fprintf(1,'\n**************************************** \n')
  fprintf(1,'LU with complete pivoting algorithm tests for sizes n = %g to %g by %g \n',nmin,nmax,ndelta)
  fprintf(1,'**************************************** \n')
else
  fprintf(1,'\n**************************************** \n')
  fprintf(1,'Undefined version number   = %g \n',routinenum)
  fprintf(1,'**************************************** \n')
  return; % terminate the tester and return to command prompt
end
%___________________________________________________________________________________
%
% Initialize error saving vectors
% 
%
AError2norm= zeros(ntests,1);
ARelError =  zeros(ntests,1);
Atrue2Norm = zeros(ntests,1);
Error2Norm = zeros(ntests,1);
RelError =  zeros(ntests,1);
xtrue2Norm = zeros(ntests,1);

% Initialize arrays to collect mean/max norm for every set with fixed n
nexperiments=size([nmin:ndelta:nmax]',1);
Ameanforthisn = zeros(nexperiments,1);
Amaxforthisn =  zeros(nexperiments,1);
meanforthisn = zeros(nexperiments,1);
maxforthisn =  zeros(nexperiments,1);
experimentindex = 0;

%___________________________________________________________________________________
%Test loop over a range of dimensions of L, U, and v

for n = nmin:ndelta:nmax

  experimentindex = experimentindex + 1;
  if (outputtype == 1) || (outputtype == 2)
    fprintf(1,'\n********************************** \n')
    fprintf(1,'Tests for size n = %g \n',n)
    fprintf(1,'********************************** \n')
  end

% for fixed n  perform ntests random problems
  for itest=1:ntests

% reinitialize data structures for a new problem with size n
    M=zeros(Mrdim,Mcdim);
    xcomp=zeros(xdim,1);
    xtrue=zeros(xdim,1);
    b=zeros(bdim,1);
% and the array to store the product
% of the computed L and U returned
% from the LU factorization
    Acomp=zeros(Mrdim,Mcdim);


%
% Generate A and/or fully expanded L and U
% along with a vector b
% and place them each of them in their
% associated work array.
%
% Note that L and U are rescaled to mitigate large condition numbers
%

% Save current state of random number generator.
% This can be used for reproducing the same sequence of
% test problems if needed. That is not done in this
% form of the code but it is often useful for debugging.
% the rng state can be restored with the statement rng(statesave);
% See the Matlab documentation for more details.
  statesave=rng;

% Select the type of random numbers  or other special case
% setting of the problems parameters.
% Also set the title of the error plot for outputtype = 3 and
% outputtype 4  by commenting/uncommenting the appropriate lines
%
  if problemtype == 1
%  
%    outputtitle = "Strictly diagonally dominant by row and column";
%    dmin=-1.0;
%    dmax=1.0;
%    Temp = dmin.*ones(n,n) + (dmax-dmin).*(rand(n,n)); 
%    Temp = Temp + 5.0*n*diag(ones(n,1));

%    outputtitle = "Lower triangular with diagonal element large so no partial pivoting should be done";
%   uses Temp from strictly diagonally dominant by row and column above or specifically defines it from a toeplitz form.
%   Note this means that the multipliers are computed and are not equal to the original elements
%   since L is unit lower triangular. U is then diagonal with elements equal to the diagonal of the
%   original lower triangular matrix.  Restricting the original matrix to a unit lower triangular
%   gives the simplest case of L = A  U = I


%    Temp = tril(Temp,0);

%    c=[10; 3; 1; 2; 5; 4;];
%    r=[10 0 0 0 0 0];
%    disp('Temp data structure initialize')
%    Temp = toeplitz(c,r)
%%    disp('condition number of tridiagonal matrix')
%    cond(Temp(1:n,1:n))

%
% Tridiagonal with large first subdiagonal  so partial pivoting should be done
% this is constructed to be Toeplitz tridiagonal for convenience only, i.e.,
% Toeplitz affects the answer but not the algorithm.

    outputtitle = "Tridiagonal with large first subdiagonal  so partial pivoting should be done";    
    c=[1; 3; 0;0;0;0;0;0;0;0];
    r=[1 2 0 0 0 0 0 0 0 0];
    disp('Temp data structure initialize')
    Temp = toeplitz(c,r)
    disp('condition number of tridiagonal matrix')
    cond(Temp(1:n,1:n))


%    outputtitle = "Lower triangular with off diagonal elements large so  partial pivoting is required";
%   


% this is constructed to be Toeplitz tridiagonal for convenience only, i.e.,
% Toeplitz affects the answer but not the algorithm.

%    c=[3; 1; 2; 5; 6; 10];
%    r=[3 0 0 0 0 0];
%    disp('Temp data structure initialize')
%    Temp = toeplitz(c,r)
%    disp('condition number of tridiagonal matrix')
%    cond(Temp(1:n,1:n))




%
% Example of exponential growth factor. Take care setting nmin, nmax and ndelta for
% this matrix family. It is suggested to only use one small example like the others
% in this problem type and use outputtype = 4
%
%Temp = eye(n,n) - tril(ones(n,n),-1);
%Ident = eye(n,n);
%Temp = Temp + ones(n,1)*Ident(n,1:n);
%Temp(n,n)=1.0;

%disp('Temp data structure initialize')
%Temp
%disp('condition number of scaled matrix')
%cond(Temp(1:n,1:n))


% Set b. Note the use of b is determined by rhstype settign
%    b(1:n)=randn(n,1); %  Replace with yours.
     b=[1:n]';


  elseif (problemtype == 2) || (problemtype == 4)
% standard uniform on (0,1) shifted to (dmin,dmax) set earlier
    strdmin=string(dmin);
    strdmax=string(dmax);
    temptitle = ["uniform distribution with shifted to (" strdmin "," strdmax ")"];
    outputtitle = join(temptitle);
    Temp = dmin.*ones(n,n) + (dmax-dmin).*(rand(n,n)); 
    b(1:n) = dmin.*ones(n,1) + (dmax-dmin).*(rand(n,1)); 
  elseif (problemtype == 3) || (problemtype == 5)
%  normal distribution with  mean and standard deviation set
%  to tempmean and tempstdev set earlier
    Temp  = tempstdev.*randn(n,n) + tempmean.*ones(n,n);
    b(1:n) = tempstdev.*randn(n,1) + tempmean.*ones(n,1);
    strmean=string(tempmean);
    strstdev=string(tempstdev);
    temptitle = ["normal element distribution with mean/deviation " strmean "and" strstdev];
    outputtitle = join(temptitle);
  else
    fprintf(1,'\n**************************************** \n')
    fprintf(1,'Undefined problem number   = %g \n',problemtype)
    fprintf(1,'**************************************** \n')
    return; % terminate the tester and return to command prompt
  end

  if (problemtype == 2) || (problemtype == 3)
%
% Create L and U before rescaling
%

  L = eye(n,n)+ tril(Temp,-1); % L is strict lower part of 
                               % the array of random numbers
                               % and a diagonal of all 1
  U = triu(Temp,1);  % U is strict upper part of 
                     % the array of random numbers. 
                     % The diagonal is set below.

%
%
%  Scale offdiagonal elements of L to control the conditioning.
%  This still lets elements with magnitude near 1 appear below the diagonal .
%  So it can be strengthened if the conditioning becomes a
%  problem for large n.  Note there is a scalar (3 now)
%  that can be adjusted to set target magnitude of rescaled maximum elements.
%  This is one of many strategies possible. Note that if the elements
%  generated in Temp are controlled in magnitude the
%  rescaling here and in U later is not needed.					       
%
      for i=2:n
       for j=1:i-1
        if abs(L(i,j)) > 1 
           lscale=max(max(abs(L(i,1:i-1))),max(abs(L(j+1:n,j))));
           lscale = max(lscale, abs(L(i,j)));
           L(i,j) = L(i,j)/(3.0*lscale);
        end
       end
      end

%
%
%  Set the diagonal of U and scale offdiagonal elements of U
%  to control the conditioning. Note there is a scalar (3 now)
%  that can be adjusted to set target magnitude of rescaled maximum elements
%

     udiag = umin.*ones(n,1) + (umax-umin).*(rand(n,1));
     udiag = udiag(1:n).*sign(randn(n,1)); % set sign randomly
     U=U+diag(udiag);
     for i=1:n-1
       for j=i+1:n
          if abs(U(i,j)) > 1 
             uscale=max(max(abs(U(i,i+1:n))),max(abs(U(1:j-1,j))));
             uscale=max(uscale,abs(U(i,j)));
             U(i,j) = U(i,j)/(3.0*uscale);
          end
       end
     end

end % end of form L and U for problemtype 2 or 3

%
%   Set A = LU or A = Temp depending on problem type
%

  if (problemtype == 2) || (problemtype == 3)
     Atrue(1:n,1:n) = L*U;     
  elseif (problemtype == 1) || (problemtype == 4) || (problemtype == 5)
     Atrue(1:n,1:n) = Temp(1:n,1:n);

  else
    fprintf(1,'\n**************************************** \n')
    fprintf(1,'Undefined problem number   = %g \n',problemtype)
    fprintf(1,'**************************************** \n')
    return; % terminate the tester and return to command prompt
  end

%
% b is used as the right-hand side vector in Ax=b
% as the default.
% Set it based on rhstype defined earlier.
%
    if (rhstype == 1) 
       xtrue(1:n)=b(1:n);  % b is used as the true solution
                          % compute new right-hand side vector b
       b(1:n) =  Atrue*xtrue(1:n);
    elseif (rhstype == 2 ) % xtrue is defined here and a new b computed
       xtrue(1:n)= ones(n,1)*10.0;
       b(1:n) =  Atrue*xtrue(1:n);
    else % b is used as the right-hand side and 
         % true sollution is generated using the built-in
         %  Matlab solve operator and store the result in
         %  a work array ztrue not involved in your subroutines.
       xtrue(1:n) = Atrue\b(1:n);
    end  % end of rhs setting if

%
% Put Atrue into the array M.
% L and U will appear in  simple compressed
% form in array M after the LU routine returns.
% Note the 0s of L and U are not stored and
% the 1's of the diagonal of L are note stored in the
% work array M(mrdim,mcdim).
% The vectors b, xcomp, and xtrue are already in work arrays.
%
    M(1:n,1:n) = Atrue(1:n,1:n);

%
% Run test number itest for the current value of n
%
    Prow=[1:n]';
    Pcol=Prow;
    [M,Prow,Pcol]=LUdense(M,n,routinenum,Prow,Pcol);
%    
%     f = Prow*b, solve Ly=f, solve Uz=y, x = Pcol*z
%
      f = b(Prow(1:n));
      ycomp = zeros(n,1);
      [ycomp]=Lvsolve_row(M,f,ycomp,n);
      [xcomp]=Uvsolve_row(M,ycomp,xcomp,n);
      temp = xcomp(Pcol(1:n));
      xcomp(1:n) = temp(1:n);
%
%    compute the product of the computed L and U

     Acomp = LUmult_row(M,Acomp,n);
%
%   Compute factorization Prow*Atrue*Pcol=LU; M(1:n,1:n)=Atrue on input
%   M(1:n,1:n)=L and U compressed form  on output
%   Prow and Pcol have general permutation forms (not incremental)
%
% Row permutation on each column in turn
  for iperm=1:n
    vtempcol(1:n,1) = Atrue(Prow(1:n),iperm);
    Atrue(1:n,iperm)=vtempcol(1:n,1);    
  end
% Column permutation on each row in turn
  for iperm=1:n
    vtemprow(1,1:n) = Atrue(iperm, Pcol(1:n));
    Atrue(iperm,1:n)=vtemprow(1,1:n);    
  end
%
% Compute error between computed and true x
%
    ErrorVector = xcomp(1:n) - xtrue(1:n);
%
%   Compute the norms for this problem and save with the others with this n
%
%   Prow Atrue Pcol (now in Atrue) is compared to Acomp
%   A new Atrue is generated for each problem
    Error2Norm(itest) = norm(ErrorVector,2);
    xtrue2Norm(itest) = norm(xtrue,2);  % xtrue is an n element vector
    RelError(itest) = Error2Norm(itest)/xtrue2Norm(itest);
    Akappa = max(Akappa,cond(Atrue));
    AError2Norm(itest) = norm(Atrue(1:n,1:n)-Acomp(1:n,1:n),2);
    Atrue2Norm(itest) = norm(Atrue,2);  % xtrue is an n element vector
    ARelError(itest) = AError2Norm(itest)/Atrue2Norm(itest);
% update running maximum relative error over test problems for all n
    maxrelerror = max(maxrelerror,RelError(itest));
    Amaxrelerror = max(Amaxrelerror,ARelError(itest));
  end % end of itest loop
%___________________________________________________________________________________
if outputtype == 1
  for iouts = 1:ntests
    fprintf(1,'n = %g Error 2-norm = %5.2e  xtrue 2-norm = %5.2e Relative Error = %5.2e\n',n,Error2Norm(iouts), xtrue2Norm(iouts), RelError(iouts))
    fprintf(1,'n = %g AError 2-norm = %5.2e  Atrue 2-norm = %5.2e A Relative Error = %5.2e\n',n,AError2Norm(iouts), Atrue2Norm(iouts), ARelError(iouts))
  end % output loop
end %  end outputtype ==1 if
meanforthisn(experimentindex) = mean(RelError(1:ntests));
maxforthisn(experimentindex) =  max(RelError(1:ntests));
Ameanforthisn(experimentindex) =  mean(ARelError(1:ntests));
Amaxforthisn(experimentindex) =  max(ARelError(1:ntests));
%___________________________________________________________________________________
if outputtype == 2
 fprintf(1,'For %g problems with n = %g  Mean Relative Error = %5.2e Maximum Relative Error = %5.2e\n',ntests,n, meanforthisn(experimentindex), maxforthisn(experimentindex))
 fprintf(1,'For %g problems with n = %g  Mean A Relative Error = %5.2e Maximum A Relative Error = %5.2e\n',ntests,n, Ameanforthisn(experimentindex), Amaxforthisn(experimentindex))
end

end % end of n loop
%___________________________________________________________________________________
% print overall maximum relative error
fprintf(1,'\n\n\nMaximum Relative Error over all tests of all sizes = %5.2e\n',maxrelerror);
fprintf(1,'\n\n\nMaximum A Relative Error over all tests of all sizes = %5.2e\n',Amaxrelerror);

% print overall maximum condition number
 fprintf(1,'\n\n\n Maximum condition number over all tests of all sizes = %5.2e\n',Akappa);
%___________________________________________________________________________________
if outputtype == 3
% protect against -infinity in log vectors
tempmean(1:nexperiments)=max(-17,log10(meanforthisn(1:nexperiments)));
tempmax(1:nexperiments)=max(-17,log10(maxforthisn(1:nexperiments)));
figx=figure; % starts figure window for x error data
plot([nmin:ndelta:nmax],tempmean(1:nexperiments),'-.k')  % plot 
hold on
plot([nmin:ndelta:nmax],tempmax(1:nexperiments),'-r')  % plot 
xlabel('dimension of matrices in set')
ylabel('log base 10 of relative x error max and mean per set')
title(outputtitle);

%_________ debug stuff
% debug printing of mean and max data
% fprintf(1,'\n\n\nMean Relative Errors vs n = %5.2e\n',meanforthisn(1:nexperiments));
% fprintf(1,'\n\n\nMaximum  Relative Errors vs n = %5.2e\n',maxforthisn(1:nexperiments));
hold off



% protect against -infinity in log vectors
tempmean(1:nexperiments)=max(-17,log10(Ameanforthisn(1:nexperiments)));
tempmax(1:nexperiments)=max(-17,log10(Amaxforthisn(1:nexperiments)));
figA=figure;  % starts new figure window for the A error data
plot([nmin:ndelta:nmax],tempmean(1:nexperiments),'-.k')  % plot 
hold on
plot([nmin:ndelta:nmax],tempmax(1:nexperiments),'-r')  % plot 
xlabel('dimension of matrices in set')
ylabel('log base 10 of relative x error max and mean per set')
title(outputtitle);


%_________ debug stuff
% debug printing of mean and max data
% fprintf(1,'\n\n\nMean  Relative Errors vs n = %5.2e\n',meanforthisn(1:nexperiments));
% fprintf(1,'\n\n\nMean A Relative Errors vs n = %5.2e\n',Ameanforthisn(1:nexperiments));
% fprintf(1,'\n\n\nMaximum   Relative Errors vs n = %5.2e\n',maxforthisn(1:nexperiments));
% fprintf(1,'\n\n\nMaximum  A Relative Errors vs n = %5.2e\n',Amaxforthisn(1:nexperiments));
hold off

%___________________________________________________________________________________

end % end of outputtype = 3 code


% outputtype = 4 special output for specific problems

if outputtype == 4
disp('output type = 4')

% This is used for structure checking for certain problemtype == 1 problem
% L, U, and PA are constructed by extraction for display and validation
% for all of the problems in the assignment. Note Atrue has
% already been permuted for comparison to Acomp in the error computation

L(1:n,1:n) = eye(n,n)+ tril(M(1:n,1:n),-1); % L is strict lower part of 
                               % the array of random numbers
                               % and a diagonal of all 1
U = triu(M(1:n,1:n));  % U is strict lower part of 

disp('M data structure after factorization')
M(1:n,1:n)

disp('L ')
L(1:n,1:n)

disp('U')
U(1:n,1:n)

disp('Permutation')
Prow(1:n)

disp('Prow*A')
Atrue(1:n,1:n)

disp('L*U')
PA=L(1:n,1:n)*U(1:n,1:n)

Aabs(1:n,1:n)=zeros(n,n);
Aabs(1:n,1:n) = absLUmult_row(M,Aabs,n);

disp('Computed growth factor')
gamma = norm(Aabs(1:n,1:n),2)/norm(Atrue(1:n,1:n),2)
end
