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