Please answer the following questions in complete sentences in a typed, clearly prepared manuscript. This homework is due by Monday, December 4th, 2017, 4:30am. However, everyone is given an automatic, full-points extension until Wednesday, December 6th, 2017 at 4am if you want your homework graded by December 7th in class. This is designed to give you maximum flexibility in order to plan your time. We highly suggest sticking to the original due date, but there is absolutely no penalty for taking more time.
Please identify anyone, whether or not they are in the class, with whom you discussed your homework. This problem is worth 1 point, but on a multiplicative scale.
Make sure you have included your source-code and prepared your solution according to the most recent Piazza note on homework submissions.
It's based on an assignment from Michael Saunders' iterative methods class.
Download the poisson2D.mat
file from Tim Davis's sparse matrix repository.
From Julia, you can use the following code
to setup and from the linear system.
using MAT
file = matread("poisson2D.mat")
A = file["Problem"]["A"]
A = (A+A')./2 # make it symmetric
b = file["Problem"]["b"]
b = vec(b)
The matrix A is symmetric-positive definite. But we are going to look at solving where . We'll declare a method converged if the relative residual has -norm less than .
Plot the relative residuals for 100 iterations of the MINRES (minres
from IterativeSolvers.jl
)
method and the GMRES (gmres
from IterativeSolvers.jl
) method restarted every 30 iterations on a semilogy scale. Describe what you observe, in particular, do the methods converge?
For convenience, here is some code to help you get started:
using IterativeSolvers
x = minres(A,b,tol = mytol, maxiter = mymaxiter)
x = gmres(A,b,restart = myrestart, tol = mytol, maxiter = mymaxiter)
If you use an incomplete cholesky ichol
or incomplete LU ilu
preconditioner with GMRES, it will converge. How many iterations does
it take? Does MINRES benefit from an incomplete Cholesky preconditioner?
We provide you with the ichol
and ilu
code.
Download the file hw11_functions.jl
from here http://www.cs.purdue.edu/homes/dgleich/cs515-2017/julia/hw11_functions.jl. The IterativeSolvers package still has some rough edges. We will help you set that up. Here is the code to help you get started
using ItertaiveSolvers
include("hw11_functions.jl")
Lchol = ichol(A)
M = MyPreconditionerType(Lchol)
x = gmres(A,b,restart = 30, tol = mytol, maxiter = mymaxiter,Pl=M)
For MINRES, it turns out that the minres
function in IterativeSolvers still doesn't accept preconditioners. We use another package to test minres
.
# note that the Krylov package is not a registered Julia package
# yet but you can still add it as follows
Pkg.clone("https://github.com/JuliaSmoothOptimizers/Krylov.jl.git")
using Krylov
using LinearOperators
L = ichol(A)
P1 = opInverse(LowerTriangular(L))
Mprecond = P1'*P1
x = Krylov.minres(A,b,M=Mprecond,itmax=mymaxiter)[1]
Show the eigenvalues of the matrix before and after preconditioning. Do you observe any clustering?
A real-symmetric matrix has an eigenvalue of of multiplicity , while the rest of the eigenvalues are all in absolute value. Describe an algorithm for finding an orthonormal basis of the 8-dimensional eigenspace corresponding to the domininant eigenvalue.