\documentclass[11pt,twoside]{article}
\newcommand{\HW}{}
\newcommand{\HWnumb}{2}
\newcommand{\HWdate}{Wednesday April 22}
\input{hw0}
\input{macros}

\def    \CG{\hbox{\small     CG}}
\def   \PCG{\hbox{\small    PCG}}
\def \CRAIG{\hbox{\small  CRAIG}}
\def  \LSQR{\hbox{\small   LSQR}}
\def\MINRES{\hbox{\small MINRES}}
\def\MINRESQLP{\hbox{\small MINRES-QLP}}
\def\SYMMLQ{\hbox{\small SYMMLQ}}
\def\MA#1#2{\hbox{\small MA#1#2}}

\thispagestyle{empty}

\begin{enumerate}
\item For the symmetric Lanczos process,
prove by induction that with exact arithmetic, the matrix
$V_k = \pmat{v_1 \ v_2 \ \ldots \ v_k}$
has orthonormal columns: \hbox{$V_k^T V_k\drop = I$}.
(Hence, $V_k$ is an orthonormal basis for the Krylov subspace
$\mathcal{K}(A,b,k) = \mathrm{range} \{b, Ab, \ldots, A^{k-1}b \}$.)

\item \url{http://www.stanford.edu/class/msande318/homework/}
and \url{http://www.stanford.edu/class/msande318/matlab/}
are the homework websites.  From links there,
download files \texttt{Si2.mat}%
\footnote{These questions use PARSEC/Si2, a symmetric indefinite sparse
  matrix $A$ of size $n=769$ from a collection of sparse-matrix data
  and software maintained by Tim Davis (University of Florida):
  \url{http://www.cise.ufl.edu/research/sparse/}.}
and \texttt{CGtest3.m}.
Script \texttt{CGtest3} shows how to load an $n \times n$ matrix $A$ into
\Matlab{} ($n=769$).  The first eigenvalue of $A$ is negative
($\lambda_1 = -0.38441$).  The remaining eigenvalues range from
$\lambda_2 = 0.24221$ to $\lambda_n = 41.3813$, so $\cond(A) \approx 170$.

The matrix $P = A + \sigma I$ with $\sigma = 0.4$ is positive
definite, and $\cond(P) \approx 2679$.
Define $x_j = 1/j$ ($j=1 \To n$) and $b = Px$.

\item[(a)]
Solve $Px = b$ using some of the
iterative solvers provided in \Matlab: pcg, minres, and symmlq.
Use the parameters \verb|tol = 1e-12| and \verb|maxit = 500|.
For each method, plot the residuals $\norm{r_k}$ for each iteration $k$.
Script \texttt{CGtest3} already does this in producing figure(1).
Add suitable \{xlabel, ylabel, legend\} and give a table of results.

\item[(b)]
The residuals for each solver are significantly different
at around iteration 33.  Did all solvers terminate at much the same point?
Should we expect them to be more different?

\item[(c)]
\verb|type pcg| allows you to see \Matlab's implementation of CG.
For each solver, find which stopping rule was used for
the results shown in figure(1).


\item The matrix $A$ is indefinite.  With the same $x = [1/j]$, define $b = Ax$.

\item[(a)]
Solve $Ax = b$ using pcg, minres, and symmlq.  
Solve $A^2x = Ab$ using pcg, and $Ax = b$ using lsqr.
Script \texttt{CGtest3} already does this in producing figure(2).
Again, add suitable \{xlabel, ylabel, legend\} and give a table of results.

\item[(b)]
Note that pcg on $Ax = b$ terminates early.  What happened?

\item[(c)]
As we see, when $A$ is reasonably well-conditioned, lsqr performs much the same as
pcg on $A^2x = Ab$.  But what is the plot for pcg really showing?
(It's not $\norm{r_k} = \norm{b - Ax_k}$.)

\item
The script plots the eigenvalues $\lambda(A)$ in figure(3).
Does there seem to be any clustering of the eigenvalues?
A better picture is given in figure(4).  Briefly describe
what the script is showing in figure(4).
Does it explain why the symmetric solvers needed
significantly fewer than $n$ iterations
to solve $Ax = b$?

\end{enumerate}

\newpage

\begin{small}
\begin{verbatim}
% CGtest3.m is a script for comparing {pcg, minres, symmlq}
% on UF sparse matrix number 1360 (PARSEC/Si2, n=769).
% See http://www.cise.ufl.edu/research/sparse/matrices/ (Tim Davis).
% The matrix A is an example from quantum chemistry.
% It is symmetric indefinite with lambda_min = -0.3844.
%
% 08 Apr 2008: First version for MS&E 318 Homework 2.
%----------------------------------------------------------

load Si2.mat
A     = Problem.A;
[m,n] = size(A);
x     = 1./(1:n)';

%----------------------------------------------------------
fprintf('\nPOSITIVE DEFINITE SYSTEM')
fprintf('\n                flag  iter   relres    error\n')
P     = A + 0.4*speye(n);
b     = P*x;
tol   = 1e-12;
maxit = 600;

[xL,flagL,relresL,iterL,resvecL] = symmlq(P,b,tol,maxit);
[xC,flagC,relresC,iterC,resvecC] = pcg   (P,b,tol,maxit);
[xM,flagM,relresM,iterM,resvecM] = minres(P,b,tol,maxit);

errL  = norm(xL-x,inf);
errC  = norm(xC-x,inf);
errM  = norm(xM-x,inf);
fprintf(' SYMMLQ   Px = b%4g %5g %8.1e %8.1e  r\n', flagL,iterL,relresL,errL)
fprintf(' CG       Px = b%4g %5g %8.1e %8.1e  b\n', flagC,iterC,relresC,errC)
fprintf(' MINRES   Px = b%4g %5g %8.1e %8.1e  g\n', flagM,iterM,relresM,errM)

figure(1)
hold off;  plot(log10(resvecL),'ro')
hold  on;  plot(log10(resvecC),'b-')
hold  on;  plot(log10(resvecM),'g.')


%----------------------------------------------------------
fprintf('\nINDEFINITE SYSTEM')
fprintf('\n                flag  iter   relres    error\n')
Afun  = @(x) A*x;       b  = A*x;      % Treat  A  as a function
Afun2 = @(x) A*(A*x);   b2 = A*b;      % Treat A*A as a function

[xL,flagL,relresL,iterL,resvecL] = symmlq(Afun ,b, tol,maxit);
[xC,flagC,relresC,iterC,resvecC] = pcg   (Afun ,b ,tol,maxit);
[xM,flagM,relresM,iterM,resvecM] = minres(Afun ,b ,tol,maxit);
[xN,flagN,relresN,iterN,resvecN] = pcg   (Afun2,b2,tol,maxit);
[xS,flagS,relresS,iterS,resvecS] = lsqr  (A    ,b ,tol,maxit);

errL  = norm(xL-x,inf);
errC  = norm(xC-x,inf);
errM  = norm(xM-x,inf);
errN  = norm(xN-x,inf);
errS  = norm(xS-x,inf);
fprintf(' SYMMLQ   Ax = b%4g %5g %8.1e %8.1e  r\n', flagL,iterL,relresL,errL)
fprintf(' CG       Ax = b%4g %5g %8.1e %8.1e  b\n', flagC,iterC,relresC,errC)
fprintf(' MINRES   Ax = b%4g %5g %8.1e %8.1e  g\n', flagM,iterM,relresM,errM)
fprintf(' CG     A^2x =Ab%4g %5g %8.1e %8.1e  k\n', flagN,iterN,relresN,errN)
fprintf(' LSQR     Ax = b%4g %5g %8.1e %8.1e  m\n', flagS,iterS,relresS,errS)

figure(2)
hold off;  plot(log10(resvecL),'ro')
hold  on;  plot(log10(resvecC),'b-')
hold  on;  plot(log10(resvecM),'g.')
hold  on;  plot(log10(resvecN),'k.')
hold  on;  plot(log10(resvecS),'m-')

%return

%----------------------------------------------------------
% Plot the eigenvalues of A.
%----------------------------------------------------------
lambda = eig(full(A));
figure(3)
hold off;    plot(lambda,'b.')
xlabel('Eigenvalue number')
ylabel('\lambda(A)')
title('Eigenvalues of A')

% Try to show if the eigenvalues are clustered.
figure(4);
hold off;    plot(lambda,10*ones(n,1),'b.')
hold on;

y1    = -5;
yn    = 44.9;
step  = 2;

for y = y1:step:yn
  y2  = y1 + step;
  nlambda = length( find(lambda>y1 & lambda<=y2) );
  bar(y1+0.5*step,nlambda)
  y1  = y2;
end

xlabel('\lambda(A)')
ylabel('No. of \lambda(A)')
title('Distribution of eigenvalues of A')
\end{verbatim}
\end{small}

\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: