matrix - How to get the Q from the QR factorization output? -

dgeqrf , sgeqrf lapack return q part of qr factorization in packed format. unpacking seems require o(k^3) steps (k low-rank products), , doesn't seem straightforward. plus, numerical stability of doing k sequential multiplications unclear me.

does lapack include subroutine unpacking q, , if not, how should it?

yes, lapack indeed offers routine retrieve q elementary reflectors (i.e. part of data returned dgeqrf), called dorgqr. describtion:

*  dorgqr generates m-by-n real matrix q orthonormal columns, *  defined first n columns of product of k elementary *  reflectors of order m * *        q  =  h(1) h(2) . . . h(k) *  returned dgeqrf. 

a complete calculation of q , r a using c-wrapper lapacke (a fortran adaption should straight forward) :

void qr( double* const _q, double* const _r, double* const _a, const size_t _m, const size_t _n) {     // maximal rank used lapacke     const size_t rank = std::min(_m, _n);       // tmp array lapacke     const std::unique_ptr<double[]> tau(new double[rank]);      // calculate qr factorisations     lapacke_dgeqrf(lapack_row_major, (int) _m, (int) _n, _a, (int) _n, tau.get());      // copy upper triangular matrix r (rank x _n) position     for(size_t row =0; row < rank; ++row) {         memset(_r+row*_n, 0, row*sizeof(double)); // set starting zeros         memcpy(_r+row*_n+row, _a+row*_n+row, (_n-row)*sizeof(double)); // copy upper triangular part lapack result.     }      // create orthogonal matrix q (in tmpa)     lapacke_dorgqr(lapack_row_major, (int) _m, (int) rank, (int) rank, _a, (int) _n, tau.get());      //copy q (_m x rank) position     if(_m == _n) {         memcpy(_q, _a, sizeof(double)*(_m*_n));     } else {         for(size_t row =0; row < _m; ++row) {             memcpy(_q+row*rank, _a+row*_n, sizeof(double)*(rank));         }     } } 

it's piece of own code, removed checks improve readability. productive use want check input valid , mind return values of lapack calls. note input a destroyed.