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