Exercise "Linear Equations"

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. In garbage-collected implementations freeing is not needed as it is done automaticall by the garbage collector.
• Fortran: one- and two-dimensional arrays as in LAPACK;
• C, C++: arrays and jagged arrays: double* v; double** m (no range checking, no storing of sizes);
• C, C++: gsl_vector and gsl_matrix structures from GSL library -- perhaps, the best choice if you don't want to build your own library;
• C++: vector and matrix class from Armadillo library;
• Python: python lists, numpy arrays;
• Java: simply 1d and 2d arrays:

  double[][] m = new double[7][2];
  System.out.println(m.length);
  System.out.println(m[0].length);
  

  or a simple matrix class.

For QR-factorization you have the following choices:

• Procedurial style:
  • Create a function void qrdec(matrix* A, matrix* R) that performs in-place Gram-Schmidt orthogonalization of matrix A (that is, A turns into Q) and fills-in the matrix R.
  • Create a function void qrbak(matrix* Q, matrix* R, vector* b, vector* result) which solves the equation QRx=b by back-substitution and puts the result in vector result.
  • Create a function double absdet(matrix* R) which calculates the absolute value of the determinant.
  • Create a function void inverse(matrix* Q, matrix* R, matrix* AI) which calculates the inverse into matrix AI.
• Object-oriented style:
  • You can simply implement static functions in the procedurial style.
  • You can try do some object-oriented stuff, like create an object QR with member variables matrix Q, R; and the constructor QR(matrix& A) which first copies A into Q and then performs the in-place orthogonalization. The functions from the procedurial style become member functions
    void qrbak(vector& b, vector& result)
double absdet()
void inverse(matrix& AI)

Tasks

1. (6 points) Gram-Schmidt orthogonalization

   • Implement a function which calculates QR-decomposition of a given matrix by modified Gram-Schmidt orthogonalization.
   • Implement a function which, given the matrices Q and R from the QR-decomposition and the right-hand-side vector b, solves the system of linear equations Ax=b by back-substitution.
   • Implement a function which, given the matrices Q and R from the QR-decomposition, calculates the (absolute value of the) determinant of matrix A.
   • Implement a function which, given the matrices Q and R from the QR-decomposition of a matrix A, calculates its inverse.
   • 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) Given's rotation

   • 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. The backsubstitution routine should not build the Q-matrix either but rather apply (smartly) the stored individual Given's rotations. This should save you a whole matrix to store as compared to your first exercise.

   • Compare the speeds of these two implementations.

3. (1 point)

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

   • Or, at your choice, Cholesky-decomposition using Cholesky-Banachiewicz or Cholesky-Crout algorithms.

4. (0 points)
   • Compare the speed of your own QR implementation with your favourite library routine, e.g.
     • (C,C++) GSL routine 'gsl_linalg_QR_decomp' (an example of its usage is here);
     • (C++) Armadillo routine 'qr';
     • (Fortran) LAPACK routine 'dgeqrf';
     • (Python) SciPy routine Sci.linalg.qr;
     by aplying them, for example, to random matrices of different sizes.

     Hint: the total number of CPU-seconds used by a program my_prog can be determined by the POSIX time utility. The command

     
     \time --format "%U" --append --output=time.txt ./my_prog > /dev/null 
     

     or, in short

     
     \time -f "%U" -ao time.txt ./my_prog > /dev/null 
     

     runs the program (disregarding the standard output) and appends the number of consumed CPU-seconds to the file time.txt. Without the --append option (or, in short, -a) the file time.txt gets overwritten. Backslash here is needed to run the actual utility rather than built-in bash command (with similar possibilities, actually).

     An example of the usage of time can be found here (see the makefile).

   • C, C++:

     Implement a garbage-collected version of the vector/matrix library and compare the memory consumption of the garbage-collect vs. malloc'ed versions. The garbage collector is installed in the proper place on molveno and in /usr/users/fedorov/include, /usr/users/fedorov/lib on lifa.