% (c) 2004/2005 Universitaet Stuttgart % Chair 'Numerische Mathematik fuer Hoechstleistungsrechner' % http://www.ians.uni-stuttgart.de/nmh % % ---------------------------------------------------------- % Permission to use this exercise for non-profit educational % purposes is hereby granted provided that this copyright % notice is included in all copies or substantial portions % of the LaTeX files. % ---------------------------------------------------------- \begin{exercise}{Unsymmetric systems}{tutorial}{10}{indir09\_eng} Let the following Matlab code be given (you may download it from the course Web site). The code loads one of the Harwell-Boeing test matrices. Then, a comparison plot between the CGS and BiCG methods is generated for solving the system $Ax=f$, $f = (1, \ldots, 1)^T$. {\small \begin{verbatim} % Script file for comparing different solvers for the west0479 test % matrix of the Harwell/Boing sparse matrix test set. data = load('west0479'); A = data.west0479; f = ones(size(A,1),1); methods = {'CGS', 'BI-CG'}; nmethods = length(methods); linestyles = {'rs-','bo-'}; domethods = 1:2; tol = 1e-6; maxit = 20; for mk = 1:nmethods; m = domethods(mk); disp(['Method: ' methods{m} ', maxit = ', num2str(maxit)]); switch m case 1 % CGS [x, flag, res, it, resvec] = cgs(A,f,tol,maxit); case 2 % BI-CG [x, flag, res, it, resvec] = bicg(A,f,tol,maxit); end disp(['Computed relative residual: ' num2str(res)]); h = semilogy(0:(length(resvec)-1),resvec/norm(f), linestyles{mk}, ... 'LineWidth', 1); hold on; end hold off; xlabel('number of iterations'); ylabel('relative residual'); legend(methods(domethods),1); \end{verbatim} } \begin{enumerate} \item Run the script in Matlab. One observes that the system cannot be solved successfully, although the matrix is regular. Suggest a possible reason for this behavior and verify it. \item Add the BiCGStab and the GMRES method to the comparison by modifying the provided code. Do these methods converge? Try to alter the maximum number of iterations and the restart value of the GMRES method. \item We will now examine the possibility of adding an incomplete LU-decomposition as preconditioner to the equation system. Use the Matlab command \texttt{help luinc} to learn more about incomplete LU decomposition in Matlab. Using a drop tolerance of $10^{-6}$, determine the condition of the matrix of the preconditioned system $(L_\text{inc}U_\text{inc})^{-1}Ax = (L_\text{inc}U_\text{inc})^{-1}f$. \item Modify the provided code so it an ILU preconditioner generated with a drop tolerance of $10^{-6}$. How many iterations are required by the four methods to reach a relative residual of $10^{-7}$? Would you user even smaller drop tolerances? Explain why or why not. \end{enumerate} \end{exercise}