Module_QR_Procedures

Copyright 2020 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.

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

LATEST REVISION : 26/11/2020

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 suboutine 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 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.

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 suboutine 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 suboutine 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 suboutine 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 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 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 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.