1 SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRF computes the factorization of a real symmetric matrix A stored
21 * in packed format using the Bunch-Kaufman diagonal pivoting method:
22 *
23 * A = U*D*U**T or A = L*D*L**T
24 *
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.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
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.
45 *
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).
49 *
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.
59 *
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.
67 *
68 * Further Details
69 * ===============
70 *
71 * 5-96 - Based on modifications by J. Lewis, Boeing Computer Services
72 * Company
73 *
74 * If UPLO = 'U', then A = U*D*U**T, where
75 * U = P(n)*U(n)* ... *P(k)U(k)* ...,
76 * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
77 * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
78 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
79 * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
80 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
81 *
82 * ( I v 0 ) k-s
83 * U(k) = ( 0 I 0 ) s
84 * ( 0 0 I ) n-k
85 * k-s s n-k
86 *
87 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
88 * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
89 * and A(k,k), and v overwrites A(1:k-2,k-1:k).
90 *
91 * If UPLO = 'L', then A = L*D*L**T, where
92 * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
93 * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
94 * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
95 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
96 * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
97 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
98 *
99 * ( I 0 0 ) k-1
100 * L(k) = ( 0 I 0 ) s
101 * ( 0 v I ) n-k-s+1
102 * k-1 s n-k-s+1
103 *
104 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
105 * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
106 * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
107 *
108 * =====================================================================
109 *
110 * .. Parameters ..
111 DOUBLE PRECISION ZERO, ONE
112 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
113 DOUBLE PRECISION EIGHT, SEVTEN
114 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL UPPER
118 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
119 $ KSTEP, KX, NPP
120 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
121 $ ROWMAX, T, WK, WKM1, WKP1
122 * ..
123 * .. External Functions ..
124 LOGICAL LSAME
125 INTEGER IDAMAX
126 EXTERNAL LSAME, IDAMAX
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL DSCAL, DSPR, DSWAP, XERBLA
130 * ..
131 * .. Intrinsic Functions ..
132 INTRINSIC ABS, MAX, SQRT
133 * ..
134 * .. Executable Statements ..
135 *
136 * Test the input parameters.
137 *
138 INFO = 0
139 UPPER = LSAME( UPLO, 'U' )
140 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
141 INFO = -1
142 ELSE IF( N.LT.0 ) THEN
143 INFO = -2
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'DSPTRF', -INFO )
147 RETURN
148 END IF
149 *
150 * Initialize ALPHA for use in choosing pivot block size.
151 *
152 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
153 *
154 IF( UPPER ) THEN
155 *
156 * Factorize A as U*D*U**T using the upper triangle of A
157 *
158 * K is the main loop index, decreasing from N to 1 in steps of
159 * 1 or 2
160 *
161 K = N
162 KC = ( N-1 )*N / 2 + 1
163 10 CONTINUE
164 KNC = KC
165 *
166 * If K < 1, exit from loop
167 *
168 IF( K.LT.1 )
169 $ GO TO 110
170 KSTEP = 1
171 *
172 * Determine rows and columns to be interchanged and whether
173 * a 1-by-1 or 2-by-2 pivot block will be used
174 *
175 ABSAKK = ABS( AP( KC+K-1 ) )
176 *
177 * IMAX is the row-index of the largest off-diagonal element in
178 * column K, and COLMAX is its absolute value
179 *
180 IF( K.GT.1 ) THEN
181 IMAX = IDAMAX( K-1, AP( KC ), 1 )
182 COLMAX = ABS( AP( KC+IMAX-1 ) )
183 ELSE
184 COLMAX = ZERO
185 END IF
186 *
187 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
188 *
189 * Column K is zero: set INFO and continue
190 *
191 IF( INFO.EQ.0 )
192 $ INFO = K
193 KP = K
194 ELSE
195 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
196 *
197 * no interchange, use 1-by-1 pivot block
198 *
199 KP = K
200 ELSE
201 *
202 ROWMAX = ZERO
203 JMAX = IMAX
204 KX = IMAX*( IMAX+1 ) / 2 + IMAX
205 DO 20 J = IMAX + 1, K
206 IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
207 ROWMAX = ABS( AP( KX ) )
208 JMAX = J
209 END IF
210 KX = KX + J
211 20 CONTINUE
212 KPC = ( IMAX-1 )*IMAX / 2 + 1
213 IF( IMAX.GT.1 ) THEN
214 JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
215 ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
216 END IF
217 *
218 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
219 *
220 * no interchange, use 1-by-1 pivot block
221 *
222 KP = K
223 ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
224 *
225 * interchange rows and columns K and IMAX, use 1-by-1
226 * pivot block
227 *
228 KP = IMAX
229 ELSE
230 *
231 * interchange rows and columns K-1 and IMAX, use 2-by-2
232 * pivot block
233 *
234 KP = IMAX
235 KSTEP = 2
236 END IF
237 END IF
238 *
239 KK = K - KSTEP + 1
240 IF( KSTEP.EQ.2 )
241 $ KNC = KNC - K + 1
242 IF( KP.NE.KK ) THEN
243 *
244 * Interchange rows and columns KK and KP in the leading
245 * submatrix A(1:k,1:k)
246 *
247 CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
248 KX = KPC + KP - 1
249 DO 30 J = KP + 1, KK - 1
250 KX = KX + J - 1
251 T = AP( KNC+J-1 )
252 AP( KNC+J-1 ) = AP( KX )
253 AP( KX ) = T
254 30 CONTINUE
255 T = AP( KNC+KK-1 )
256 AP( KNC+KK-1 ) = AP( KPC+KP-1 )
257 AP( KPC+KP-1 ) = T
258 IF( KSTEP.EQ.2 ) THEN
259 T = AP( KC+K-2 )
260 AP( KC+K-2 ) = AP( KC+KP-1 )
261 AP( KC+KP-1 ) = T
262 END IF
263 END IF
264 *
265 * Update the leading submatrix
266 *
267 IF( KSTEP.EQ.1 ) THEN
268 *
269 * 1-by-1 pivot block D(k): column k now holds
270 *
271 * W(k) = U(k)*D(k)
272 *
273 * where U(k) is the k-th column of U
274 *
275 * Perform a rank-1 update of A(1:k-1,1:k-1) as
276 *
277 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
278 *
279 R1 = ONE / AP( KC+K-1 )
280 CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
281 *
282 * Store U(k) in column k
283 *
284 CALL DSCAL( K-1, R1, AP( KC ), 1 )
285 ELSE
286 *
287 * 2-by-2 pivot block D(k): columns k and k-1 now hold
288 *
289 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
290 *
291 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
292 * of U
293 *
294 * Perform a rank-2 update of A(1:k-2,1:k-2) as
295 *
296 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
297 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
298 *
299 IF( K.GT.2 ) THEN
300 *
301 D12 = AP( K-1+( K-1 )*K / 2 )
302 D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
303 D11 = AP( K+( K-1 )*K / 2 ) / D12
304 T = ONE / ( D11*D22-ONE )
305 D12 = T / D12
306 *
307 DO 50 J = K - 2, 1, -1
308 WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
309 $ AP( J+( K-1 )*K / 2 ) )
310 WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
311 $ AP( J+( K-2 )*( K-1 ) / 2 ) )
312 DO 40 I = J, 1, -1
313 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
314 $ AP( I+( K-1 )*K / 2 )*WK -
315 $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1
316 40 CONTINUE
317 AP( J+( K-1 )*K / 2 ) = WK
318 AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
319 50 CONTINUE
320 *
321 END IF
322 *
323 END IF
324 END IF
325 *
326 * Store details of the interchanges in IPIV
327 *
328 IF( KSTEP.EQ.1 ) THEN
329 IPIV( K ) = KP
330 ELSE
331 IPIV( K ) = -KP
332 IPIV( K-1 ) = -KP
333 END IF
334 *
335 * Decrease K and return to the start of the main loop
336 *
337 K = K - KSTEP
338 KC = KNC - K
339 GO TO 10
340 *
341 ELSE
342 *
343 * Factorize A as L*D*L**T using the lower triangle of A
344 *
345 * K is the main loop index, increasing from 1 to N in steps of
346 * 1 or 2
347 *
348 K = 1
349 KC = 1
350 NPP = N*( N+1 ) / 2
351 60 CONTINUE
352 KNC = KC
353 *
354 * If K > N, exit from loop
355 *
356 IF( K.GT.N )
357 $ GO TO 110
358 KSTEP = 1
359 *
360 * Determine rows and columns to be interchanged and whether
361 * a 1-by-1 or 2-by-2 pivot block will be used
362 *
363 ABSAKK = ABS( AP( KC ) )
364 *
365 * IMAX is the row-index of the largest off-diagonal element in
366 * column K, and COLMAX is its absolute value
367 *
368 IF( K.LT.N ) THEN
369 IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
370 COLMAX = ABS( AP( KC+IMAX-K ) )
371 ELSE
372 COLMAX = ZERO
373 END IF
374 *
375 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
376 *
377 * Column K is zero: set INFO and continue
378 *
379 IF( INFO.EQ.0 )
380 $ INFO = K
381 KP = K
382 ELSE
383 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
384 *
385 * no interchange, use 1-by-1 pivot block
386 *
387 KP = K
388 ELSE
389 *
390 * JMAX is the column-index of the largest off-diagonal
391 * element in row IMAX, and ROWMAX is its absolute value
392 *
393 ROWMAX = ZERO
394 KX = KC + IMAX - K
395 DO 70 J = K, IMAX - 1
396 IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
397 ROWMAX = ABS( AP( KX ) )
398 JMAX = J
399 END IF
400 KX = KX + N - J
401 70 CONTINUE
402 KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
403 IF( IMAX.LT.N ) THEN
404 JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
405 ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
406 END IF
407 *
408 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
409 *
410 * no interchange, use 1-by-1 pivot block
411 *
412 KP = K
413 ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
414 *
415 * interchange rows and columns K and IMAX, use 1-by-1
416 * pivot block
417 *
418 KP = IMAX
419 ELSE
420 *
421 * interchange rows and columns K+1 and IMAX, use 2-by-2
422 * pivot block
423 *
424 KP = IMAX
425 KSTEP = 2
426 END IF
427 END IF
428 *
429 KK = K + KSTEP - 1
430 IF( KSTEP.EQ.2 )
431 $ KNC = KNC + N - K + 1
432 IF( KP.NE.KK ) THEN
433 *
434 * Interchange rows and columns KK and KP in the trailing
435 * submatrix A(k:n,k:n)
436 *
437 IF( KP.LT.N )
438 $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
439 $ 1 )
440 KX = KNC + KP - KK
441 DO 80 J = KK + 1, KP - 1
442 KX = KX + N - J + 1
443 T = AP( KNC+J-KK )
444 AP( KNC+J-KK ) = AP( KX )
445 AP( KX ) = T
446 80 CONTINUE
447 T = AP( KNC )
448 AP( KNC ) = AP( KPC )
449 AP( KPC ) = T
450 IF( KSTEP.EQ.2 ) THEN
451 T = AP( KC+1 )
452 AP( KC+1 ) = AP( KC+KP-K )
453 AP( KC+KP-K ) = T
454 END IF
455 END IF
456 *
457 * Update the trailing submatrix
458 *
459 IF( KSTEP.EQ.1 ) THEN
460 *
461 * 1-by-1 pivot block D(k): column k now holds
462 *
463 * W(k) = L(k)*D(k)
464 *
465 * where L(k) is the k-th column of L
466 *
467 IF( K.LT.N ) THEN
468 *
469 * Perform a rank-1 update of A(k+1:n,k+1:n) as
470 *
471 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
472 *
473 R1 = ONE / AP( KC )
474 CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
475 $ AP( KC+N-K+1 ) )
476 *
477 * Store L(k) in column K
478 *
479 CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
480 END IF
481 ELSE
482 *
483 * 2-by-2 pivot block D(k): columns K and K+1 now hold
484 *
485 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
486 *
487 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
488 * of L
489 *
490 IF( K.LT.N-1 ) THEN
491 *
492 * Perform a rank-2 update of A(k+2:n,k+2:n) as
493 *
494 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
495 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
496 *
497 * where L(k) and L(k+1) are the k-th and (k+1)-th
498 * columns of L
499 *
500 D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 )
501 D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
502 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
503 T = ONE / ( D11*D22-ONE )
504 D21 = T / D21
505 *
506 DO 100 J = K + 2, N
507 WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
508 $ AP( J+K*( 2*N-K-1 ) / 2 ) )
509 WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
510 $ AP( J+( K-1 )*( 2*N-K ) / 2 ) )
511 *
512 DO 90 I = J, N
513 AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
514 $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
515 $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1
516 90 CONTINUE
517 *
518 AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
519 AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
520 *
521 100 CONTINUE
522 END IF
523 END IF
524 END IF
525 *
526 * Store details of the interchanges in IPIV
527 *
528 IF( KSTEP.EQ.1 ) THEN
529 IPIV( K ) = KP
530 ELSE
531 IPIV( K ) = -KP
532 IPIV( K+1 ) = -KP
533 END IF
534 *
535 * Increase K and return to the start of the main loop
536 *
537 K = K + KSTEP
538 KC = KNC + N - K + 2
539 GO TO 60
540 *
541 END IF
542 *
543 110 CONTINUE
544 RETURN
545 *
546 * End of DSPTRF
547 *
548 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRF computes the factorization of a real symmetric matrix A stored
21 * in packed format using the Bunch-Kaufman diagonal pivoting method:
22 *
23 * A = U*D*U**T or A = L*D*L**T
24 *
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.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
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.
45 *
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).
49 *
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.
59 *
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.
67 *
68 * Further Details
69 * ===============
70 *
71 * 5-96 - Based on modifications by J. Lewis, Boeing Computer Services
72 * Company
73 *
74 * If UPLO = 'U', then A = U*D*U**T, where
75 * U = P(n)*U(n)* ... *P(k)U(k)* ...,
76 * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
77 * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
78 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
79 * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
80 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
81 *
82 * ( I v 0 ) k-s
83 * U(k) = ( 0 I 0 ) s
84 * ( 0 0 I ) n-k
85 * k-s s n-k
86 *
87 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
88 * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
89 * and A(k,k), and v overwrites A(1:k-2,k-1:k).
90 *
91 * If UPLO = 'L', then A = L*D*L**T, where
92 * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
93 * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
94 * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
95 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
96 * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
97 * that if the diagonal block D(k) is of order s (s = 1 or 2), then
98 *
99 * ( I 0 0 ) k-1
100 * L(k) = ( 0 I 0 ) s
101 * ( 0 v I ) n-k-s+1
102 * k-1 s n-k-s+1
103 *
104 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
105 * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
106 * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
107 *
108 * =====================================================================
109 *
110 * .. Parameters ..
111 DOUBLE PRECISION ZERO, ONE
112 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
113 DOUBLE PRECISION EIGHT, SEVTEN
114 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
115 * ..
116 * .. Local Scalars ..
117 LOGICAL UPPER
118 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
119 $ KSTEP, KX, NPP
120 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
121 $ ROWMAX, T, WK, WKM1, WKP1
122 * ..
123 * .. External Functions ..
124 LOGICAL LSAME
125 INTEGER IDAMAX
126 EXTERNAL LSAME, IDAMAX
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL DSCAL, DSPR, DSWAP, XERBLA
130 * ..
131 * .. Intrinsic Functions ..
132 INTRINSIC ABS, MAX, SQRT
133 * ..
134 * .. Executable Statements ..
135 *
136 * Test the input parameters.
137 *
138 INFO = 0
139 UPPER = LSAME( UPLO, 'U' )
140 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
141 INFO = -1
142 ELSE IF( N.LT.0 ) THEN
143 INFO = -2
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'DSPTRF', -INFO )
147 RETURN
148 END IF
149 *
150 * Initialize ALPHA for use in choosing pivot block size.
151 *
152 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
153 *
154 IF( UPPER ) THEN
155 *
156 * Factorize A as U*D*U**T using the upper triangle of A
157 *
158 * K is the main loop index, decreasing from N to 1 in steps of
159 * 1 or 2
160 *
161 K = N
162 KC = ( N-1 )*N / 2 + 1
163 10 CONTINUE
164 KNC = KC
165 *
166 * If K < 1, exit from loop
167 *
168 IF( K.LT.1 )
169 $ GO TO 110
170 KSTEP = 1
171 *
172 * Determine rows and columns to be interchanged and whether
173 * a 1-by-1 or 2-by-2 pivot block will be used
174 *
175 ABSAKK = ABS( AP( KC+K-1 ) )
176 *
177 * IMAX is the row-index of the largest off-diagonal element in
178 * column K, and COLMAX is its absolute value
179 *
180 IF( K.GT.1 ) THEN
181 IMAX = IDAMAX( K-1, AP( KC ), 1 )
182 COLMAX = ABS( AP( KC+IMAX-1 ) )
183 ELSE
184 COLMAX = ZERO
185 END IF
186 *
187 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
188 *
189 * Column K is zero: set INFO and continue
190 *
191 IF( INFO.EQ.0 )
192 $ INFO = K
193 KP = K
194 ELSE
195 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
196 *
197 * no interchange, use 1-by-1 pivot block
198 *
199 KP = K
200 ELSE
201 *
202 ROWMAX = ZERO
203 JMAX = IMAX
204 KX = IMAX*( IMAX+1 ) / 2 + IMAX
205 DO 20 J = IMAX + 1, K
206 IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
207 ROWMAX = ABS( AP( KX ) )
208 JMAX = J
209 END IF
210 KX = KX + J
211 20 CONTINUE
212 KPC = ( IMAX-1 )*IMAX / 2 + 1
213 IF( IMAX.GT.1 ) THEN
214 JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
215 ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
216 END IF
217 *
218 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
219 *
220 * no interchange, use 1-by-1 pivot block
221 *
222 KP = K
223 ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
224 *
225 * interchange rows and columns K and IMAX, use 1-by-1
226 * pivot block
227 *
228 KP = IMAX
229 ELSE
230 *
231 * interchange rows and columns K-1 and IMAX, use 2-by-2
232 * pivot block
233 *
234 KP = IMAX
235 KSTEP = 2
236 END IF
237 END IF
238 *
239 KK = K - KSTEP + 1
240 IF( KSTEP.EQ.2 )
241 $ KNC = KNC - K + 1
242 IF( KP.NE.KK ) THEN
243 *
244 * Interchange rows and columns KK and KP in the leading
245 * submatrix A(1:k,1:k)
246 *
247 CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
248 KX = KPC + KP - 1
249 DO 30 J = KP + 1, KK - 1
250 KX = KX + J - 1
251 T = AP( KNC+J-1 )
252 AP( KNC+J-1 ) = AP( KX )
253 AP( KX ) = T
254 30 CONTINUE
255 T = AP( KNC+KK-1 )
256 AP( KNC+KK-1 ) = AP( KPC+KP-1 )
257 AP( KPC+KP-1 ) = T
258 IF( KSTEP.EQ.2 ) THEN
259 T = AP( KC+K-2 )
260 AP( KC+K-2 ) = AP( KC+KP-1 )
261 AP( KC+KP-1 ) = T
262 END IF
263 END IF
264 *
265 * Update the leading submatrix
266 *
267 IF( KSTEP.EQ.1 ) THEN
268 *
269 * 1-by-1 pivot block D(k): column k now holds
270 *
271 * W(k) = U(k)*D(k)
272 *
273 * where U(k) is the k-th column of U
274 *
275 * Perform a rank-1 update of A(1:k-1,1:k-1) as
276 *
277 * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
278 *
279 R1 = ONE / AP( KC+K-1 )
280 CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
281 *
282 * Store U(k) in column k
283 *
284 CALL DSCAL( K-1, R1, AP( KC ), 1 )
285 ELSE
286 *
287 * 2-by-2 pivot block D(k): columns k and k-1 now hold
288 *
289 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
290 *
291 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
292 * of U
293 *
294 * Perform a rank-2 update of A(1:k-2,1:k-2) as
295 *
296 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
297 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
298 *
299 IF( K.GT.2 ) THEN
300 *
301 D12 = AP( K-1+( K-1 )*K / 2 )
302 D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
303 D11 = AP( K+( K-1 )*K / 2 ) / D12
304 T = ONE / ( D11*D22-ONE )
305 D12 = T / D12
306 *
307 DO 50 J = K - 2, 1, -1
308 WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
309 $ AP( J+( K-1 )*K / 2 ) )
310 WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
311 $ AP( J+( K-2 )*( K-1 ) / 2 ) )
312 DO 40 I = J, 1, -1
313 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
314 $ AP( I+( K-1 )*K / 2 )*WK -
315 $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1
316 40 CONTINUE
317 AP( J+( K-1 )*K / 2 ) = WK
318 AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
319 50 CONTINUE
320 *
321 END IF
322 *
323 END IF
324 END IF
325 *
326 * Store details of the interchanges in IPIV
327 *
328 IF( KSTEP.EQ.1 ) THEN
329 IPIV( K ) = KP
330 ELSE
331 IPIV( K ) = -KP
332 IPIV( K-1 ) = -KP
333 END IF
334 *
335 * Decrease K and return to the start of the main loop
336 *
337 K = K - KSTEP
338 KC = KNC - K
339 GO TO 10
340 *
341 ELSE
342 *
343 * Factorize A as L*D*L**T using the lower triangle of A
344 *
345 * K is the main loop index, increasing from 1 to N in steps of
346 * 1 or 2
347 *
348 K = 1
349 KC = 1
350 NPP = N*( N+1 ) / 2
351 60 CONTINUE
352 KNC = KC
353 *
354 * If K > N, exit from loop
355 *
356 IF( K.GT.N )
357 $ GO TO 110
358 KSTEP = 1
359 *
360 * Determine rows and columns to be interchanged and whether
361 * a 1-by-1 or 2-by-2 pivot block will be used
362 *
363 ABSAKK = ABS( AP( KC ) )
364 *
365 * IMAX is the row-index of the largest off-diagonal element in
366 * column K, and COLMAX is its absolute value
367 *
368 IF( K.LT.N ) THEN
369 IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
370 COLMAX = ABS( AP( KC+IMAX-K ) )
371 ELSE
372 COLMAX = ZERO
373 END IF
374 *
375 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
376 *
377 * Column K is zero: set INFO and continue
378 *
379 IF( INFO.EQ.0 )
380 $ INFO = K
381 KP = K
382 ELSE
383 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
384 *
385 * no interchange, use 1-by-1 pivot block
386 *
387 KP = K
388 ELSE
389 *
390 * JMAX is the column-index of the largest off-diagonal
391 * element in row IMAX, and ROWMAX is its absolute value
392 *
393 ROWMAX = ZERO
394 KX = KC + IMAX - K
395 DO 70 J = K, IMAX - 1
396 IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
397 ROWMAX = ABS( AP( KX ) )
398 JMAX = J
399 END IF
400 KX = KX + N - J
401 70 CONTINUE
402 KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
403 IF( IMAX.LT.N ) THEN
404 JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
405 ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
406 END IF
407 *
408 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
409 *
410 * no interchange, use 1-by-1 pivot block
411 *
412 KP = K
413 ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
414 *
415 * interchange rows and columns K and IMAX, use 1-by-1
416 * pivot block
417 *
418 KP = IMAX
419 ELSE
420 *
421 * interchange rows and columns K+1 and IMAX, use 2-by-2
422 * pivot block
423 *
424 KP = IMAX
425 KSTEP = 2
426 END IF
427 END IF
428 *
429 KK = K + KSTEP - 1
430 IF( KSTEP.EQ.2 )
431 $ KNC = KNC + N - K + 1
432 IF( KP.NE.KK ) THEN
433 *
434 * Interchange rows and columns KK and KP in the trailing
435 * submatrix A(k:n,k:n)
436 *
437 IF( KP.LT.N )
438 $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
439 $ 1 )
440 KX = KNC + KP - KK
441 DO 80 J = KK + 1, KP - 1
442 KX = KX + N - J + 1
443 T = AP( KNC+J-KK )
444 AP( KNC+J-KK ) = AP( KX )
445 AP( KX ) = T
446 80 CONTINUE
447 T = AP( KNC )
448 AP( KNC ) = AP( KPC )
449 AP( KPC ) = T
450 IF( KSTEP.EQ.2 ) THEN
451 T = AP( KC+1 )
452 AP( KC+1 ) = AP( KC+KP-K )
453 AP( KC+KP-K ) = T
454 END IF
455 END IF
456 *
457 * Update the trailing submatrix
458 *
459 IF( KSTEP.EQ.1 ) THEN
460 *
461 * 1-by-1 pivot block D(k): column k now holds
462 *
463 * W(k) = L(k)*D(k)
464 *
465 * where L(k) is the k-th column of L
466 *
467 IF( K.LT.N ) THEN
468 *
469 * Perform a rank-1 update of A(k+1:n,k+1:n) as
470 *
471 * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
472 *
473 R1 = ONE / AP( KC )
474 CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
475 $ AP( KC+N-K+1 ) )
476 *
477 * Store L(k) in column K
478 *
479 CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
480 END IF
481 ELSE
482 *
483 * 2-by-2 pivot block D(k): columns K and K+1 now hold
484 *
485 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
486 *
487 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
488 * of L
489 *
490 IF( K.LT.N-1 ) THEN
491 *
492 * Perform a rank-2 update of A(k+2:n,k+2:n) as
493 *
494 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
495 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
496 *
497 * where L(k) and L(k+1) are the k-th and (k+1)-th
498 * columns of L
499 *
500 D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 )
501 D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
502 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
503 T = ONE / ( D11*D22-ONE )
504 D21 = T / D21
505 *
506 DO 100 J = K + 2, N
507 WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
508 $ AP( J+K*( 2*N-K-1 ) / 2 ) )
509 WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
510 $ AP( J+( K-1 )*( 2*N-K ) / 2 ) )
511 *
512 DO 90 I = J, N
513 AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
514 $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
515 $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1
516 90 CONTINUE
517 *
518 AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
519 AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
520 *
521 100 CONTINUE
522 END IF
523 END IF
524 END IF
525 *
526 * Store details of the interchanges in IPIV
527 *
528 IF( KSTEP.EQ.1 ) THEN
529 IPIV( K ) = KP
530 ELSE
531 IPIV( K ) = -KP
532 IPIV( K+1 ) = -KP
533 END IF
534 *
535 * Increase K and return to the start of the main loop
536 *
537 K = K + KSTEP
538 KC = KNC + N - K + 2
539 GO TO 60
540 *
541 END IF
542 *
543 110 CONTINUE
544 RETURN
545 *
546 * End of DSPTRF
547 *
548 END