LUPCH

Updates the RT R Cholesky factorization of a real symmetric positive definite matrix after a rank-one matrix is added.

Required Arguments

R — N by N upper triangular matrix containing the upper triangular factor to be updated. (Input)

Only the upper triangle of R is referenced.

Only the upper triangle of R is referenced.

X — Vector of length N determining the rank-one matrix to be added to the factorization RT R. (Input)

RNEW — N by N upper triangular matrix containing the updated triangular factor of RT R + XXT. (Output)

Only the upper triangle of RNEW is referenced. If R is not needed, R and RNEW can share the same storage locations.

Only the upper triangle of RNEW is referenced. If R is not needed, R and RNEW can share the same storage locations.

Optional Arguments

N — Order of the matrix. (Input)

Default: N = size (R,2).

Default: N = size (R,2).

LDR — Leading dimension of R exactly as specified in the dimension statement of the calling program. (Input)

Default: LDR = size (R,1).

Default: LDR = size (R,1).

LDRNEW — Leading dimension of RNEW exactly as specified in the dimension statement of the calling program. (Input)

Default: LDRNEW = size (RNEW,1).

Default: LDRNEW = size (RNEW,1).

CS — Vector of length N containing the cosines of the rotations. (Output)

SN — Vector of length N containing the sines of the rotations. (Output)

FORTRAN 90 Interface

Generic: CALL LUPCH (R, X, RNEW [, …])

Specific: The specific interface names are S_LUPCH and D_LUPCH.

FORTRAN 77 Interface

Single: CALL LUPCH (N, R, LDR, X, RNEW, LDRNEW, CS, SN)

Double: The double precision name is DLUPCH.

Description

The routine LUPCH is based on the LINPACK routine SCHUD; see Dongarra et al. (1979).

The Cholesky factorization of a matrix is A = RT R, where R is an upper triangular matrix. Given this factorization, LUPCH computes the factorization

In the program

is called RNEW.

LUPCH determines an orthogonal matrix U as the product GN…G1 of Givens rotations, such that

By multiplying this equation by its transpose, and noting that UT U = I, the desired result

is obtained.

Each Givens rotation, Gi, is chosen to zero out an element in xT. The matrix

Gi is (N + 1) × (N + 1) and has the form

Gi is (N + 1) × (N + 1) and has the form

Where Ik is the identity matrix of order k and ci = cosθi = CS(I), si = sinθi = SN(I) for some θi.

Example

A linear system Az = b is solved using the Cholesky factorization of A. This factorization is then updated and the system (A + xxT) z = b is solved using this updated factorization.

USE IMSL_LIBRARIES

! Declare variables

INTEGER LDA, LDFACT, N

PARAMETER (LDA=3, LDFACT=3, N=3)

REAL A(LDA,LDA), FACT(LDFACT,LDFACT), FACNEW(LDFACT,LDFACT), &

X(N), B(N), CS(N), SN(N), Z(N)

!

! Set values for A

! A = ( 1.0 -3.0 2.0)

! ( -3.0 10.0 -5.0)

! ( 2.0 -5.0 6.0)

!

DATA A/1.0, -3.0, 2.0, -3.0, 10.0, -5.0, 2.0, -5.0, 6.0/

!

! Set values for X and B

DATA X/3.0, 2.0, 1.0/

DATA B/53.0, 20.0, 31.0/

! Factor the matrix A

CALL LFTDS (A, FACT)

! Solve the original system

CALL LFSDS (FACT, B, Z)

! Print the results

CALL WRRRN (’FACT’, FACT, ITRING=1)

CALL WRRRN (’Z’, Z, 1, N, 1)

! Update the factorization

CALL LUPCH (FACT, X, FACNEW)

! Solve the updated system

CALL LFSDS (FACNEW, B, Z)

! Print the results

CALL WRRRN (’FACNEW’, FACNEW, ITRING=1)

CALL WRRRN (’Z’, Z, 1, N, 1)

!

END

Output

FACT

1 2 3

1 1.000 -3.000 2.000

2 1.000 1.000

3 1.000

Z

1 2 3

1860.0 433.0 -254.0

FACNEW

1 2 3

1 3.162 0.949 1.581

2 3.619 -1.243

3 -1.719

Z

1 2 3

4.000 1.000 2.000