Exercise "QR-decomposition"

• 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** A;
• C, C++: gsl_vector and gsl_matrix structures from GSL library;
• C++: vector and matrix class from Armadillo library;
• Python: numpy;

• Procedurial style:
  • Create a function void qrdec(matrix* A, matrix* R) that performs in-place Gram-Schmidt orthogonalization of matrix A 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:
  • 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 now member functions
    void qrbak(vector& b, vector& result)
double absdet()
void inverse(matrix& 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 matrices and checks that everything works as intended.
2. (3 points)
   • Implement (if you have not done it in part "A") your own vector and matrix library and for the qrdec, qrbak, inverse functions. An example of such library can be found here. The header file vector.h might looke like this,

     #ifndef HAVE_VECTOR_H
     typedef struct {int size, stride; double *data;} vector;
     vector* vector_alloc(int n);
     double vector_get(vector* v, int i);
     void vector_set(vector* v, int i, double x);
     /* whatever functions you deem necessary go here... */
     void vector_free(vector* v);
     #define HAVE_VECTOR_H
     #endif
     

     The header file matrix.h might looke like this,

     #ifndef HAVE_VECTOR_H
     #include<vector.h>
     #endif
     #ifndef HAVE_MATRIX_H
     typedef struct {int size1, size2; double *data;} matrix;
     matrix* matrix_alloc (int n1, int n2);
     double matrix_get(matrix* m, int i, int j);
     void matrix_set(matrix* m, int i, int j, double x);
     /* whatever functions you deem necessary go here... */
     void matrix_free(matrix* a);
     #define HAVE_MATRIX_H
     #endif
     

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

3. (1 points)

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

4. (0 point)

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