LQRRR
Computes the QR decomposition, AP = QR, using Householder transformations.
Required Arguments
A — Real NRA by NCA matrix containing the matrix whose QR factorization is to be computed. (Input)
QR — Real NRA by NCA matrix containing information required for the QR factorization. (Output)
The upper trapezoidal part of QR contains the upper trapezoidal part of R with its diagonal elements ordered in decreasing magnitude. The strict lower trapezoidal part of QR contains information to recover the orthogonal matrix Q of the factorization. Arguments A and QR can occupy the same storage locations. In this case, A will not be preserved on output.
The upper trapezoidal part of QR contains the upper trapezoidal part of R with its diagonal elements ordered in decreasing magnitude. The strict lower trapezoidal part of QR contains information to recover the orthogonal matrix Q of the factorization. Arguments A and QR can occupy the same storage locations. In this case, A will not be preserved on output.
QRAUX — Real vector of length NCA containing information about the orthogonal part of the decomposition in the first min(NRA, NCA) position. (Output)
Optional Arguments
NRA — Number of rows of A. (Input)
Default: NRA = size (A,1).
Default: NRA = size (A,1).
NCA — Number of columns of A. (Input)
Default: NCA = size (A,2).
Default: NCA = size (A,2).
LDA — Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)
Default: LDA = size (A,1).
Default: LDA = size (A,1).
PIVOT — Logical variable. (Input)
PIVOT = .TRUE. means column pivoting is enforced.
PIVOT = .FALSE. means column pivoting is not done.
Default: PIVOT = .TRUE.
PIVOT = .TRUE. means column pivoting is enforced.
PIVOT = .FALSE. means column pivoting is not done.
Default: PIVOT = .TRUE.
IPVT — Integer vector of length NCA containing information that controls the final order of the columns of the factored matrix A. (Input/Output)
On input, if IPVT(K) > 0, then the K-th column of A is an initial column. If IPVT(K) = 0, then the K-th column of A is a free column. If IPVT(K) < 0, then the K-th column of A is a final column. See the Comments section below. On output, IPVT(K) contains the index of the column of A that has been interchanged into the K-th column. This defines the permutation matrix P. The array IPVT is referenced only if PIVOT is equal to .TRUE.
Default: IPVT = 0.
On input, if IPVT(K) > 0, then the K-th column of A is an initial column. If IPVT(K) = 0, then the K-th column of A is a free column. If IPVT(K) < 0, then the K-th column of A is a final column. See the Comments section below. On output, IPVT(K) contains the index of the column of A that has been interchanged into the K-th column. This defines the permutation matrix P. The array IPVT is referenced only if PIVOT is equal to .TRUE.
Default: IPVT = 0.
LDQR — Leading dimension of QR exactly as specified in the dimension statement of the calling program. (Input)
Default: LDQR = size (QR,1).
Default: LDQR = size (QR,1).
CONORM — Real vector of length NCA containing the norms of the columns of the input matrix. (Output)
If this information is not needed, CONORM and QRAUX can share the same storage locations.
If this information is not needed, CONORM and QRAUX can share the same storage locations.
FORTRAN 90 Interface
Generic: CALL LQRRR (A, QR, QRAUX [, …])
Specific: The specific interface names are S_LQRRR and D_LQRRR.
FORTRAN 77 Interface
Single: CALL LQRRR (NRA, NCA, A, LDA, PIVOT, IPVT, QR, LDQR, QRAUX, CONORM)
Double: The double precision name is DLQRRR.
ScaLAPACK Interface
Generic: CALL LQRRR (A0, QR0, QRAUX0 [, …])
Specific: The specific interface names are S_LQRRR and D_LQRRR.
See the ScaLAPACK Usage Notes below for a description of the arguments for distributed computing.
Description
The routine LQRRR computes the QR decomposition of a matrix using Householder transformations. The underlying code is based on either LINPACK , LAPACK, or ScaLAPACK code depending upon which supporting libraries are used during linking. For a detailed explanation see Using ScaLAPACK, LAPACK, LINPACK, and EISPACK in the Introduction section of this manual.
LQRRR determines an orthogonal matrix Q, a permutation matrix P, and an upper trapezoidal matrix R with diagonal elements of nonincreasing magnitude, such that AP = QR. The Householder transformation for column k is of the form
for k = 1, 2, …, min(NRA, NCA), where u has zeros in the first k − 1 positions. The matrix Q is not produced directly by LQRRR . Instead the information needed to reconstruct the Householder transformations is saved. If the matrix Q is needed explicitly, the subroutine LQERR can be called after LQRRR. This routine accumulates Q from its factored form.
Before the decomposition is computed, initial columns are moved to the beginning of the array A and the final columns to the end. Both initial and final columns are frozen in place during the computation. Only free columns are pivoted. Pivoting, when requested, is done on the free columns of largest reduced norm.
Comments
1. Workspace may be explicitly provided, if desired, by use of L2RRR/DL2RRR. The reference is:
CALL L2RRR (NRA, NCA, A, LDA, PIVOT, IPVT, QR, LDQR, QRAUX, CONORM, WORK)
The additional argument is
WORK — Work vector of length 2NCA − 1. Only NCA − 1 locations of WORK are referenced if PIVOT = .FALSE. .
2. LQRRR determines an orthogonal matrix Q, permutation matrix P, and an upper trapezoidal matrix R with diagonal elements of nonincreasing magnitude, such that AP = QR. The Householder transformation for column k, k = 1, …, min(NRA, NCA) is of the form
where u has zeros in the first k − 1 positions. If the explicit matrix Q is needed, the user can call routine LQERR after calling LQRRR. This routine accumulates Q from its factored form.
3. Before the decomposition is computed, initial columns are moved to the beginning and the final columns to the end of the array A. Both initial and final columns are not moved during the computation. Only free columns are moved. Pivoting, if requested, is done on the free columns of largest reduced norm.
4. When pivoting has been selected by having entries of IPVT initialized to zero, an estimate of the condition number of A can be obtained from the output by computing the magnitude of the number QR(1, 1)/QR(K, K), where K = MIN(NRA, NCA). This estimate can be used to select the number of columns, KBASIS, used in the solution step computed with routine LQRSL.
ScaLAPACK Usage Notes
The arguments which differ from the standard version of this routine are:
A0 — MXLDA by MXCOL local matrix containing the local portions of the distributed matrix A. A contains the matrix whose QR factorization is to be computed. (Input)
QR0 — MXLDA by MXCOL local matrix containing the local portions of the distributed matrix QR. QR contains the information required for the QR factorization. (Output)
The upper trapezoidal part of QR contains the upper trapezoidal part of R with its diagonal elements ordered in decreasing magnitude. The strict lower trapezoidal part of QR contains information to recover the orthogonal matrix Q of the factorization. Arguments A and QR can occupy the same storage locations. In this case, A will not be preserved on output.
The upper trapezoidal part of QR contains the upper trapezoidal part of R with its diagonal elements ordered in decreasing magnitude. The strict lower trapezoidal part of QR contains information to recover the orthogonal matrix Q of the factorization. Arguments A and QR can occupy the same storage locations. In this case, A will not be preserved on output.
QRAUX0 — Real vector of length MXCOL containing the local portions of the distributed matrix QRAUX. QRAUX contains information about the orthogonal part of the decomposition in the first MIN(NRA, NCA) position. (Output)
IPVT0 — Integer vector of length MXLDB containing the local portions of the distributed vector IPVT. IPVT contains the information that controls the final order of the columns of the factored matrix A. (Input/Output)
On input, if IPVT(K) > 0, then the K-th column of A is an initial column. If IPVT(K) = 0, then the K-th column of A is a free column. If IPVT(K) < 0, then the K-th column of A is a final column. See Comments.
On output, IPVT(K) contains the index of the column of A that has been interchanged into the K‑th column. This defines the permutation matrix P. The array IPVT is referenced only if PIVOT is equal to .TRUE.
Default: IPVT = 0.
On input, if IPVT(K) > 0, then the K-th column of A is an initial column. If IPVT(K) = 0, then the K-th column of A is a free column. If IPVT(K) < 0, then the K-th column of A is a final column. See Comments.
On output, IPVT(K) contains the index of the column of A that has been interchanged into the K‑th column. This defines the permutation matrix P. The array IPVT is referenced only if PIVOT is equal to .TRUE.
Default: IPVT = 0.
All other arguments are global and are the same as described for the standard version of the routine. In the argument descriptions above, MXLDA, MXLDB, and MXCOL can be obtained through a call to SCALAPACK_GETDIM (see Utilities) after a call to SCALAPACK_SETUP (see Utilities) has been made. See the ScaLAPACK Example below.
Examples
Example
In various statistical algorithms it is necessary to compute q = xT(AT A)-1x, where A is a rectangular matrix of full column rank. By using the QR decomposition, q can be computed without forming ATA. Note that
AT A = (QRP-1)T(QRP-1) = P-T RT (QTQ)RP-1 = P RTRPT
since Q is orthogonal (QTQ = I) and P is a permutation matrix. Let
where R1 is an upper triangular nonsingular matrix. Then
In the following program, first the vector t = P-1 x is computed. Then
t := R1-Tt
Finally,
q = ∥t∥2
USE IMSL_LIBRARIES
! Declare variables
INTEGER LDA, LDQR, NCA, NRA
PARAMETER (NCA=3, NRA=4, LDA=NRA, LDQR=NRA)
! SPECIFICATIONS FOR PARAMETERS
INTEGER LDQ
PARAMETER (LDQ=NRA)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IPVT(NCA), NOUT
REAL CONORM(NCA), Q, QR(LDQR,NCA), QRAUX(NCA), T(NCA)
LOGICAL PIVOT
REAL A(LDA,NCA), X(NCA)
!
! Set values for A
!
! A = ( 1 2 4 )
! ( 1 4 16 )
! ( 1 6 36 )
! ( 1 8 64 )
!
DATA A/4*1.0, 2.0, 4.0, 6.0, 8.0, 4.0, 16.0, 36.0, 64.0/
!
! Set values for X
!
! X = ( 1 2 3 )
DATA X/1.0, 2.0, 3.0/
!
! QR factorization
PIVOT = .TRUE.
IPVT=0
CALL LQRRR (A, QR, QRAUX, pivot=pivot, IPVT=IPVT)
! Set t = inv(P)*x
CALL PERMU (X, IPVT, T, IPATH=1)
! Compute t = inv(trans(R))*t
CALL LSLRT (QR, T, T, IPATH=4)
! Compute 2-norm of t, squared.
Q = SDOT(NCA,T,1,T,1)
! Print result
CALL UMACH (2, NOUT)
WRITE (NOUT,*) ’Q = ’, Q
!
END
Output
Q = 0.840624
ScaLAPACK Example
The previous example is repeated here as a distributed computing example. In various statistical algorithms it is necessary to compute q = xT(AT A)-1x, where A is a rectangular matrix of full column rank. By using the QR decomposition, q can be computed without forming AT A. Note that
ATA = (QRP-1)T(QRP-1) = P-TRT(QTQ)RP-1 = P RTRPT
since Q is orthogonal (QTQ = I) and P is a permutation matrix. Let
where R1 is an upper triangular nonsingular matrix. Then
In the following program, first the vector t = P-1 x is computed. Then
t := R1-Tt
Finally,
q = ∥t∥2
SCALAPACK_MAP and SCALAPACK_UNMAP are IMSL utility routines (see Utilities) used to map and unmap arrays to and from the processor grid. They are used here for brevity. DESCINIT is a ScaLAPACK tools routine which initializes the descriptors for the local arrays.
USE MPI_SETUP_INT
USE LQRRR_INT
USE PERMU_INT
USE LSLRT_INT
USE UMACH_INT
USE SCALAPACK_SUPPORT
IMPLICIT NONE
INCLUDE ‘mpif.h’
! Declare variables
INTEGER LDA, LDQR, NCA, NRA, DESCA(9), DESCB(9), DESCL(9)
INTEGER INFO, MXCOL, MXLDA, MXLDB, MXCOLB, NOUT
INTEGER, ALLOCATABLE :: IPVT(:), IPVT0(:)
LOGICAL PIVOT
REAL Q
REAL, ALLOCATABLE :: A(:,:), X(:), T(:)
REAL, ALLOCATABLE :: A0(:,:), T0(:), QR0(:,:), QRAUX0(:)
REAL, (KIND(1E0))SDOT
PARAMETER (NRA=4, NCA=3, LDA=NRA, LDQR=NRA)
! Set up for MPI
MP_NPROCS = MP_SETUP()
IF(MP_RANK .EQ. 0) THEN
ALLOCATE (A(LDA,NCA), X(NCA), T(NCA), IPVT(NCA))
! Set values for A and the righthand side
A(1,:) = (/ 1.0, 2.0, 4.0/)
A(2,:) = (/ 1.0, 4.0, 16.0/)
A(3,:) = (/ 1.0, 6.0, 36.0/)
A(4,:) = (/ 1.0, 8.0, 64.0/)
!
X = (/ 1.0, 2.0, 3.0/)
!
IPVT = 0
ENDIF
! Set up a 1D processor grid and define
! its context ID, MP_ICTXT
CALL SCALAPACK_SETUP(NRA, NCA, .TRUE., .TRUE.)
! Get the array descriptor entities MXLDA,
! MXCOL, MXLDB, MXCOLB
CALL SCALAPACK_GETDIM(NRA, NCA, MP_MB, MP_NB, MXLDA, MXCOL)
CALL SCALAPACK_GETDIM(NCA, 1, MP_NB, 1, MXLDB, MXCOLB)
! Set up the array descriptors
CALL DESCINIT(DESCA, NRA, NCA, MP_MB, MP_NB, 0, 0, MP_ICTXT, MXLDA, &
INFO)
CALL DESCINIT(DESCL, 1, NCA, 1, MP_NB, 0, 0, MP_ICTXT, 1, INFO)
CALL DESCINIT(DESCB, NCA, 1, MP_NB, 1, 0, 0, MP_ICTXT, MXLDB, &
INFO)
! Allocate space for the local arrays
ALLOCATE (A0(MXLDA,MXCOL), QR0(MXLDA,MXCOL), QRAUX0(MXCOL), &
IPVT0(MXCOL), T0(MXLDB))
! Map input array to the processor grid
CALL SCALAPACK_MAP(A, DESCA, A0)
PIVOT = .TRUE.
CALL SCALAPACK_MAP(IPVT, DESCL, IPVT0)
! QR factorization
CALL LQRRR (A0, QR0, QRAUX0, PIVOT=PIVOT, IPVT=IPVT0)
! Unmap the results from the distributed
! array back to a non-distributed array.
! After the unmap, only Rank=0 has the full
! array.
CALL SCALAPACK_UNMAP(IPVT0, DESCL, IPVT, NCA, .FALSE.)
IF(MP_RANK .EQ. 0) CALL PERMU (X, IPVT, T, IPATH=1)
CALL SCALAPACK_MAP(T, DESCB, T0)
CALL LSLRT (QR0, T0, T0, IPATH=4)
CALL SCALAPACK_UNMAP(T0, DESCB, T)
! Print results.
! Only Rank=0 has the solution.
IF(MP_RANK .EQ. 0)THEN
Q = SDOT(NCA, T, 1, T, 1)
CALL UMACH (2, NOUT)
WRITE (NOUT, *) ‘Q = ‘, Q
ENDIF
! Exit ScaLAPACK usage
CALL SCALAPACK_EXIT(MP_ICTXT)
! Shut down MPI
MP_NPROCS = MP_SETUP(‘FINAL’)
END
Output
Q = 0.840624
