LUPQR

Computes an updated QR factorization after the rank-one matrix αxyT is added.

Required Arguments

ALPHA — Scalar determining the rank-one update to be added. (Input)

W — Vector of length NROW determining the rank-one matrix to be added. (Input)
The updated matrix is A + αxyT. If I = 0 then W contains the vector x. If I = 1 then W contains the vector QTx.

Y — Vector of length NCOL determining the rank-one matrix to be added. (Input)

R — Matrix of order NROW by NCOL containing the R matrix from the QR factorization. (Input)
Only the upper trapezoidal part of R is referenced.

IPATH — Flag used to control the computation of the QR update. (Input)
IPATH has the decimal expansion IJ such that: I = 0 means W contains the vector x.
I= 1 means W contains the vector QTx.
J = 0 means do not update the matrix Q. J = 1 means update the matrix Q. For example, if IPATH = 10 then, I = 1 and J = 0.

RNEW — Matrix of order NROW by NCOL containing the updated R matrix in the QR factorization. (Output)
Only the upper trapezoidal part of RNEW is updated. R and RNEW may be the same.

Optional Arguments

NROW — Number of rows in the matrix A = Q * R. (Input)
Default: NROW = size (W,1).

NCOL — Number of columns in the matrix A = Q * R. (Input)
Default: NCOL = size (Y,1).

Q — Matrix of order NROW containing the Q matrix from the QR factorization. (Input)
Ignored if IPATH = 0.
Default: Q is 1x1 and un-initialized.

LDQ — Leading dimension of Q exactly as specified in the dimension statement of the calling program. (Input)
Ignored if IPATH = 0.
Default: LDQ = size (Q,1).

LDR — Leading dimension of R exactly as specified in the dimension statement of the calling program. (Input)
Default: LDR = size (R,1).

QNEW — Matrix of order NROW containing the updated Q matrix in the QR factorization. (Output)
Ignored if J = 0. See IPATH for a definition of J.

LDQNEW — Leading dimension of QNEW exactly as specified in the dimension statement of the calling program. (Input)
Ignored if J = 0. See IPATH for a definition of J.
Default: LDQNEW = size (QNEW,1).

LDRNEW — Leading dimension of RNEW exactly as specified in the dimension statement of the calling program. (Input)
Default: LDRNEW = size (RNEW,1).

FORTRAN 90 Interface

Generic: CALL LUPQR (ALPHA, W, Y, R, IPATH, RNEW [, …])

Specific: The specific interface names are S_LUPQR and D_LUPQR.

FORTRAN 77 Interface

Single: CALL LUPQR (NROW, NCOL, ALPHA, W, Y, Q, LDQ, R, LDR, IPATH, QNEW, LDQNEW, RNEW, LDRNEW)

Double: The double precision name is DLUPQR.

Description

Let A be an m × n matrix and let A = QR be its QR decomposition. (In the program, m is called NROW and n is called NCOL) Then

A + αxyT = QR + αxyT = Q(R + αQTxyT) = Q(R + αwyT)

where w = QT x. An orthogonal transformation J can be constructed, using a sequence of m − 1 Givens rotations, such that Jw = ωe1, where ω = ±∥w∥2 and e1 = (1, 0, …, 0)T. Then

A + αxyT = (QJT)(JR + αωe1yT)

Since JR is an upper Hessenberg matrix, H = JR + αωe1yT is also an upper Hessenberg matrix. Again using m − 1 Givens rotations, an orthogonal transformation G can be constructed such that GH is an upper triangular matrix. Then

 

is orthogonal and

 

is upper triangular.

If the last k components of w are zero, then the number of Givens rotations needed to construct J or G is m − k − 1 instead of m − 1.

For further information, see Dennis and Schnabel (1983, pages 55-58 and 311-313), or Golub and Van Loan (1983, pages 437− 439).

Comments

1. Workspace may be explicitly provided, if desired, by use of L2PQR/DL2PQR. The reference is:

CALL L2PQR (NROW, NCOL, ALPHA, W, Y, Q, LDQ, R, LDR, IPATH, QNEW, LDQNEW, RNEW, LDRNEW, Z, WORK)

The additional arguments are as follows:

Z — Work vector of length NROW.

WORK — Work vector of length MIN(NROW − 1, NCOL).

Example

The QR factorization of A is found. It is then used to find the QR factorization of A + xyT. Since pivoting is used, the QR factorization routine finds AP = QR, where P is a permutation matrix determined by IPVT. We compute

 

The IMSL routine PERMU (see Utilities) is used to compute Py. As a check

 

is computed and printed. It can also be obtained from A + xyT by permuting its columns using the order given by IPVT.

 

USE IMSL_LIBRARIES

! Declare variables

INTEGER LDA, LDAQR, LDQ, LDQNEW, LDQR, LDR, LDRNEW, NCOL, NROW

PARAMETER (NCOL=3, NROW=4, LDA=NROW, LDAQR=NROW, LDQ=NROW, &

LDQNEW=NROW, LDQR=NROW, LDR=NROW, LDRNEW=NROW)

!

INTEGER IPATH, IPVT(NCOL), J, MIN0

REAL A(LDA,NCOL), ALPHA, AQR(LDAQR,NCOL), CONORM(NCOL), &

Q(LDQ,NROW), QNEW(LDQNEW,NROW), QR(LDQR,NCOL), &

QRAUX(NCOL), R(LDR,NCOL), RNEW(LDRNEW,NCOL), W(NROW), &

Y(NCOL)

LOGICAL PIVOT

INTRINSIC MIN0

!

! 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 W and Y

DATA W/1., 2., 3., 4./

DATA Y/3., 2., 1./

!

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

! Permute Y

CALL PERMU (Y, IPVT, Y)

! R is the upper trapezoidal part of QR

R = 0.0E0

DO 10 J=1, NCOL

CALL SCOPY (MIN0(J,NROW), QR(:,J), 1, R(:,J), 1)

10 CONTINUE

! Update Q and R

ALPHA = 1.0

IPATH = 01

CALL LUPQR (ALPHA, W, Y, R, IPATH, RNEW, Q=Q, QNEW=QNEW)

! Compute AQR = Q*R

CALL MRRRR (QNEW, RNEW, AQR)

! Print results

CALL WRIRN (’IPVT’, IPVT, 1, NCOL,1)

CALL WRRRN (’QNEW’, QNEW)

CALL WRRRN (’RNEW’, RNEW)

CALL WRRRN (’QNEW*RNEW’, AQR)

END

Output

 

IPVT

1 2 3

3 2 1

 

QNEW

1 2 3 4

1 -0.0620 -0.5412 0.8082 -0.2236

2 -0.2234 -0.6539 -0.2694 0.6708

3 -0.4840 -0.3379 -0.4490 -0.6708

4 -0.8438 0.4067 0.2694 0.2236

 

RNEW

1 2 3

1 -80.59 -21.34 -17.62

2 0.00 -4.94 -4.83

3 0.00 0.00 0.36

4 0.00 0.00 0.00

 

QNEW*RNEW

1 2 3

1 5.00 4.00 4.00

2 18.00 8.00 7.00

3 39.00 12.00 10.00

4 68.00 16.00 13.00