Exercise "Linear Equations"

Objective

Implement functions for solving linear equations, calculating matrix inverse and matrix determinant.

Introduction

To work with vectors and matrices you can use:

• Your own structures and the associated functions, like here:

  • python/matrix
  • c/matrix
  • java/matrix

  At the bare minimum, you only need to implement four functions: matrix_allocate, matrix_set_element, matrix_get_element, matrix_free.

  On top of that you can implement other useful functions like add, times, scale, transpose, print and others.

  Building your own vector/matrix library should help you learn your language better.

• C, C++: gsl_vector and gsl_matrix structures from GSL library.
• C++: vector and matrix class from Armadillo library.
• Python: python lists, numpy arrays.

Tasks

1. (6 points) QR-decomposition by modified Gram-Schmidt orthogonalization.

   • Implement a function, say

     void qr_gs_decomp(matrix* A, matrix* R)

     that performs in-place modified Gram-Schmidt orthogonalization of matrix A (that is, A turns into Q) and fills the matrix R.

   • Implement a function, say

     void qr_gs_solve(const matrix* Q, const matrix* R, const vector* b, vector* x)

     that, given the matrices Q and R from qr_gs_decomp, solves the equation QRx=b by applying Q^T on the right-hand side b and then performing back-substitution.

   • Implement a function, say

     double qr_gs_absdet(const matrix* R)

     that, given the matrix R from qr_gs_decomp, calculates the absolute value of the determinant of matrix A.

   • Implement a function, say

     void qr_gs_invert(const matrix* Q, const matrix* R, matrix* Ainverse)

     that, given the matrices Q and R from qr_gs_decomp, calculates the inverse of the matrix A into the matrix Ainverse.

   • Implement the main function which applies the implemented functions on some matrices and checks that everything works as intended.

   • Make sure your functions work as they should on non-square matrices -- we shall need this for least-squares.

2. (3 points) QR-decomposition by Given's rotations.

   • Do the same, but this time use Given's rotation algorithm instead of Gram-Schmidt algorithm.

     Devise a smart method of storing the Q-matrix: do not actually calculate and store the Q-matrix but rather store the rotation θ-angles in the places of the corresponding zeroed elements of matrix A. Neither should the back-substitution routine build the Q-matrix explicitly, but but rather apply consecutively the stored individual Given's rotations.

3. (1 point)

   • Implement LU-decomposition (with or without partial pivoting), back-substitution, and inverse.

   • Compare the relative speed of these three implementations.

4. (0 points)

   • Implement Cholesky-decomposition using Cholesky-Banachiewicz or Cholesky-Crout algorithms.

   • Make up an interesting exercise on linear equations.