a logo for the course

Matrix Computations

David Gleich

Purdue University

Fall 2017

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Forney B124

Homework 6

Homework 6

Please answer the following questions in complete sentences in a typed, clearly prepared manuscript and submit the solution by the due date on Blackboard (Monday, October 23rd, 2017, around 4-5 am)

Problem 0: Homework checklist

Problem 1: Direct Methods for Tridiagonal systems

In class, we said the definition of a sparse matrix was one with enough zeros to take advantage of it. In this problem, we'll show how to take advantage of tridiagonal structure in a matrix.

Let be a symmetric, tridiagonal, positive definite matrix:

  1. Let is the leading principal minor of . Suppose you are given the Choleksy factorization of . Determine the structure of (hint, must be tridiagonal!).

  2. Now, show how to compute the Cholesky factorization of from in a constant number of operations.

  3. Use the result of the first two parts to write down an algorithm that computes the Cholesky factorization of in a linear amount of work (starting from scratch!)

Problem 2: Sparse matrices in Julia

If you wish to use another language, please see me

In this problem, we'll return to our favorite linear system, Poisson's equation. (I told you that we'd see a lot of this system.)

We'll assume that the boundary conditions are zero everywhere. Consequently, we'll use the form of the system where we have removed all rows that are set to zero. Some of you did this already, so you may have an easier time on this question.
For the case, we have: where and .

  1. You can change an entry of a sparse matrix in Julia with . Write a Julia routine to create the matrix above for a general by changing each element one at a time.
A = spzeros((N-1)*(N-1),(N-1)*(N-1)); # this creates an all-zero sparse matrix
<fill in the rest>
  1. As shown in class, another way to create a sparse matrix is to use the coordinate form. Create the arrays I,J,andV for the matrix above. Hint there are non-zeros. Use this information to initialize the arrays I,J, and V. Write a routine to create the matrix above for a general using the coordinate form
nz = (N-1)^2 + 4*(N-1)*(N-2);
I = zeros(nz);
J = zeros(nz);
V = zeros(nz);
<fill in the rest>
A = sparse(i,j,v,(N-1)*(N-1),(N-1)*(N-1));
  1. Which way is faster? Try creating matrices with .

  2. Let . Is solving the linear system using backslash faster or slower with a sparse matrix? Does it depend on the value of ? Investigate and report.

  3. Again let . Implement the Jacobi method for this problem. You should start from the iterate , so start with the guess of the right hand side. How many iterations does it take to converge to a relative residual of ? Use the same values of from part 3.

  4. Does the number of iterations of your Jacobi solver depend on the the right hand side? Try comparing , , , . Use the same values of from part 3.

Problem 3: Fun with sparse matrices

In Julia, calling y = A*x on a sparse matrix A, calls a special function to perform this multiplication faster. To generate a sparse random matrix in julia, you can call A = sprand(M,N,d) where M and N are the dimensions of A, and d is the density of A. 1. Generate a random sparse matrix A and transform it to dense matrix B. Use a dense vector x of suitable dimensions and compute A*x and B*x. Vary M, N, and d and report performance changes.

  1. Generate a random square sparse matrix A, and a dense vector x of suitable dimension and compute A*x and A'*x. Is one faster than the other? Vary the dimensions and the density and report your results. Provide an analysis of your results.

  2. Write a sparse matrix vector routine where both, the matrix and the vector are sparse. You can generate a sparse vector in Julia using x = sprand(N,d), where N is the dimension and d is the density. Test your method and report performance results.

function spmatvec{T,F}(A::SparseMatrixCSC{T,Int64},x::SparseVector{F,Int64})
       y = <fill in>
       for j = <fill in>
           for k = A.colptr[j]:(A.colptr[j+1]-1)
               <fill in>
       return y