Update an Augmented Cholesky Decomposition

Calling Sequence

CALL DCHUD	(DA, LDA, N, DX, DZ, LDZ, NZ, DY, DRHO, DCOS, DSIN)
CALL SCHUD	(SA, LDA, N, SX, SZ, LDZ, NZ, SY, SRHO, SCOS, SSIN)
CALL ZCHUD	(ZA, LDA, N, ZX, ZZ, LDZ, NZ, ZY, DRHO, DCOS, ZSIN)
CALL CCHUD	(CA, LDA, N, CX, CZ, LDZ, NZ, CY, SRHO, SCOS, CSIN)
void dchud	(double *da, int lda, int n, double *dx, double *dz, int ldz, int nz, double *dy, double *drho, double *dcos, double *dsin)
void schud	(float *sa, int lda, int n, float *sx, float *sz, int ldz, int nz, float *sy, float *srho, float *scos, float *ssin)
void zchud	(doublecomplex *za, int lda, int n, doublecomplex *zx, doublecomplex *zz, int ldz, int nz, doublecomplex *zy, double *drho, double *dcos, doublecomplex *zsin)
void cchud	(complex *ca, int lda, int n, complex *cx, complex *cz, int ldz, int nz, complex *cy, float *srho, float *scos, complex *csin)

Arguments

xA	On entry, the upper triangular matrix A. On exit, A has been updated. The strict lower triangle of A is not referenced.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1,N).
N	Order of the matrix A. N 0.
xX	Row to be added to A.
xZ	Vectors to be updated with A.
LDZ	Leading dimension on the array Z as specified in a dimension or type statement. LDZ max(1,N).
NZ	Number of vectors to be updated with A. NZ 0. If NZ = 0 then Z, Y, and RHO are not used.
xY	Scalars for updating the vectors in Z.
xRHO	On entry, the norms of the residual vectors that are to be updated. On exit, RHO has been updated. If RHO(i) is negative on entry then it is not changed.
xCOS	Cosines of the transforming rotations.
xSIN	Sines of the transforming rotations.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, N, NOPIV, NZ

PARAMETER (N = 4)

PARAMETER (NOPIV = 0)

PARAMETER (NZ = 0)

PARAMETER (LDA = N)

C

DOUBLE PRECISION A(LDA,N), ANULL, C(N), S(N), WORK(N), X(N)

INTEGER I, ICOL, INFO, IPIVOT(N), IROW, JOB, NULL

C

C Initialize the arrays A and Z to store the matrices A and Z

C shown below and initialize X and Y to store the vectors x

C and y shown below.

C

C 4 3 2 1 1

C A = 3 4 3 2 x = 1

C 2 3 4 3 1

C 1 2 3 4 1

C

DATA A / 4.0D0, 3*8D8, 3.0D0, 4.0D0, 2*8D8, 2.0D0, 3.0D0,

$ 4.0D0, 8D8, 1.0D0, 2.0D0, 3.0D0, 4.0D0 /

C

PRINT 1000

PRINT 1010, A(1,1), A(1,2), A(1,3), A(1,4)

PRINT 1010, A(1,2), A(2,2), A(2,3), A(2,4)

PRINT 1010, A(1,3), A(2,3), A(3,3), A(3,4)

PRINT 1010, (A(IROW,4), IROW = 1, N)

PRINT 1020

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

JOB = NOPIV

CALL DCHDC (A, LDA, N, WORK, IPIVOT, JOB, INFO)

IF (INFO .EQ. N) THEN

PRINT 1030

PRINT 1010, A(1,1), A(1,2), A(1,3), A(1,4)

PRINT 1040, A(2,2), A(2,3), A(2,4)

PRINT 1050, A(3,3), A(3,4)

PRINT 1060, A(4,4)

ANULL = 0.0D0

NULL = 1

CALL DCHUD (A, LDA, N, X, ANULL, NULL, NZ, ANULL, ANULL,

$ C, S)

PRINT 1070

PRINT 1080, (C(I), S(I), I = 1, N)

ELSE

PRINT 1090

END IF

C

1000 FORMAT (1X, 'A in full form:')

1010 FORMAT (4(3X, F7.3))

1020 FORMAT (/1X, 'A in symmetric form (* in unused entries)')

1030 FORMAT (/1X, 'Upper Cholesky factor:')

1040 FORMAT (10X, 3(3X, F7.3))

1050 FORMAT (20X, 2(3X, F7.3))

1060 FORMAT (30X, 1(3X, F7.3))

1070 FORMAT (1X, 'Cosine', 3X, ' Sine')

1080 FORMAT (1X, F6.3, 3X, F6.3)

1090 FORMAT (/1X, 'A is not positive definite.')

C

END

Sample Output

A in full form:
4.000 3.000 2.000 1.000
3.000 4.000 3.000 2.000
2.000 3.000 4.000 3.000
1.000 2.000 3.000 4.000
A in symmetric form (* in unused entries)
4.000 3.000 2.000 1.000
******* 4.000 3.000 2.000
******* ******* 4.000 3.000
******* ******* ******* 4.000
Upper Cholesky factor:
2.000 1.500 1.000 0.500
1.323 1.134 0.945
1.309 1.091
1.291
Cosine Sine
1.000 0.000
1.000 0.000
1.000 0.000
1.000 0.000