Exercise "QR-decomposition"

• Fortran: two-dimensional arrays as in LAPACK;
• C: jagged arrays: int m, n; double** A;
• C: struct gsl_matrix from GSL library;
• C++: matrix class from Armadillo library;
• Python: numpy.matrix;

• Procedurial style:
  • Create a function void qrdec(mat A, mat R) that performs in-place Gram-Schmidt orthogonalization of matrix A and fills-in the matrix R.
  • Create a function void qrbak(mat Q, mat R, vec b, vec result) which solves the equation QRx=b by back-substitution and puts the result in vector result.
  • Create a function double absdet(mat R) which calculates the absolute value of the determinant.
  • Create a function void inverse(mat Q, mat R, mat AI) which calculates the inverse into matrix AI.
• Object-oriented style:
  • Create an object QR with member variables mat Q, R; and the constructor QR(mat A) which first copies A into Q and then performs the in-place orthogonalization.
  • The functions from the procedurial style become now member functions
    void qrbak(vec b, vec result)
double absdet()
void inverse(mat AI)

1. (6 points)
   • 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 data and checks that everything works as intended.
2. (3 points)
   Compare the speed of your QR-routine with a 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 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 "time = %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 (or, in short, -a) option 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).

3. (1 point) Check that the time to perform a QR-decomposition of an (n rows)x(m cols) matrix scales as mn².
4. (0 points) Implement LU-decomposition and back-substitution.