Exercise "Linear Equations"

Introduction

To work with vectors and matrices you can use

• your own structures and the associated functions;
• Fortran: one- and two-dimensional arrays as in LAPACK;
• C, C++: arrays and jagged arrays: double* v; double** m;
• C, C++: gsl_vector and gsl_matrix structures from GSL 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 vector b of the right-hand-sides, 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.

   Never tried it myself yet, but it looks like it should work ;)

3. (1 point) The speed

   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).

4. (0 points)
   • All:

     Implement LU-decomposition, back-substitution, and inverse.

   • 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.