%---------------------------------------------------------------------------
%
% Driver script for cg.m solving the linear system Ax = b, where A is
% an nxn matrix. This sets up a symmetric and positive definite 
% (spd) matrix, calls cg, then plots the results. The two ways of setting up
% a matrix are 'curve' or 'normalequations'.
%
% It also writes out the linear system to ASCII files:
%
%   coeffmatrix: the matrix A
%   initialx   : starting vector for the algorithm
%   rhs        : the vector b
%   solutionx  : the true solution vector
%
% Each one of the files has on its first line
%   nrows ncols nnz
% where nnz is the number of nonzeros in A. To see the sparsity pattern,
% a spy plot of A is shown. Use the files as input to your CG code.
% Shifting to a sparse version is easy, and at most would require 
% changing the matrix*vector function that computes y = A*x. For now,
% all of the entries of A are written out, including ones that are
% equal to zero.
%
% Although A is n x n, both of the integers nrows and ncols are 
% required - the parallel version will have ncols = n but nrows = n/p
% (approximately, the last process with rank p-1 will be smaller/larger
% when n/p is not an integer).
%
%--------------------------------------------
%
%    Randall Bramley
%    Department of Computer Science
%    Indiana University, Bloomington
%    bramley somewhere around cs dot indiana dot (allegedly) edu                             
%
%    Started: Fri Jun  8 11:21:20 EDT 2007
%    Modified: 06 Apr 2009 to have different matrix types
%    Last Modified: Tue 07 Apr 2009, 11:45 AM to call setA for matrix set-up
%
%
%---------------------------------------------------------------------------

%----------------------------------------------------------
% Parameters. n is the order of the matrix A
%             maxits = maximum number of iterations allowed
%             tol = stopping criterion on norm(residual)
%----------------------------------------------------------
n      = 100;
maxits = 1.5*n;
tol    = 1.0e-8;

%-------------------------------------------------------------
% Choose power large to have non-convergent systems, and small
% but still positive to get quickly convergent ones. 0.5 is a
% good starting point.
%-------------------------------------------------------------
power = 0.5;

[A, b] = setA(n, 'curve', power);

xinit = randn(n, 1);

figure;
spy(A)
title('Sparsity pattern of coefficient matrix');
%---------------------------------------------------------------------
% Output the quantities, each in its own file. Must use 17 significant
% digits to get reproducibility. First line in the file for the arrays
% has the value of n and the number of nonzeros. All succeeding lines
% have three values: row number, column number, and value. So a line
% with 
%
%           i  j  valij
%
% specifies A(i,j) = valij.  A vector will have all of the values
% j = 1, and nnz = n. For a dense matrix, nnz = n^2.
%
%---------------------------------------------------------------------

t0 = clock;
writearray(A, 'coeffmatrix');
writearray(b, 'rhs');
writearray(xinit, 'initialx');
t1 = clock;
writetime = etime(t1, t0);
disp(sprintf('Time to write out the files: %e seconds', writetime));

%----------------------------------------------
% The next line  runs the Matlab version of CG:
%----------------------------------------------
[x, its, pltret] = cg(A, b, xinit, tol, maxits);

figure;
semilogy(pltret(:,1), pltret(:,2), 'b+', ...
         pltret(:,1), pltret(:,2), 'r-');

title('Residual norms from CG');
xlabel('Iteration Number');
ylabel('norm(r)');


disp(sprintf('Actual iterations: %d', its))
disp(sprintf('Final residual norm: %e', pltret(its, 2)))

writearray(x, 'solutionx');