Module_QR_Procedures¶

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 AND LQ DECOMPOSITIONS.

LATEST REVISION : 31/10/2018

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

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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 ) <= ( 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 ), ( 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 ), ( 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

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.

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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 (complete) orthogonal factorization of the array section SYST(1:m,1:n) where n<=size(SYST,2) and m=size(SYST,1).

The routine first computes a QR factorization with column pivoting :

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

, where P is 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.

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

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

Then, the rank of SYST(1:m,1: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 the effective rank of SYST(1:m,1:n).

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal matrix) and R22 is considered to be negligible.

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(1:m,1: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 and has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

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(:,n+1:).

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

KFIX (INPUT) integer(i4b)

KFIX=k, implies that the first k columns of SYST(:,1:n) are to be forced into the basis. Pivoting is performed on the last n-k columns of SYST(:,1:n). KFIX<=0 can be used when pivoting is desired on all columns of SYST(:,1: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(:,1:n) is known to be 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(:,1:n), i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of SYST(:,1:n).

MIN_NORM (INPUT) logical(lgl)

MIN_NORM=true indicates that a complete orthogonal factorization of SYST(:,1:n) must be computed. MIN_NORM=false indicates that only a QR factorization (with column pivoting) of SYST(:,1: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 T11 are stored in DIAGR(1:KRANK). see Further Details.

The size of DIAGR must verify:

size( DIAGR ) >= 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( size(SYST,1) , KRANK ) .

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

On exit, if MIN_NORM=true, the scalars factors of the elementary reflectors defining Z are stored in H(1:KRANK). see Further Details.

The size of H must verify size(H) >= 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(:,1:n) are performed. Then, TOL is used to determine the effective rank of SYST(:,1:n), which is defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of SYST(:,1:n), whose estimated condition number < 1/TOL. If TOL=0 is specified the numerical rank of SYST(:,1: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(:,1:n) * P was the k-th column of SYST(:,1:n). IP must be present if KFIX<min( SYST,1), KRANK )=min( m, n).

The size of IP must verify size(IP) >= 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).

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(:,1:n) is computed and, on exit, the elements above the diagonal of the array SYST(:,1:n) contain the corresponding elements of the triangular matrix R. The elements of the diagonal of R are stored in the array DIAGR.

If MIN_NORM=true, a complete orthogonal factorization of SYST(:,1: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 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 SYST, such that the elements of z(k) are in SYST(k,KRANK+1:n).

On exit, the elements above the diagonal of the array section SYST(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK).

The computations are parallelized if OPENMP is used.

For further details, see

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

G.H. Golub and C.F. Van Loan, 1996: Matrix Computations. The Johns

Hopkins University Press, Baltimore, Maryland.

`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 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 to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal matrix).

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), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

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 T11 are stored in DIAGR(1:KRANK).

The size of DIAGR must be 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 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 size(MAT,2)=n.

KRANK (INPUT/OUTPUT) integer(i4b)

On input, 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. When KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank. KRANK<=0 can be used when pivoting is desired on all columns of MAT.

On exit, KRANK contains the effective rank of MAT, i.e., the order of the submatrix R11. 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, only a QR factorization with column pivoting of MAT is computed.

On exit, the scalars factors of the elementary reflectors defining Z. See Further Details.

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

Further Details¶

If it is possible that MAT may not be 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.

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

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.

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 contain the corresponding elements of the triangular matrix R. The elements of the diagonal of R are stored in the array DIAGR.

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 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 elements above the diagonal of the array section MAT(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK).

The computations are parallelized if OPENMP is used.

For further details, see

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

G.H. Golub and C.F. Van Loan, 1996: Matrix Computations. The Johns

Hopkins University Press, Baltimore, Maryland.

`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 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 ), ( 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

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.

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 ), ( 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

Lawson, C.L., and Hanson, R.J., 1974: Solving least square problems.

Prentice-Hall.

Golub, G.H., and van Loan, C.F., 1996: Matrix Computations, 3rd ed.

The Johns Hopkins University Press, Baltimore.

Module_QR_Procedures¶

subroutine lq_cmp ( mat, diagl, tau, use_qr )¶

Purpose¶

Arguments¶

Further Details¶

subroutine ortho_gen_lq ( mat, tau, use_qr )¶

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

subroutine qr_cmp ( mat, diagr, beta )¶

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

subroutine ortho_gen_qr ( mat, beta )¶

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

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

Purpose¶

Arguments¶

Further Details¶

`subroutine lq_cmp ( mat, diagl, tau, use_qr )`¶

`subroutine ortho_gen_lq ( mat, tau, use_qr )`¶

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

`subroutine apply_q_lq ( mat, tau, c, trans )`¶

`subroutine qr_cmp ( mat, diagr, beta )`¶

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

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

`subroutine ortho_gen_qr ( mat, beta )`¶

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

`subroutine apply_q_qr ( mat, beta, c, trans )`¶