Module_QR_Procedures

Copyright 2022 IRD

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR COMPUTING QR, LQ DECOMPOSITIONS AND RELATED DECOMPOSITIONS.

LATEST REVISION : 18/07/2022

subroutine lq_cmp ( mat, diagl, tau, use_qr )

Purpose

LQ_CMP computes a LQ factorization of a real m-by-n matrix MAT :

MAT = L * Q

where Q is orthogonal and L is lower trapezoidal (lower triangular if m<=n).

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real matrix to be decomposed.

On exit, the elements below the diagonal of the array contain the corresponding elements of L. The elements on and above the diagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors, see Further Details.

DIAGL (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix L.

The size of DIAGL must be min( size(MAT,1) , size(MAT,2) ).

TAU (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors. see Further Details.

The size of TAU must be min( size(MAT,1) , size(MAT,2) ).

USE_QR (INPUT, OPTIONAL) logical(lgl)

If the optional argument USE_QR is set to true, the input matrix is transposed, a fast QR decomposition is used for computing the LQ decomposition and the results are transposed again in the input matrix.

The default is true.

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(k) * … * H(2) * H(1), where k = min( size(MAT,1) , size(MAT,2) )

Each H(i) has the form

H(i) = I + tau * ( v * v’ ) ,

where tau is a real scalar and v is a real n-element vector with v(1:i-1) = 0. v(i:n) is stored on exit in MAT(i,i:n) and tau in TAU(i).

A blocked algorithm is used for computing the LQ factorization. Furthermore, the computations are parallelized if OPENMP is used.

On exit of LQ_CMP, the orthogonal matrix Q (or its first rows) can be computed explicitly by a call to subroutine ORTHO_GEN_LQ with arguments MAT and TAU.

For further details on the LQ factorization and its use or the blocked algorithm used here, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine ortho_gen_lq ( mat, tau, use_qr )

Purpose

ORTHO_GEN_LQ generates an m-by-n real matrix with orthonormal rows, which is defined as the first m rows of a product of k elementary reflectors of order n

Q = H(k) * … * H(2) * H(1)

as returned by LQ_CMP.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by LQ_CMP in the first k rows of its array argument MAT.

On exit, the first m rows of Q.

The shape of MAT must verify: size( MAT, 1 ) <= size( MAT, 2 ).

TAU (INPUT) real(stnd), dimension(:)

TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by LQ_CMP. The size of TAU determines the number k of elementary reflectors whose product defines the matrix Q.

The size of TAU must verify: size( TAU ) <= size( MAT, 1 ) .

USE_QR (INPUT, OPTIONAL) logical(lgl)

If the optional argument USE_QR is set to true, the input matrix is transposed, a fast QR type algorithm is used for computing the first m rows of Q and the results are transposed again in the input matrix.

The default is true.

Further Details

This subroutine used a blocked algorithm for agregating the Householder transformations (e.g. the elementary reflectors) stored in the upper triangle of MAT and generating the orthogonal matrix Q of the LQ factorization.

Furthermore, the computations are parallelized if OPENMP is used.

For further details on the LQ factorization and its use or the blocked algorithm, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine apply_q_lq ( mat, tau, c, left, trans )

Purpose

APPLY_Q_LQ overwrites the general real m-by-n matrix C with

Q * C if LEFT = true and TRANS = false, or

Q’ * C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ if LEFT = false and TRANS = true,

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

Q = H(k) * … * H(2) * H(1)

as returned by LQ_CMP. Q is of order m if LEFT = true and of order n if LEFT = false.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)
On entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by LQ_CMP in the first k rows of its array argument MAT. MAT is not modified by the routine.
TAU (INPUT) real(stnd), dimension(:)

TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by LQ_CMP. The size of TAU determines the number k of elementary reflectors whose product defines the matrix Q.

The size of TAU must verify: size( TAU ) <= min( size(MAT,1) , size(MAT,2) ) .

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

  • if LEFT = true, size( C, 1 ) = size( MAT, 2 )
  • if LEFT = false, size( C, 2 ) = size( MAT, 2 ) .
LEFT (INPUT) logical(lgl)

If:

  • LEFT = true : apply Q or Q’ from the left
  • LEFT = false : apply Q or Q’ from the right
TRANS (INPUT) logical(lgl)

If:

  • TRANS = false : apply Q (no transpose)
  • TRANS = true : apply Q’ (transpose)

Further Details

This subroutine is adapted from the routine DORML2 in LAPACK.

This subroutine used a blocked algorithm for agregating the Householder transformations (e.g. the elementary reflectors) stored in the upper triangle of MAT and applying the orthogonal matrix Q of the LQ factorization to the real m-by-n matrix C.

Furthermore, the computations are parallelized if OPENMP is used.

For further details on the LQ factorization and its use or the blocked algorithm used here, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine apply_q_lq ( mat, tau, c, trans )

Purpose

APPLY_Q_LQ overwrites the general real m vector C with

Q * C if TRANS = false, or

Q’ * C if TRANS = true,

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

Q = H(k) * … * H(2) * H(1)

as returned by LQ_CMP. Q is of order m.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)
On entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by LQ_CMP in the first k rows of its array argument MAT. MAT is not modified by the routine.
TAU (INPUT) real(stnd), dimension(:)

TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by LQ_CMP. The size of TAU determines the number k of elementary reflectors whose product defines the matrix Q.

The size of TAU must verify: size( TAU ) <= min( size(MAT,1) , size(MAT,2) ) .

C (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m vector C.

On exit, C is overwritten by Q * C or Q’ * C.

The shape of C must verify: size( C ) = size( MAT, 2 ).

TRANS (INPUT) logical(lgl)

If:

  • TRANS = false : apply Q (no transpose)
  • TRANS = true : apply Q’ (transpose)

Further Details

This subroutine is adapted from the routine DORML2 in LAPACK.

For further details on the LQ factorization and its use or the algorithm used here, see

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

subroutine qr_cmp ( mat, diagr, beta )

Purpose

QR_CMP computes a QR factorization of a real m-by-n matrix MAT :

MAT = Q * R

where Q is orthogonal and R is upper trapezoidal (upper triangular if m>=n).

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real matrix to be decomposed.

On exit, the elements above the diagonal of the array contain the corresponding elements of R. The elements on and below the diagonal, with the array BETA, represent the orthogonal matrix Q as a product of elementary reflectors, see Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R.

The size of DIAGR must verify: size( DIAGR ) <= min( size(MAT,1) , size(MAT,2) )

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors.

The size of BETA must verify:

  • size( BETA ) = size( DIAGR ) <= min( size(MAT,1) , size(MAT,2) )

See Further Details.

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(k), where k = size( BETA ) = size( DIAGR )

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

A blocked algorithm is used for computing the QR factorization. Furthermore, the computations are parallelized if OPENMP is used.

On exit of QR_CMP, the orthogonal matrix Q (or its first columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

For further details on the QR factorization and its use or the blocked algorithm used here, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine qrfac ( name_proc, syst, kfix, krank, min_norm, diagr, beta, h, tol, ip )

Purpose

QRFAC is a low-level subroutine for computing a QR or complete orthogonal factorization of the array section SYST(:m,:n) where m=size(SYST,1) and n<=size(SYST,2).

The routine first computes a QR factorization with column pivoting of SYST(:m,:n):

SYST(:m,:n) * P = Q * R

, where P is an n-by-n permutation matrix, R is an upper triangular or trapezoidal (i.e., if n>m) matrix and Q is a m-by-m orthogonal matrix.

The orthogonal transformation Q is then applied to the array section SYST(:m,n+1:):

SYST(:m,n+1:) = Q * B

Then, the rank of SYST(:m,:n) is determined by finding the submatrix R11 of R which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL or such that abs(R11[j,j])>0 if TOL is absent. The order of R11, krank, is an estimate of the rank of SYST(:m,:n).

This leads implicitly to the following partition of Q:

[ Q1 Q2 ]

where Q1 is an m-by-krank orthonormal matrix and Q2 is m-by-(m-krank) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a krank-by-krank triangular matrix, R21 is zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a full (m-krank)-by-(n-krank) matrix.

Q1 and Q2 are, respectively, orthonormal bases of the range and null space of SYST(:m,:n).

If MIN_NORM=true, R22 is considered to be negligible and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

SYST(:m,:n) * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-m orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a m-by-n matrix, which has the form:

[ T11 T12 ]

[ T21 T22 ]

Here T21 (=R21) and T12 are all zero, T22 (=R22) is considered to be negligible and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the vector argument IP, krank is stored in the scalar argument KRANK and if:

  • MIN_NORM=false, QRFAC computes Q and submatrices R11 and R12. Submatrices R11, and R12 are stored in the array section SYST(:krank,:n) and the array argument DIAGR. Q is stored compactly in factored form in the array section SYST(:m,:n) (in its lower triangle) and the array argument BETA.
  • MIN_NORM=true, QRFAC computes Q, Z and submatrice T11. Submatrice T11 is stored in the array section SYST(:krank,:krank) and the array argument DIAGR. Q is stored compactly in factored form in the array section SYST(:m,:n) (in its lower triangle) and the array argument BETA. Z is stored compactly in factored form in the array sections SYST(:krank,krank+1:n) and H(:krank).

See Further Details for more information.

Arguments

NAME_PROC (INPUT) character(len=*)
Name of the subroutine calling QRFAC.
SYST (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed and, eventually, the the right hand side matrix of an associated least squares problem in the array section SYST(:m,n+1:).

On exit, SYST(:m,1:n) has been overwritten by details of its QR or (complete) orthogonal factorization and B is stored in SYST(:m,n+1:).

See Further Details.

KFIX (INPUT) integer(i4b)

On entry:

  • KFIX=k, implies that the first k columns of SYST(:m,:n) are to be forced into the basis. Pivoting is performed only on the last n-k columns of SYST(:m,:n) if k<min(m,n).
  • KFIX<=0 can be used when pivoting is desired on all columns of SYST(:m,:n).
  • If KFIX<min(m,n) then the optional argument IP must be present to store the permutation matrix P.
  • When KFIX>=min(m,n) is used, pivoting is not performed. This is appropriate when SYST(:m,:n) is known to be of full rank.
KRANK (INPUT/OUTPUT) integer(i4b)

On entry, KRANK=n . KRANK must verify: KRANK <= size( SYST, 2 ).

On exit, KRANK contains the effective rank of SYST(:m,:n), i.e., the order of the submatrix R11, krank. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of SYST(:m,:n).

MIN_NORM (INPUT) logical(lgl)

On entry:

  • MIN_NORM=true indicates that a complete orthogonal factorization of SYST(:m,:n) must be computed.
  • MIN_NORM=false indicates that only a QR factorization (eventually with column pivoting if KFIX<min(m,n)) of SYST(:m,:n) must be computed.
DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R if MIN_NORM=false or the diagonal elements of the matrix T11 if MIN_NORM=true.

The diagonal elements of R are stored in DIAGR(:min(m,n)) and the diagonal elements of T11 are stored in DIAGR(:krank) and krank is stored in the real argument KRANK on exit.

See Further Details.

The size of DIAGR must verify:

  • size( DIAGR ) >= min( m , n ) = min( size(SYST,1) , KRANK ).
BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must verify:

  • size( BETA ) >= min( m , n ) = min( size(SYST,1) , KRANK ).
H (OUTPUT) real(stnd), dimension(:)

On exit, if MIN_NORM=true, the scalars factors of the elementary reflectors defining Z in the complete orthogonal factorization of MAT are stored in H(:krank). On exit, krank is output in the real argument KRANK.

See Further Details.

The size of H must verify: size( H ) >= n = KRANK .

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

If TOL is present and is in [0,1[, then:

  • On entry, the calculations to determine the condition number of of SYST(:m,:n) in the 1-norm are performed. Then, TOL is used to determine the effective rank of SYST(:m,:n), which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of SYST(:m,:n), whose estimated condition number is less than 1/TOL. If TOL=0 is specified the numerical rank of SYST(:m,:n) is determined.
  • On exit, the reciprocal of the condition number is returned in TOL.

If TOL is not specified or is outside [0,1[ :

  • The calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are done to determine the rank of MAT. If TOL is present, it is not changed.
IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of SYST(:m,:n) * P was the k-th column of SYST(:m,:n).

IP must be present if KFIX < min( m , n ) = min( SYST,1), KRANK ).

The size of IP must verify: size( IP ) >= n = KRANK.

Further Details

QRFAC is called by the subroutines QR_CMP2, LLSQ_QR_SOLVE, LLSQ_QR_SOLVE2 and SOLVE_LLSQ in modules QR_Procedures and LLSQ_Procedures.

Since QRFAC is a low level subroutine, no checking of the correctness of the dimensions of the array arguments is performed inside of the subroutine and such checking must be done before calling QRFAC.

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(k), where k = min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in SYST(i:m,i) and beta in BETA(i).

On exit of QRFAC, the orthonormal matrix Q (or its first n columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

The matrix P is represented in the array IP (if present) as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

If MIN_NORM=false, a QR factorization with column pivoting of SYST(:m,:n) is computed and, on exit:

  • The elements above the diagonal of the array SYST(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
  • The elements of the diagonal of R11 are stored in the array DIAGR.
  • The submatrix R12 is stored in SYST(:krank,krank+1:n).
  • krank is stored in the real argument KRANK.

If MIN_NORM=true, a complete orthogonal factorization of SYST(:m,:n) is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R, is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if MIN_NORM=true:

  • The scalar tau defining T(k) is returned in the kth element of H and the vector u(k) in the kth row of SYST, such that the elements of z(k) are in SYST(k,krank+1:n). The other elements of u(k) are not stored.
  • The elements above the diagonal of the array section SYST(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
  • krank is stored in the real argument KRANK.

The computations are parallelized if OPENMP is used. QRFAC uses a blocked “BLAS3” algorithm to compute the QR factorization of the columns of SYST(:m,:n), which are forced to be included in the basis as specified with the KFIX argument. However, if column pivoting is requested, QRFAC uses a standard “BLAS2” algorithm without any blocking on the columns which must be pivoted and is thus not optimized for computing a full QR factorization with column pivoting of very large matrices.

QRFAC is an updated version of a subroutine with the same name provided in the reference (1) with improvements suggested in the reference (3).

For further details, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Drmac, Z., and Bujanovic, Z., 2006:
    On the failure of rank revealing QR factorization software - a case study LAPACK Working Note 176.

subroutine qr_cmp2 ( mat, diagr, beta, ip, krank, tol, tau )

Purpose

QR_CMP2 computes a (complete) orthogonal factorization of a real m-by-n matrix MAT. MAT may be rank-deficient. The routine first computes a QR factorization with column pivoting of MAT:

MAT * P = Q * R

, here P is an n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, krank, is the effective rank of MAT.

This leads implicitly to the following partition of Q:

[ Q1 Q2 ]

where Q1 is an m-by-krank orthonormal matrix and Q2 is m-by-(m-krank) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a krank-by-krank triangular matrix, R21 is zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a full (m-krank)-by-(n-krank) matrix.

Q1 and Q2 are, respectively, orthonormal bases of the range and null space of MAT.

In a final step, if TAU is present, R22 is considered to be negligible and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

MAT * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-m orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a m-by-n matrix and has the form:

[ T11 T12 ]

[ T21 T22 ]

Here T21 (=R21) and T12 are all zero, T22 (=R22) is considered to be negligible and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the vector argument IP and if:

  • TAU is absent, QR_CMP2 computes Q and submatrices R11, R12. Submatrices R11 and R12 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA.
  • TAU is present, QR_CMP2 computes Q, Z and submatrice T11. Submatrice T11 is stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA. Z is stored compactly in factored form in the array arguments MAT and TAU.

See Further Details for more information on how the partial QR or orthogonal decomposition is stored in MAT on exit.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its (complete) orthogonal factorization.

See Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R if TAU is absent or the diagonal elements of the matrix T11 if TAU is present.

The diagonal elements of R or T11 are stored in DIAGR(1:krank) and krank is stored in the real argument KRANK on exit.

See Further Details.

The size of DIAGR must be equal to min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must be equal to min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (INPUT/OUTPUT) integer(i4b)

On entry:

  • KRANK=k, implies that the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.
  • KRANK<=0 can be used when pivoting is desired on all columns of MAT.
  • When KRANK>=min(m,n) is used, pivoting is not performed at all. This is appropriate when MAT is known to be full rank.

On exit, KRANK contains the effective rank of MAT, i.e., the order of the submatrix R11, krank. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

If TOL is present and is in [0,1[, then :

  • On entry, the calculations to determine the condition number of MAT are performed. Then, TOL is used to determine the effective rank of MAT, which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of MAT, whose estimated condition number < 1/TOL. If TOL=0 is specified the numerical rank of MAT is determined.
  • On exit, the reciprocal of the condition number is returned in TOL.

If TOL is not specified or is outside [0,1[ :

  • The calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are done to determine the numerical rank of MAT. If TOL is present, it is not changed.
TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a complete orthogonal factorization of MAT is computed. Otherwise, a QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the complete orthogonal factorization of MAT.

See Further Details.

The size of TAU must be equal to min( size(MAT,1) , size(MAT,2) ).

Further Details

  1. If it is possible that MAT may not be of full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be determined by using KRANK=0 and TOL=relative precision of the elements in MAT. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

  2. The matrix Q is represented as a product of elementary reflectors

    Q = H(1) * H(2) * … * H(k), where k = min( size(MAT,1) , size(MAT,2) )

    Each H(i) has the form

    H(i) = I + beta * ( v * v’ ) ,

    where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

    On exit of QR_CMP2, the orthonormal matrix Q (or its first n columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

  3. The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

  4. If TAU is absent, a QR factorization with column pivoting of MAT is computed and, on exit:

  • The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
  • The elements of the diagonal of R11 are stored in the array DIAGR.
  • The submatrix R12 is stored in MAT(:krank,krank+1:n).
  • krank is stored in the real argument KRANK.

  1. If TAU is present, a complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R12), is given in the form

    [ I 0 ]

    [ 0 T(k) ]

    where

    T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

    tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

    The Z n-by-n orthogonal matrix is given by

    Z = Z(1) * Z(2) * … * Z(krank)

    On exit, the scalar tau is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). On exit:

  • The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
  • The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
  • krank is stored in the real argument KRANK.

The computations are parallelized if OPENMP is used. QR_CMP2 uses a blocked “BLAS3” algorithm to compute the QR factorization of the columns of MAT, which are forced to be included in the basis as specified with the KRANK argument on entry. However, if column pivoting is requested, QR_CMP2 uses a standard “BLAS2” algorithm without any blocking on the columns which must be pivoted and is thus not optimized for computing a full QR factorization with column pivoting of very large matrices.

For computing (partial) QR factorization with column pivoting of very large matrices, subroutines PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 in module Random are a better choice. These two subroutines are much faster than QR_CMP2 on large matrices because of the use of randomized and blocked “BLAS3” algorithms instead of a standard “BLAS2” algorithm as in QR_CMP2. See references (4), (5) and (6) for further information on randomized QR algorithms.

For further details, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Drmac, Z., and Bujanovic, Z., 2006:
    On the failure of rank revealing QR factorization software - a case study LAPACK Working Note 176.
  4. Duersch, J.A., and Gu, M., 2017:
    Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.
  5. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:
    Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.
  6. Duersch, J.A., and Gu, M., 2020:
    Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine partial_qr_cmp ( mat, diagr, beta, ip, krank, tol, tau )

Purpose

PARTIAL_QR_CMP computes a (partial) QR factorization with column pivoting or complete orthogonal factorization of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (i.e., if n>m) matrix and Q is a m-by-m orthogonal matrix.

At the user option, the QR factorization can be only partial, e.g., the subroutine ends when the numbers of selected columns of Q is equal to a predefined value equals to kpartial = size( DIAGR ) = size( BETA ).

This leads implicitly to the following partition of Q:

[ Q1 Q2 ]

where Q1 is a m-by-kpartial orthonormal matrix and Q2 is a m-by-(m-kpartial) orthonormal matrix orthogonal to Q1, and to the following corresponding partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R11 is a kpartial-by-kpartial triangular matrix, R21 is zero by construction, R12 is a full kpartial-by-(n-kpartial) matrix and R22 is a full (m-kpartial)-by-(n-kpartial) matrix.

Then, if the optional scalar argument TOL is present and:

  • is in ]0,1[, the rank of R11 is determined by finding the submatrix of R11 which is defined as the largest leading submatrix whose estimated condition number, in the 1-norm, is less than 1/TOL. The order of this submatrix, krank, is the effective rank of R11 (and MAT if krank is less than kpartial or if krank=kpartial=min(m,n)).
  • is equal to 0, the rank of R11, krank, is determined by finding the largest submatrix of R11 such that abs(R11[j,j])>0.

In both cases, the order of this submatrix, krank, is the effective rank of R11 (and MAT if krank is less than kpartial or if krank=kpartial=min(m,n)).

If TOL is absent or outside [0,1[, the rank of R11 is not checked and is assumed to be equal to kpartial.

If krank is less than kpartial, then MAT is not of full rank (i.e., certain columns of MAT(:m,:kpartial) are linear combinations of other columns of MAT(:m,:kpartial)) and krank is also an estimate of the rank of MAT.

This leads to a redefinition of the partition of Q = [ Q1 Q2 ], where Q1 and Q2 are now m-by-krank and m-by-(m-krank) orthonormal matrices, and a corresponding redefinition of the associated partition of R, where R11 is now a krank-by-krank triangular matrix, R21 is again zero by construction, R12 is a full krank-by-(n-krank) matrix and R22 is a (m-krank)-by-(n-krank) matrix.

In a final step, if TAU is present, R22 is considered to be negligible and R12 is annihilated by orthogonal transformations from the right, arriving at the partial or complete orthogonal factorization:

MAT * P = Q * T * Z

, where P is a n-by-n permutation matrix, Q is a m-by-m orthogonal matrix, Z is a n-by-n orthogonal matrix and T is a m-by-n matrix and has the form:

[ T11 T12 ]

[ T21 T22 ]

Here T21 (=R21) and T12 are all zero, T22 (=R22) is considered to be negligible and T11 is a krank-by-krank upper triangular matrix.

On exit, P is stored compactly in the vector argument IP, krank is stored in the scalar argument KRANK and if:

  • TAU is absent, PARTIAL_QR_CMP computes Q and submatrices R11, R12 and R22. Submatrices R11, R12 and R22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA.
  • TAU is present, PARTIAL_QR_CMP computes Q, Z and submatrices T11 and T22 (=R22). Submatrices T11 and T22 are stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA. Z is stored compactly in factored form in the array arguments MAT and TAU.

See Further Details for more information.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its (partial) QR factorization.

See Further Details.

DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R11 or T11.

See Further Details.

The size of DIAGR must verify:

  • size( DIAGR ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).
BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must verify:

  • size( BETA ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).
IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of R11, i.e., krank, which is the order of this submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT and is also the rank of MAT if krank is less than kpartial = size( BETA ).

If the computed pseudo-rank, krank, is less than kpartial = size( BETA ), BETA(krank+1:kpartial) and, eventually, TAU(krank+1:kpartial) are set to zero and MAT(krank+1:m,krank+1:n) (e.g., R22=T22) is updated on exit. In other words, the subroutine outputs a partial QR factorization of rank krank instead of rank kpartial.

In all cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, if TOL is present then:

  • if TOL is in ]0,1[, the calculations to determine the condition number of R11 are performed. Then, TOL is used to determine the effective pseudo-rank of R11, which is defined as the order of the largest leading triangular submatrix in the partial QR factorization with column pivoting of MAT, whose estimated condition number in the 1-norm is less than 1/TOL. On exit, the reciprocal of the condition number is returned in TOL.
  • if TOL=0 is specified, the calculations to determine the condition number of R11 are not performed and crude tests on R(j,j) are done to determine the numerical pseudo-rank of R11. On exit, TOL is not changed.

If TOL is not specified or is outside [0,1[, the calculations to determine the rank of R11 are not performed and this rank is assumed to be equal to kpartial.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a (partial) complete orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the (partial) complete orthogonal factorization of MAT.

See Further Details.

The size of TAU must verify:

  • size( TAU ) = kpartial <= min( size(MAT,1) , size(MAT,2) ).

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank = size( BETA ) <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_QR_CMP, the orthonormal matrix Q (or its first n columns) can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA.

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

  • The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R11.
  • The elements of the diagonal of R11 are stored in the array DIAGR.
  • The submatrix R12 is stored in MAT(:krank,krank+1:n).
  • The submatrix R22 is stored in MAT(krank+1:m,krank+1:n).
  • krank is stored in the real argument KRANK.

If TAU is present, a partial or complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R12), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

  • The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
  • The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array DIAGR.
  • The submatrix T22 (=R22) is stored in MAT(krank+1:m,krank+1:n).
  • krank is stored in the real argument KRANK.

In both cases, norm(MAT(krank+1:m,krank+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

If it is possible that MAT may not be of full rank (i.e., certain columns of MAT are linear combinations of other columns), then the eventual linearly dependent columns in the partial QR decomposition of MAT, which is sought, can be determined by using TOL=relative precision of the elements in MAT and the partial QR or complete orthogonal factorization of MAT is adjusted accordingly. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

The computations are parallelized if OPENMP is used. However, note that PARTIAL_QR_CMP uses a standard “BLAS2” algorithm without any blocking and is thus not optimized for computing a partial QR or complete orthogonal factorization with column pivoting of very large matrices. For large matrices, subroutines PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 in module Random, which use randomized blocked “BLAS3” algorithms described in the references (4), (5) and (6), are a much better choice.

For further details, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Drmac, Z., and Bujanovic, Z., 2006:
    On the failure of rank revealing QR factorization software - a case study LAPACK Working Note 176.
  4. Duersch, J.A., and Gu, M., 2017:
    Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.
  5. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:
    Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.
  6. Duersch, J.A., and Gu, M., 2020:
    Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine partial_qr_cmp_fixed_precision ( mat, relerr, diagr, beta, ip, krank, tau )

Purpose

PARTIAL_QR_CMP_FIXED_PRECISION computes a partial QR factorization with column pivoting (or orthogonal factorization) of a m-by-n matrix MAT:

MAT * P = Q * R

, where P is a n-by-n permutation matrix, R is a krank-by-n (upper trapezoidal) matrix and Q is a m-by-krank matrix with orthogonal columns. This leads to the following matrix approximation of MAT of rank krank:

MAT = Q * (R * P’)

krank is the target rank of the matrix approximation, which is sought, and this partial factorization must have an approximation error which fulfills:

|| MAT - Q * ( R * P’ ) ||_F <= ||MAT||_F * relerr

|| ||_F is the Frobenius norm and relerr is a prescribed accuracy tolerance for the relative error of the computed matrix approximation, specified in the input argument RELERR.

PARTIAL_QR_CMP_FIXED_PRECISION searches incrementally the best (e.g., of smallest rank) Q * (R * P’) approximation, which fulfills the prescribed accuracy tolerance for the relative error. More precisely, the rank of the matrix approximation is increased progressively until the prescribed accuracy tolerance is satisfied.

In other words, the rank, krank, of the matrix approximation is not known in advance and is determined in the subroutine. krank is stored in the argument KRANK and the relative error of the computed partial matrix approximation in the argument RELERR on exit.

The computed partial matrix approximation leads implicitly to the following partition of R:

[ R1 R2 ]

where R1 is a krank-by-krank triangular matrix and R2 is a full krank-by-(n-krank) matrix.

In a final step, if TAU is present, R2 is annihilated by orthogonal transformations from the right, arriving at the partial orthogonal factorization:

MAT * P = Q * T1 * Z

, where P is a n-by-n permutation matrix, Q is a m-by-krank matrix with orthonormal columns Z is a krank-by-n matrix with orthonormal rows and T1 is a krank-by-krank upper triangular matrix.

Note, however, that this final step does not change the matrix approximation and its relative error, only the output format of this matrix approximation, which is now composed of four factors instead of three.

On exit, P is stored compactly in the vector argument IP, krank is stored in the scalar argument KRANK and if:

  • TAU is absent, PARTIAL_QR_CMP_FIXED_PRECISION computes Q and submatrices R1 and R2. Submatrices R1 and R2 are stored in the array arguments MAT and DIAGR(:KRANK). Q is stored compactly in factored form in the array arguments MAT and BETA(:KRANK).
  • TAU is present, PARTIAL_QR_CMP_FIXED_PRECISION computes Q, Z and submatrice T1. Submatrice T1 is stored in the array arguments MAT and DIAGR. Q is stored compactly in factored form in the array arguments MAT and BETA(:KRANK). Z is stored compactly in factored form in the array arguments MAT and TAU(:KRANK).

In all cases, the relative error of the computed matrix approximation is output in argument RELERR.

See Further Details for more information.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real m-by-n matrix to be decomposed.

On exit, MAT has been overwritten by details of its partial QR or complete orthogonal factorization.

See Further Details.

RELERR (INPUT/OUTPUT) real(stnd)

On entry, the requested accuracy tolerance for the relative error of the computed partial matrix approximation.

The preset value for RELERR must be greater than 4*epsilon( RELERR ) and less than one.

On exit, RELERR contains the relative error of the computed partial matrix approximation:

  • RELERR = || MAT - Q * ( R * P’ ) ||_F / ||MAT||_F
DIAGR (OUTPUT) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R1 or T1 are stored in the array section DIAGR(:KRANK). Other elements of DIAGR are set to zero on exit.

See Further Details.

The size of DIAGR must verify:

  • size( DIAGR ) = min( size(MAT,1) , size(MAT,2) ).
BETA (OUTPUT) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q are stored in the array section BETA(:KRANK). Other elements of BETA are set to zero on exit.

See Further Details.

The size of BETA must verify:

  • size( BETA ) = min( size(MAT,1) , size(MAT,2) ).
IP (OUTPUT) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be equal to size( MAT, 2 ) = n.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the rank of R1, i.e., krank, which is the order of this submatrix R1. This is the same as the order of the submatrix T1 in the “partial” complete orthogonal factorization of MAT and is also the rank of the computed matrix approximation.

In all cases, norm(MAT(KRANK+1:m,KRANK+1:n)) gives the error of the associated matrix approximation in the Frobenius norm, on exit.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a partial complete orthogonal factorization of MAT is computed. Otherwise, a simple QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining the orthogonal matrix Z in the partial complete orthogonal factorization of MAT are stored in the array section TAU(:KRANK). Other elements of TAU are set to zero on exit.

See Further Details.

The size of TAU must verify:

  • size( TAU ) = min( size(MAT,1) , size(MAT,2) ).

Further Details

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(krank), where krank <= min( m , n ).

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

On exit of PARTIAL_QR_CMP_FIXED_PRECISION, the matrix Q can be computed explicitly by a call to subroutine ORTHO_GEN_QR with arguments MAT and BETA(:KRANK).

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

On exit, if the optional argument TAU is absent:

  • The elements above the diagonal of the array MAT(:krank,:krank) contain the corresponding elements of the triangular matrix R1.
  • The elements of the diagonal of R1 are stored in the array DIAGR.
  • The submatrix R2 is stored in MAT(:krank,krank+1:n).
  • krank is stored in the real argument KRANK.

If TAU is present, a “partial” complete orthogonal factorization of MAT is computed. The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R (e.g., in R2), is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-krank) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R2.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(krank)

On exit, if the optional argument TAU is present:

  • The scalar tau defining T(k) is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,krank+1:n). The other elements of u(k) are not stored.
  • The elements above the diagonal of the array section MAT(:krank,:krank) contain the corresponding elements of the triangular matrix T1. The elements of the diagonal of T1 are stored in the array DIAGR.
  • krank is stored in the real argument KRANK.

In both cases, norm(MAT(KRANK+1:m,KRANK+1:n)) gives the error of the associated partial matrix approximation in the Frobenius norm, and argument RELERR stores the relative error in the Frobenius norm of the matrix approximation on exit.

The computations are parallelized if OPENMP is used. However, note that PARTIAL_QR_CMP_FIXED_PRECISION uses a standard “BLAS2” algorithm without any blocking and is thus not optimized for computing a partial QR or complete orthogonal factorization with column pivoting of very large matrices. For large matrices, subroutine PARTIAL_RQR_CMP_FIXED_PRECISION in module Random, which uses a randomized blocked “BLAS3” algorithm described in the references (4), (5) and (6) is a much better choice.

For further details, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Drmac, Z., and Bujanovic, Z., 2006:
    On the failure of rank revealing QR factorization software - a case study LAPACK Working Note 176.
  4. Duersch, J.A., and Gu, M., 2017:
    Randomized QR with column pivoting. SIAM J. Sci. Comput., Volume 39, Issue 4, C263-C291.
  5. Martinsson, P.G., Quintana-Orti, G., Heavner, N., and Van de Geijn, R., 2017:
    Householder QR factorization with randomization for column pivoting (HQRRP). SIAM J. Sci. Comput., Volume 39, Issue 2, C96-C115.
  6. Duersch, J.A., and Gu, M., 2020:
    Randomized projection for rank-revealing matrix factorizations and low-rank approximations. SIAM Review, Volume 62, Issue 3, 661-682.

subroutine ortho_gen_qr ( mat, beta )

Purpose

ORTHO_GEN_QR generates an m-by-n real matrix 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 QR_CMP or QR_CMP2.

Arguments

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by QR_CMP (or QR_CMP2) in the first k columns of its array argument MAT.

On exit, the first n columns of Q.

The shape of MAT must verify: size( MAT, 2 ) <= size( MAT, 1 ).

BETA (INPUT) real(stnd), dimension(:)

BETA(i) must contain the scalar factor of the elementary reflector H(i), as returned by QR_CMP (or QR_CMP2). The size of BETA determines the number k of elementary reflectors whose product defines the matrix Q.

The size of BETA must verify: size( BETA ) <= ( MAT, 2 ) .

Further Details

This subroutine used a blocked algorithm for agregating the Householder transformations (e.g. the elementary reflectors) stored in the lower triangle of MAT and generating the orthogonal matrix Q of the QR factorization.

Furthermore, the computations are parallelized if OPENMP is used.

For further details on the QR factorization and its use or the blocked algorithm, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine apply_q_qr ( mat, beta, c, left, trans )

Purpose

APPLY_Q_QR overwrites the general real m-by-n matrix C with

Q * C if LEFT = true and TRANS = false, or

Q’ * C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ if LEFT = false and TRANS = true,

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

Q = H(1) * H(2) * … * H(k)

as returned by QR_CMP or QR_CMP2. Q is of order m if LEFT = true and of order n if LEFT = false.

Arguments

MAT (INPUT) real(stnd), dimension(:,:)
On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by QR_CMP (or QR_CMP2) in the first k columns of its array argument MAT. MAT is not modified by the routine.
BETA (INPUT) real(stnd), dimension(:)

BETA(i) must contain the scalar factor of the elementary reflector H(i), as returned by QR_CMP (or QR_CMP2). The size of BETA determines the number k of elementary reflectors whose product defines the matrix Q.

The size of BETA must verify:

  • size( BETA ) <= min( size(MAT,1) , size(MAT,2) ) .
C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

  • if LEFT = true, size( C, 1 ) = size( MAT, 1 )
  • if LEFT = false, size( C, 2 ) = size( MAT, 1 ) .
LEFT (INPUT) logical(lgl)

If:

  • LEFT = true : apply Q or Q’ from the left
  • LEFT = false : apply Q or Q’ from the right
TRANS (INPUT) logical(lgl)

If:

  • TRANS = false : apply Q (no transpose)
  • TRANS = true : apply Q’ (transpose)

Further Details

This subroutine is adapted from the routine DORM2R in LAPACK.

This subroutine used a blocked algorithm for agregating the Householder transformations (e.g. the elementary reflectors) stored in the lower triangle of MAT and applying the orthogonal matrix Q of the QR factorization to the real m-by-n matrix C.

Furthermore, the computations are parallelized if OPENMP is used.

For further details on the QR factorization and its use or the blocked algorithm used here, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
  3. Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:
    Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.
  4. Walker, H.F., 1988:
    Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

subroutine apply_q_qr ( mat, beta, c, trans )

Purpose

APPLY_Q_QR overwrites the real m vector C with

Q * C if TRANS = false, or

Q’ * C if TRANS = true

where Q is a real orthogonal matrix defined as the product of k elementary reflectors

Q = H(1) * H(2) * … * H(k)

as returned by QR_CMP or QR_CMP2. Q is of order m .

Arguments

MAT (INPUT) real(stnd), dimension(:,:)
On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,…,k, as returned by QR_CMP (or QR_CMP2) in the first k columns of its array argument MAT. MAT is not modified by the routine.
BETA (INPUT) real(stnd), dimension(:)

BETA(i) must contain the scalar factor of the elementary reflector H(i), as returned by QR_CMP (or QR_CMP2). The size of BETA determines the number k of elementary reflectors whose product defines the matrix Q.

The size of BETA must verify: size( BETA ) <= min( size(MAT,1) , size(MAT,2) ) .

C (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m vector C.

On exit, C is overwritten by Q * C or Q’ * C .

The shape of C must verify: size( C ) = size( MAT, 1 ) .

TRANS (INPUT) logical(lgl)

If:

  • TRANS = false : apply Q (no transpose)
  • TRANS = true : apply Q’ (transpose)

Further Details

This subroutine is adapted from the routine DORM2R in LAPACK.

For further details on the QR factorization and its use or the algorithm used here, see:

  1. Lawson, C.L., and Hanson, R.J., 1974:
    Solving least square problems. Prentice-Hall.
  2. Golub, G.H., and Van Loan, C.F., 1996:
    Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.
Flag Counter