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.

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

NCA — Number of columns of A. (Input)
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).

PIVOT — Logical variable. (Input)
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.

LDQR — Leading dimension of QR exactly as specified in the dimension statement of the calling program. (Input)
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.

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.

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.

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 functions (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