How to get the Q from the QR factorization output?

Question

DGEQRF and SGEQRF from LAPACK return the Q part of the QR factorization in a packed format. Unpacking it seems to require O(k^3) steps (k low-rank products), and doesn't seem to be very straightforward. Plus, the numerical stability of doing k sequential multiplications is unclear to me.

Does LAPACK include a subroutine for unpacking Q, and if not, how should I do it?

Answer 1

Yes, LAPACK indeed offers a routine to retrieve Q from the elementary reflectors (i.e. the part of data returned by DGEQRF), it is called DORGQR. From the describtion:

*  DORGQR generates an M-by-N real matrix Q with orthonormal columns,
*  which is defined as the first N columns of a product of K elementary
*  reflectors of order M
*
*        Q  =  H(1) H(2) . . . H(k)
*  as returned by DGEQRF.

A complete calculation of Q and R from A using the C-wrapper LAPACKE could look like this (a Fortran adaption should be straight forward) :

void qr( double* const _Q, double* const _R, double* const _A, const size_t _m, const size_t _n) {
    // Maximal rank is used by Lapacke
    const size_t rank = std::min(_m, _n); 

    // Tmp Array for 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 the upper triangular Matrix R (rank x _n) into 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 from 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) into 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 a piece of my own code, where I removed all checks to improve readability. For productive use you want to check that the input is valid and also mind the return values of the LAPACK calls. Note that the input A is destroyed.