1
SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
3
* -- LAPACK routine (version 2.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
8
* .. Scalar Arguments ..
12
* .. Array Arguments ..
14
DOUBLE PRECISION AP( * )
20
* DSPTRF computes the factorization of a real symmetric matrix A stored
21
* in packed format using the Bunch-Kaufman diagonal pivoting method:
23
* A = U*D*U**T or A = L*D*L**T
25
* where U (or L) is a product of permutation and unit upper (lower)
26
* triangular matrices, and D is symmetric and block diagonal with
27
* 1-by-1 and 2-by-2 diagonal blocks.
32
* UPLO (input) CHARACTER*1
33
* = 'U': Upper triangle of A is stored;
34
* = 'L': Lower triangle of A is stored.
37
* The order of the matrix A. N >= 0.
39
* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
40
* On entry, the upper or lower triangle of the symmetric matrix
41
* A, packed columnwise in a linear array. The j-th column of A
42
* is stored in the array AP as follows:
43
* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
44
* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
46
* On exit, the block diagonal matrix D and the multipliers used
47
* to obtain the factor U or L, stored as a packed triangular
48
* matrix overwriting A (see below for further details).
50
* IPIV (output) INTEGER array, dimension (N)
51
* Details of the interchanges and the block structure of D.
52
* If IPIV(k) > 0, then rows and columns k and IPIV(k) were
53
* interchanged and D(k,k) is a 1-by-1 diagonal block.
54
* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
55
* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
56
* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
57
* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
58
* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
60
* INFO (output) INTEGER
61
* = 0: successful exit
62
* < 0: if INFO = -i, the i-th argument had an illegal value
63
* > 0: if INFO = i, D(i,i) is exactly zero. The factorization
64
* has been completed, but the block diagonal matrix D is
65
* exactly singular, and division by zero will occur if it
66
* is used to solve a system of equations.
71
* If UPLO = 'U', then A = U*D*U', where
72
* U = P(n)*U(n)* ... *P(k)U(k)* ...,
73
* i.e., U is a product of terms P(k)*U(k), where k decreases from n to
74
* 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
75
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
76
* defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
77
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
84
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
85
* If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
86
* and A(k,k), and v overwrites A(1:k-2,k-1:k).
88
* If UPLO = 'L', then A = L*D*L', where
89
* L = P(1)*L(1)* ... *P(k)*L(k)* ...,
90
* i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
91
* n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
92
* and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
93
* defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
94
* that if the diagonal block D(k) is of order s (s = 1 or 2), then
101
* If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
102
* If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
103
* and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
105
* =====================================================================
108
DOUBLE PRECISION ZERO, ONE
109
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
110
DOUBLE PRECISION EIGHT, SEVTEN
111
PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
113
* .. Local Scalars ..
115
INTEGER IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC, KSTEP,
117
DOUBLE PRECISION ABSAKK, ALPHA, C, COLMAX, R1, R2, ROWMAX, S, T
119
* .. External Functions ..
122
EXTERNAL LSAME, IDAMAX
124
* .. External Subroutines ..
125
EXTERNAL DLAEV2, DROT, DSCAL, DSPR, DSWAP, XERBLA
127
* .. Intrinsic Functions ..
128
INTRINSIC ABS, MAX, SQRT
130
* .. Executable Statements ..
132
* Test the input parameters.
135
UPPER = LSAME( UPLO, 'U' )
136
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
138
ELSE IF( N.LT.0 ) THEN
142
CALL XERBLA( 'DSPTRF', -INFO )
146
* Initialize ALPHA for use in choosing pivot block size.
148
ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
152
* Factorize A as U*D*U' using the upper triangle of A
154
* K is the main loop index, decreasing from N to 1 in steps of
158
KC = ( N-1 )*N / 2 + 1
162
* If K < 1, exit from loop
168
* Determine rows and columns to be interchanged and whether
169
* a 1-by-1 or 2-by-2 pivot block will be used
171
ABSAKK = ABS( AP( KC+K-1 ) )
173
* IMAX is the row-index of the largest off-diagonal element in
174
* column K, and COLMAX is its absolute value
177
IMAX = IDAMAX( K-1, AP( KC ), 1 )
178
COLMAX = ABS( AP( KC+IMAX-1 ) )
183
IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
185
* Column K is zero: set INFO and continue
191
IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
193
* no interchange, use 1-by-1 pivot block
198
* JMAX is the column-index of the largest off-diagonal
199
* element in row IMAX, and ROWMAX is its absolute value
203
KX = IMAX*( IMAX+1 ) / 2 + IMAX
204
DO 20 J = IMAX + 1, K
205
IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
206
ROWMAX = ABS( AP( KX ) )
211
KPC = ( IMAX-1 )*IMAX / 2 + 1
213
JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
214
ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
217
IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
219
* no interchange, use 1-by-1 pivot block
222
ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
224
* interchange rows and columns K and IMAX, use 1-by-1
230
* interchange rows and columns K-1 and IMAX, use 2-by-2
243
* Interchange rows and columns KK and KP in the leading
244
* submatrix A(1:k,1:k)
246
CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
248
DO 30 J = KP + 1, KK - 1
251
AP( KNC+J-1 ) = AP( KX )
255
AP( KNC+KK-1 ) = AP( KPC+KP-1 )
257
IF( KSTEP.EQ.2 ) THEN
259
AP( KC+K-2 ) = AP( KC+KP-1 )
264
* Update the leading submatrix
266
IF( KSTEP.EQ.1 ) THEN
268
* 1-by-1 pivot block D(k): column k now holds
272
* where U(k) is the k-th column of U
274
* Perform a rank-1 update of A(1:k-1,1:k-1) as
276
* A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
278
R1 = ONE / AP( KC+K-1 )
279
CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
281
* Store U(k) in column k
283
CALL DSCAL( K-1, R1, AP( KC ), 1 )
286
* 2-by-2 pivot block D(k): columns k and k-1 now hold
288
* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
290
* where U(k) and U(k-1) are the k-th and (k-1)-th columns
293
* Perform a rank-2 update of A(1:k-2,1:k-2) as
295
* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
296
* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
298
* Convert this to two rank-1 updates by using the eigen-
299
* decomposition of D(k)
301
CALL DLAEV2( AP( KC-1 ), AP( KC+K-2 ), AP( KC+K-1 ), R1,
305
CALL DROT( K-2, AP( KNC ), 1, AP( KC ), 1, C, S )
306
CALL DSPR( UPLO, K-2, -R1, AP( KNC ), 1, AP )
307
CALL DSPR( UPLO, K-2, -R2, AP( KC ), 1, AP )
309
* Store U(k) and U(k-1) in columns k and k-1
311
CALL DSCAL( K-2, R1, AP( KNC ), 1 )
312
CALL DSCAL( K-2, R2, AP( KC ), 1 )
313
CALL DROT( K-2, AP( KNC ), 1, AP( KC ), 1, C, -S )
317
* Store details of the interchanges in IPIV
319
IF( KSTEP.EQ.1 ) THEN
326
* Decrease K and return to the start of the main loop
334
* Factorize A as L*D*L' using the lower triangle of A
336
* K is the main loop index, increasing from 1 to N in steps of
345
* If K > N, exit from loop
351
* Determine rows and columns to be interchanged and whether
352
* a 1-by-1 or 2-by-2 pivot block will be used
354
ABSAKK = ABS( AP( KC ) )
356
* IMAX is the row-index of the largest off-diagonal element in
357
* column K, and COLMAX is its absolute value
360
IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
361
COLMAX = ABS( AP( KC+IMAX-K ) )
366
IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
368
* Column K is zero: set INFO and continue
374
IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
376
* no interchange, use 1-by-1 pivot block
381
* JMAX is the column-index of the largest off-diagonal
382
* element in row IMAX, and ROWMAX is its absolute value
386
DO 50 J = K, IMAX - 1
387
IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
388
ROWMAX = ABS( AP( KX ) )
393
KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
395
JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
396
ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
399
IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
401
* no interchange, use 1-by-1 pivot block
404
ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
406
* interchange rows and columns K and IMAX, use 1-by-1
412
* interchange rows and columns K+1 and IMAX, use 2-by-2
422
$ KNC = KNC + N - K + 1
425
* Interchange rows and columns KK and KP in the trailing
426
* submatrix A(k:n,k:n)
429
$ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
432
DO 60 J = KK + 1, KP - 1
435
AP( KNC+J-KK ) = AP( KX )
439
AP( KNC ) = AP( KPC )
441
IF( KSTEP.EQ.2 ) THEN
443
AP( KC+1 ) = AP( KC+KP-K )
448
* Update the trailing submatrix
450
IF( KSTEP.EQ.1 ) THEN
452
* 1-by-1 pivot block D(k): column k now holds
456
* where L(k) is the k-th column of L
460
* Perform a rank-1 update of A(k+1:n,k+1:n) as
462
* A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
465
CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
468
* Store L(k) in column K
470
CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
474
* 2-by-2 pivot block D(k): columns K and K+1 now hold
476
* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
478
* where L(k) and L(k+1) are the k-th and (k+1)-th columns
483
* Perform a rank-2 update of A(k+2:n,k+2:n) as
485
* A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
486
* = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
488
* Convert this to two rank-1 updates by using the eigen-
489
* decomposition of D(k)
491
CALL DLAEV2( AP( KC ), AP( KC+1 ), AP( KNC ), R1, R2,
495
CALL DROT( N-K-1, AP( KC+2 ), 1, AP( KNC+1 ), 1, C,
497
CALL DSPR( UPLO, N-K-1, -R1, AP( KC+2 ), 1,
499
CALL DSPR( UPLO, N-K-1, -R2, AP( KNC+1 ), 1,
502
* Store L(k) and L(k+1) in columns k and k+1
504
CALL DSCAL( N-K-1, R1, AP( KC+2 ), 1 )
505
CALL DSCAL( N-K-1, R2, AP( KNC+1 ), 1 )
506
CALL DROT( N-K-1, AP( KC+2 ), 1, AP( KNC+1 ), 1, C,
512
* Store details of the interchanges in IPIV
514
IF( KSTEP.EQ.1 ) THEN
521
* Increase K and return to the start of the main loop
added
removed
Utilities/vxl/v3p/netlib/lapack/double/dsptrf.f
