GSWEP
Performs a generalized sweep of a row of a nonnegative definite matrix.
Required Arguments
KROW — Row/column number to be swept. (Input)
A — N by N nonnegative definite matrix whose row KROW is to be swept. (Input/Output)
Only the upper triangle of A is referenced.
Optional Arguments
N — Order of the matrix to be swept. (Input)
Default: N = size (A,2).
LDA — Leading dimension of A exactly as specified in the dimension statement in the calling program. (Input)
Default: LDA = size (A,1).
IREV — Reversibility option. (Input)
Default: IREV = 0.
IREV | Action When Linear Dependence Is Declared |
|---|
0 | Elements of row and column KROW of A are set to 0.0. Reversibility of the generalized sweep operator is lost. |
1 | Elements of row and column KROW of A are left unchanged. Reversibility of the generalized sweep operator is maintained, but some post processing by the user is required. See Comments. |
TOL — Tolerance used in determining linear dependence. (Input)
TOL = 100
* AMACH(4) is a common choice. See documentation for routine
AMACH in the
Reference Material.
Default:
TOL = 1.e-5 for single precision and 2.d –14 for double precision.
SCALE — Vector of length N containing the diagonal scaling matrix used in the tolerance check. (Input)
A common choice for SCALE(I) is the I-th diagonal element of A before any calls to GSWEP have been made. If TOL = 0.0, SCALE is not referenced and can be a vector of length one.
SWEPT — Vector of length N with information to indicate what has and has not been swept. (Input/Output)
On the first call to GSWEP all elements must equal ‑1.0. On output, SWEPT(KROW) = 1.0 if the sweep was successful. If a linear dependence is declared, SWEPT(KROW) remains equal to ‑1.0.
FORTRAN 90 Interface
Generic: CALL GSWEP (KROW, A [, …])
Specific: The specific interface names are S_GSWEP and D_GSWEP.
FORTRAN 77 Interface
Single: CALL GSWEP (KROW, N, A, LDA, IREV, TOL, SCALE, SWEPT)
Double: The double precision name is DGSWEP.
Description
Routine GSWEP computes an upper triangular generalized sweep of a nonnegative definite matrix. The versatility of the SWEEP operator for statistical computations, in particular for regression computations, is discussed by Goodnight (1979).
Routine GSWEP is based on UTG2SWEEP and RUTG2SWEEP described by Goodnight (1979, pages 157-158). (A misprint appears twice in “Step 5”, page 157 of Goodnight’s article. The “aij” should be replaced by “aik.”) The test for linear dependence is the same as that given by Clarke (1982).
Comments
Say we wish to sweep k different rows of the matrix A. For purposes of discussion, let these be rows 1, 2, …, k of A. Partition A into its first k rows and columns and the remainder,
For a nonsingular A11, successive invocations of GSWEP with A and KROW equal to 1, 2, …, k yields
Only the elements in the upper triangle of A are referenced. Thus, the elements in the lower triangles of the symmetric matrices
are not returned. For a singular A11and IREV equal to zero, a symmetric g2 inverse of A11, denoted by
is used. For a singular A11and IREV not equal to zero, the first k rows of the swept A are not the same as for the IREV equal to one case. However,
can be obtained from the output A as follows:
and
H is the Hermite canonical form (also referred to as the Hermite normal form or a rowechelon form) of A11.
Example
We consider the correlation matrix for the first three regressors from the example used by Berk (1976) and discussed by Frane (1977). The matrix is “nearly” singular. The rows of the correlation matrix are swept sequentially with KROW equal 1, 2, 3. With a tolerance of 0.001, the sweeps for 1 and 2 are successful. When a sweep on row 3 is attempted a linear dependence is declared. This is because
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDA, N
PARAMETER (N=3, LDA=N)
!
INTEGER ISETNG, KROW
REAL A(LDA,N), SCALE(N), SQRT, SWEPT(N), TOL
INTRINSIC SQRT
!
A(1,1) = 1.0
A(1,2) = SQRT(0.99)
A(1,3) = 0.1*SQRT(0.99)
A(2,2) = 1.0
A(2,3) = 0.0
A(3,3) = 1.0
TOL = 0.001
! Copy diagonal of A to SCALE.
CALL SCOPY (N, A(1:,1), LDA+1, SCALE, 1)
! Initialize elements of SWEPT to -1.
SWEPT = -1.0
ISETNG = 4
CALL WROPT (-6, ISETNG, 1)
CALL WRRRN ('A', A, ITRING=1)
CALL WRRRN ('SWEPT', SWEPT)
SWEPT = -1.0
DO 10 KROW=1, 3
CALL GSWEP (KROW, A, tol=tol, scale=scale, swept=swept)
CALL WRRRN ('A', A, ITRING=1)
CALL WRRRN ('SWEPT', SWEPT)
10 CONTINUE
END
Output
A
1 2 3
1 1.00000 0.99499 0.09950
2 1.00000 0.00000
3 1.00000
SWEPT
1 -1.00000
2 -1.00000
3 -1.00000
A
1 2 3
1 1.00000 0.99499 0.09950
2 0.01000 -0.09900
3 0.99010
SWEPT
1 1.00000
2 -1.00000
3 -1.00000
A
1 2 3
1 100.000 -99.499 9.950
2 100.000 -9.900
3 0.010
SWEPT
1 1.00000
2 1.00000
3 -1.00000
A
1 2 3
1 100.000 -99.499 0.000
2 100.000 0.000
3 0.000
SWEPT
1 1.00000
2 1.00000
3 -1.00000
