LQERR
Accumulates the orthogonal matrix Q from its factored form given the QR factorization of a rectangular matrix A.
Required Arguments
QR — Real NRQR by NCQR matrix containing the factored form of the matrix Q in the first min(NRQR, NCQR) columns of the strict lower trapezoidal part of QR as output from subroutine LQRRR/DLQRRR. (Input)
QRAUX — Real vector of length NCQR containing information about the orthogonal part of the decomposition in the first min(NRQR, NCQR) position as output from routine LQRRR/DLQRRR. (Input)
Q — Real NRQR by NRQR matrix containing the accumulated orthogonal matrix Q; Q and QR can share the same storage locations if QR is not needed. (Output)
Optional Arguments
NRQR — Number of rows in QR. (Input)
Default: NRQR = size (QR,1).
Default: NRQR = size (QR,1).
NCQR — Number of columns in QR. (Input)
Default: NCQR = size (QR,2).
Default: NCQR = size (QR,2).
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).
LDQ — Leading dimension of Q exactly as specified in the dimension statement of the calling program. (Input)
Default: LDQ = size (Q,1).
Default: LDQ = size (Q,1).
FORTRAN 90 Interface
Generic: CALL LQERR (QR, QRAUX, Q [, …])
Specific: The specific interface names are S_LQERR and D_LQERR.
FORTRAN 77 Interface
Single: CALL LQERR (NRQR, NCQR, QR, LDQR, QRAUX, Q, LDQ)
Double: The double precision name is DLQERR.
ScaLAPACK Interface
Generic: CALL LQERR (QR0, QRAUX0, Q0 [, …])
Specific: The specific interface names are S_LQERR and D_LQERR.
See the ScaLAPACK Usage Notes below for a description of the arguments for distributed computing.
Description
The routine LQERR accumulates the Householder transformations computed by IMSL routine LQRRR to produce the orthogonal matrix Q.
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.
Comments
1. Workspace may be explicitly provided, if desired, by use of L2ERR/DL2ERR. The reference is:
CALL L2ERR (NRQR, NCQR, QR, LDQR, QRAUX, Q, LDQ, WK)
The additional argument is
WK — Work vector of length 2 * NRQR.
ScaLAPACK Usage Notes
The arguments which differ from the standard version of this routine are:
QR0 — MXLDA by MXCOL local matrix containing the local portions of the distributed matrix QR. QR contains the factored form of the matrix Q in the first min(NRQR, NCQR) columns of the strict lower trapezoidal part of QR as output from subroutine LQRRR/DLQRRR. (Input)
QRAUX0 — Real vector of length MXCOL containing the local portions of the distributed matrix QRAUX. QRAUX contains the information about the orthogonal part of the decomposition in the first min(NRA, NCA) positions as output from subroutine LQRRR/DLQRRR. (Input)
Q0 — MXLDA by MXLDA local matrix containing the local portions of the distributed matrix Q. Q contains the accumulated orthogonal matrix ; Q and QR can share the same storage locations if QR is not needed. (Output)
All other arguments are global and are the same as described for the standard version of the routine. In the argument descriptions above, MXLDA 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 1
In this example, the orthogonal matrix Q in the QR decomposition of a matrix A is computed. The product X = QR is also computed. Note that X can be obtained from A by reordering the columns of A according to IPVT.
USE IMSL_LIBRARIES
! Declare variables
INTEGER LDA, LDQ, LDQR, NCA, NRA
PARAMETER (NCA=3, NRA=4, LDA=NRA, LDQ=NRA, LDQR=NRA)
!
INTEGER IPVT(NCA), J
REAL A(LDA,NCA), CONORM(NCA), Q(LDQ,NRA), QR(LDQR,NCA), &
QRAUX(NCA), R(NRA,NCA), X(NRA,NCA)
LOGICAL PIVOT
!
! 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/
!
! QR factorization
! Set IPVT = 0 (all columns free)
IPVT = 0
PIVOT = .TRUE.
CALL LQRRR (A, QR, QRAUX, IPVT=IPVT, PIVOT=PIVOT)
! Accumulate Q
CALL LQERR (QR, QRAUX, Q)
! R is the upper trapezoidal part of QR
R = 0.0E0
DO 10 J=1, NCA
CALL SCOPY (J, QR(:,J), 1, R(:,J), 1)
10 CONTINUE
! Compute X = Q*R
CALL MRRRR (Q, R, X)
! Print results
CALL WRIRN (’IPVT’, IPVT, 1, NCA, 1)
CALL WRRRN (’Q’, Q)
CALL WRRRN (’R’, R)
CALL WRRRN (’X = Q*R’, X)
!
END
Output
IPVT
1 2 3
3 2 1
Q
1 2 3 4
1 -0.0531 -0.5422 0.8082 -0.2236
2 -0.2126 -0.6574 -0.2694 0.6708
3 -0.4783 -0.3458 -0.4490 -0.6708
4 -0.8504 0.3928 0.2694 0.2236
R
1 2 3
1 -75.26 -10.63 -1.59
2 0.00 -2.65 -1.15
3 0.00 0.00 0.36
4 0.00 0.00 0.00
X = Q*R
1 2 3
1 4.00 2.00 1.00
2 16.00 4.00 1.00
3 36.00 6.00 1.00
4 64.00 8.00 1.00
ScaLAPACK Example
In this example, the orthogonal matrix Q in the QR decomposition of a matrix A is computed. 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 LQERR_INT
USE WRRRN_INT
USE SCALAPACK_SUPPORT
IMPLICIT NONE
INCLUDE ‘mpif.h’
! Declare variables
INTEGER LDA, LDQR, NCA, NRA, DESCA(9), DESCL(9), DESCQ(9)
INTEGER INFO, MXCOL, MXLDA, LDQ
INTEGER, ALLOCATABLE :: IPVT(:), IPVT0(:)
LOGICAL PIVOT
REAL, ALLOCATABLE :: A(:,:), QR(:,:), Q(:,:), QRAUX(:)
REAL, ALLOCATABLE :: A0(:,:), QR0(:,:), Q0(:,:), QRAUX0(:)
PARAMETER (NRA=4, NCA=3, LDA=NRA, LDQR=NRA, LDQ=NRA)
! Set up for MPI
MP_NPROCS = MP_SETUP()
IF(MP_RANK .EQ. 0) THEN
ALLOCATE (A(NRA,NCA), Q(NRA,NRA), QR(NRA,NCA), &
QRAUX(NCA), IPVT(NCA))
! Set values for A and the righthand sides
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/)
!
IPVT = 0
ENDIF
! Set up a 1D processor grid and define
! its context ID, MP_ICTXT
CALL SCALAPACK_SETUP(NRA, NCA, .FALSE., .TRUE.)
! Get the array descriptor entities MXLDA,
! and MXCOL
CALL SCALAPACK_GETDIM(NRA, NCA, MP_MB, MP_NB, MXLDA, MXCOL)
! 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(DESCQ, NRA, NRA, MP_MB, MP_NB, 0, 0, MP_ICTXT, MXLDA, &
INFO)
! Allocate space for the local arrays
ALLOCATE (A0(MXLDA,MXCOL), QR0(MXLDA,MXCOL), QRAUX0(MXCOL), &
IPVT0(MXCOL), Q0(MXLDA,MXLDA))
! 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)
CALL LQERR (QR0, QRAUX0, Q0)
! 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(Q0, DESCQ, Q)
! Print results.
! Only Rank=0 has the solution, Q.
IF(MP_RANK .EQ. 0) CALL WRRRN (’Q’, Q)
! Exit Scalapack usage
CALL SCALAPACK_EXIT(MP_ICTXT)
! Shut down MPI
MP_NPROCS = MP_SETUP(‘FINAL’)
END
