![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytf2_rk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zsytf2_rk.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: *> \brief \b ZSYTF2_RK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 unblocked algorithm).
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZSYTF2_RK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2_rk.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2_rk.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2_rk.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), E ( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *> ZSYTF2_RK computes the factorization of a complex symmetric matrix A
38: *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
39: *>
40: *> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
41: *>
42: *> where U (or L) is unit upper (or lower) triangular matrix,
43: *> U**T (or L**T) is the transpose of U (or L), P is a permutation
44: *> matrix, P**T is the transpose of P, and D is symmetric and block
45: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46: *>
47: *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48: *> For more information see Further Details section.
49: *> \endverbatim
50: *
51: * Arguments:
52: * ==========
53: *
54: *> \param[in] UPLO
55: *> \verbatim
56: *> UPLO is CHARACTER*1
57: *> Specifies whether the upper or lower triangular part of the
58: *> symmetric matrix A is stored:
59: *> = 'U': Upper triangular
60: *> = 'L': Lower triangular
61: *> \endverbatim
62: *>
63: *> \param[in] N
64: *> \verbatim
65: *> N is INTEGER
66: *> The order of the matrix A. N >= 0.
67: *> \endverbatim
68: *>
69: *> \param[in,out] A
70: *> \verbatim
71: *> A is COMPLEX*16 array, dimension (LDA,N)
72: *> On entry, the symmetric matrix A.
73: *> If UPLO = 'U': the leading N-by-N upper triangular part
74: *> of A contains the upper triangular part of the matrix A,
75: *> and the strictly lower triangular part of A is not
76: *> referenced.
77: *>
78: *> If UPLO = 'L': the leading N-by-N lower triangular part
79: *> of A contains the lower triangular part of the matrix A,
80: *> and the strictly upper triangular part of A is not
81: *> referenced.
82: *>
83: *> On exit, contains:
84: *> a) ONLY diagonal elements of the symmetric block diagonal
85: *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
86: *> (superdiagonal (or subdiagonal) elements of D
87: *> are stored on exit in array E), and
88: *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
89: *> If UPLO = 'L': factor L in the subdiagonal part of A.
90: *> \endverbatim
91: *>
92: *> \param[in] LDA
93: *> \verbatim
94: *> LDA is INTEGER
95: *> The leading dimension of the array A. LDA >= max(1,N).
96: *> \endverbatim
97: *>
98: *> \param[out] E
99: *> \verbatim
100: *> E is COMPLEX*16 array, dimension (N)
101: *> On exit, contains the superdiagonal (or subdiagonal)
102: *> elements of the symmetric block diagonal matrix D
103: *> with 1-by-1 or 2-by-2 diagonal blocks, where
104: *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
105: *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
106: *>
107: *> NOTE: For 1-by-1 diagonal block D(k), where
108: *> 1 <= k <= N, the element E(k) is set to 0 in both
109: *> UPLO = 'U' or UPLO = 'L' cases.
110: *> \endverbatim
111: *>
112: *> \param[out] IPIV
113: *> \verbatim
114: *> IPIV is INTEGER array, dimension (N)
115: *> IPIV describes the permutation matrix P in the factorization
116: *> of matrix A as follows. The absolute value of IPIV(k)
117: *> represents the index of row and column that were
118: *> interchanged with the k-th row and column. The value of UPLO
119: *> describes the order in which the interchanges were applied.
120: *> Also, the sign of IPIV represents the block structure of
121: *> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
122: *> diagonal blocks which correspond to 1 or 2 interchanges
123: *> at each factorization step. For more info see Further
124: *> Details section.
125: *>
126: *> If UPLO = 'U',
127: *> ( in factorization order, k decreases from N to 1 ):
128: *> a) A single positive entry IPIV(k) > 0 means:
129: *> D(k,k) is a 1-by-1 diagonal block.
130: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
131: *> interchanged in the matrix A(1:N,1:N);
132: *> If IPIV(k) = k, no interchange occurred.
133: *>
134: *> b) A pair of consecutive negative entries
135: *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
136: *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
137: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
138: *> 1) If -IPIV(k) != k, rows and columns
139: *> k and -IPIV(k) were interchanged
140: *> in the matrix A(1:N,1:N).
141: *> If -IPIV(k) = k, no interchange occurred.
142: *> 2) If -IPIV(k-1) != k-1, rows and columns
143: *> k-1 and -IPIV(k-1) were interchanged
144: *> in the matrix A(1:N,1:N).
145: *> If -IPIV(k-1) = k-1, no interchange occurred.
146: *>
147: *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
148: *>
149: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
150: *>
151: *> If UPLO = 'L',
152: *> ( in factorization order, k increases from 1 to N ):
153: *> a) A single positive entry IPIV(k) > 0 means:
154: *> D(k,k) is a 1-by-1 diagonal block.
155: *> If IPIV(k) != k, rows and columns k and IPIV(k) were
156: *> interchanged in the matrix A(1:N,1:N).
157: *> If IPIV(k) = k, no interchange occurred.
158: *>
159: *> b) A pair of consecutive negative entries
160: *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
161: *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
162: *> (NOTE: negative entries in IPIV appear ONLY in pairs).
163: *> 1) If -IPIV(k) != k, rows and columns
164: *> k and -IPIV(k) were interchanged
165: *> in the matrix A(1:N,1:N).
166: *> If -IPIV(k) = k, no interchange occurred.
167: *> 2) If -IPIV(k+1) != k+1, rows and columns
168: *> k-1 and -IPIV(k-1) were interchanged
169: *> in the matrix A(1:N,1:N).
170: *> If -IPIV(k+1) = k+1, no interchange occurred.
171: *>
172: *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
173: *>
174: *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
175: *> \endverbatim
176: *>
177: *> \param[out] INFO
178: *> \verbatim
179: *> INFO is INTEGER
180: *> = 0: successful exit
181: *>
182: *> < 0: If INFO = -k, the k-th argument had an illegal value
183: *>
184: *> > 0: If INFO = k, the matrix A is singular, because:
185: *> If UPLO = 'U': column k in the upper
186: *> triangular part of A contains all zeros.
187: *> If UPLO = 'L': column k in the lower
188: *> triangular part of A contains all zeros.
189: *>
190: *> Therefore D(k,k) is exactly zero, and superdiagonal
191: *> elements of column k of U (or subdiagonal elements of
192: *> column k of L ) are all zeros. The factorization has
193: *> been completed, but the block diagonal matrix D is
194: *> exactly singular, and division by zero will occur if
195: *> it is used to solve a system of equations.
196: *>
197: *> NOTE: INFO only stores the first occurrence of
198: *> a singularity, any subsequent occurrence of singularity
199: *> is not stored in INFO even though the factorization
200: *> always completes.
201: *> \endverbatim
202: *
203: * Authors:
204: * ========
205: *
206: *> \author Univ. of Tennessee
207: *> \author Univ. of California Berkeley
208: *> \author Univ. of Colorado Denver
209: *> \author NAG Ltd.
210: *
211: *> \date December 2016
212: *
213: *> \ingroup complex16SYcomputational
214: *
215: *> \par Further Details:
216: * =====================
217: *>
218: *> \verbatim
219: *> TODO: put further details
220: *> \endverbatim
221: *
222: *> \par Contributors:
223: * ==================
224: *>
225: *> \verbatim
226: *>
227: *> December 2016, Igor Kozachenko,
228: *> Computer Science Division,
229: *> University of California, Berkeley
230: *>
231: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
232: *> School of Mathematics,
233: *> University of Manchester
234: *>
235: *> 01-01-96 - Based on modifications by
236: *> J. Lewis, Boeing Computer Services Company
237: *> A. Petitet, Computer Science Dept.,
238: *> Univ. of Tenn., Knoxville abd , USA
239: *> \endverbatim
240: *
241: * =====================================================================
242: SUBROUTINE ZSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
243: *
244: * -- LAPACK computational routine (version 3.7.0) --
245: * -- LAPACK is a software package provided by Univ. of Tennessee, --
246: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
247: * December 2016
248: *
249: * .. Scalar Arguments ..
250: CHARACTER UPLO
251: INTEGER INFO, LDA, N
252: * ..
253: * .. Array Arguments ..
254: INTEGER IPIV( * )
255: COMPLEX*16 A( LDA, * ), E( * )
256: * ..
257: *
258: * =====================================================================
259: *
260: * .. Parameters ..
261: DOUBLE PRECISION ZERO, ONE
262: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
263: DOUBLE PRECISION EIGHT, SEVTEN
264: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
265: COMPLEX*16 CONE, CZERO
266: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
267: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
268: * ..
269: * .. Local Scalars ..
270: LOGICAL UPPER, DONE
271: INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
272: $ P, II
273: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
274: COMPLEX*16 D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
275: * ..
276: * .. External Functions ..
277: LOGICAL LSAME
278: INTEGER IZAMAX
279: DOUBLE PRECISION DLAMCH
280: EXTERNAL LSAME, IZAMAX, DLAMCH
281: * ..
282: * .. External Subroutines ..
283: EXTERNAL ZSCAL, ZSWAP, ZSYR, XERBLA
284: * ..
285: * .. Intrinsic Functions ..
286: INTRINSIC ABS, MAX, SQRT, DIMAG, DBLE
287: * ..
288: * .. Statement Functions ..
289: DOUBLE PRECISION CABS1
290: * ..
291: * .. Statement Function definitions ..
292: CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
293: * ..
294: * .. Executable Statements ..
295: *
296: * Test the input parameters.
297: *
298: INFO = 0
299: UPPER = LSAME( UPLO, 'U' )
300: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
301: INFO = -1
302: ELSE IF( N.LT.0 ) THEN
303: INFO = -2
304: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
305: INFO = -4
306: END IF
307: IF( INFO.NE.0 ) THEN
308: CALL XERBLA( 'ZSYTF2_RK', -INFO )
309: RETURN
310: END IF
311: *
312: * Initialize ALPHA for use in choosing pivot block size.
313: *
314: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
315: *
316: * Compute machine safe minimum
317: *
318: SFMIN = DLAMCH( 'S' )
319: *
320: IF( UPPER ) THEN
321: *
322: * Factorize A as U*D*U**T using the upper triangle of A
323: *
324: * Initilize the first entry of array E, where superdiagonal
325: * elements of D are stored
326: *
327: E( 1 ) = CZERO
328: *
329: * K is the main loop index, decreasing from N to 1 in steps of
330: * 1 or 2
331: *
332: K = N
333: 10 CONTINUE
334: *
335: * If K < 1, exit from loop
336: *
337: IF( K.LT.1 )
338: $ GO TO 34
339: KSTEP = 1
340: P = K
341: *
342: * Determine rows and columns to be interchanged and whether
343: * a 1-by-1 or 2-by-2 pivot block will be used
344: *
345: ABSAKK = CABS1( A( K, K ) )
346: *
347: * IMAX is the row-index of the largest off-diagonal element in
348: * column K, and COLMAX is its absolute value.
349: * Determine both COLMAX and IMAX.
350: *
351: IF( K.GT.1 ) THEN
352: IMAX = IZAMAX( K-1, A( 1, K ), 1 )
353: COLMAX = CABS1( A( IMAX, K ) )
354: ELSE
355: COLMAX = ZERO
356: END IF
357: *
358: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN
359: *
360: * Column K is zero or underflow: set INFO and continue
361: *
362: IF( INFO.EQ.0 )
363: $ INFO = K
364: KP = K
365: *
366: * Set E( K ) to zero
367: *
368: IF( K.GT.1 )
369: $ E( K ) = CZERO
370: *
371: ELSE
372: *
373: * Test for interchange
374: *
375: * Equivalent to testing for (used to handle NaN and Inf)
376: * ABSAKK.GE.ALPHA*COLMAX
377: *
378: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
379: *
380: * no interchange,
381: * use 1-by-1 pivot block
382: *
383: KP = K
384: ELSE
385: *
386: DONE = .FALSE.
387: *
388: * Loop until pivot found
389: *
390: 12 CONTINUE
391: *
392: * Begin pivot search loop body
393: *
394: * JMAX is the column-index of the largest off-diagonal
395: * element in row IMAX, and ROWMAX is its absolute value.
396: * Determine both ROWMAX and JMAX.
397: *
398: IF( IMAX.NE.K ) THEN
399: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ),
400: $ LDA )
401: ROWMAX = CABS1( A( IMAX, JMAX ) )
402: ELSE
403: ROWMAX = ZERO
404: END IF
405: *
406: IF( IMAX.GT.1 ) THEN
407: ITEMP = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
408: DTEMP = CABS1( A( ITEMP, IMAX ) )
409: IF( DTEMP.GT.ROWMAX ) THEN
410: ROWMAX = DTEMP
411: JMAX = ITEMP
412: END IF
413: END IF
414: *
415: * Equivalent to testing for (used to handle NaN and Inf)
416: * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
417: *
418: IF( .NOT.( CABS1( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ))
419: $ THEN
420: *
421: * interchange rows and columns K and IMAX,
422: * use 1-by-1 pivot block
423: *
424: KP = IMAX
425: DONE = .TRUE.
426: *
427: * Equivalent to testing for ROWMAX .EQ. COLMAX,
428: * used to handle NaN and Inf
429: *
430: ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
431: *
432: * interchange rows and columns K+1 and IMAX,
433: * use 2-by-2 pivot block
434: *
435: KP = IMAX
436: KSTEP = 2
437: DONE = .TRUE.
438: ELSE
439: *
440: * Pivot NOT found, set variables and repeat
441: *
442: P = IMAX
443: COLMAX = ROWMAX
444: IMAX = JMAX
445: END IF
446: *
447: * End pivot search loop body
448: *
449: IF( .NOT. DONE ) GOTO 12
450: *
451: END IF
452: *
453: * Swap TWO rows and TWO columns
454: *
455: * First swap
456: *
457: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
458: *
459: * Interchange rows and column K and P in the leading
460: * submatrix A(1:k,1:k) if we have a 2-by-2 pivot
461: *
462: IF( P.GT.1 )
463: $ CALL ZSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
464: IF( P.LT.(K-1) )
465: $ CALL ZSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ),
466: $ LDA )
467: T = A( K, K )
468: A( K, K ) = A( P, P )
469: A( P, P ) = T
470: *
471: * Convert upper triangle of A into U form by applying
472: * the interchanges in columns k+1:N.
473: *
474: IF( K.LT.N )
475: $ CALL ZSWAP( N-K, A( K, K+1 ), LDA, A( P, K+1 ), LDA )
476: *
477: END IF
478: *
479: * Second swap
480: *
481: KK = K - KSTEP + 1
482: IF( KP.NE.KK ) THEN
483: *
484: * Interchange rows and columns KK and KP in the leading
485: * submatrix A(1:k,1:k)
486: *
487: IF( KP.GT.1 )
488: $ CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
489: IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) )
490: $ CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
491: $ LDA )
492: T = A( KK, KK )
493: A( KK, KK ) = A( KP, KP )
494: A( KP, KP ) = T
495: IF( KSTEP.EQ.2 ) THEN
496: T = A( K-1, K )
497: A( K-1, K ) = A( KP, K )
498: A( KP, K ) = T
499: END IF
500: *
501: * Convert upper triangle of A into U form by applying
502: * the interchanges in columns k+1:N.
503: *
504: IF( K.LT.N )
505: $ CALL ZSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
506: $ LDA )
507: *
508: END IF
509: *
510: * Update the leading submatrix
511: *
512: IF( KSTEP.EQ.1 ) THEN
513: *
514: * 1-by-1 pivot block D(k): column k now holds
515: *
516: * W(k) = U(k)*D(k)
517: *
518: * where U(k) is the k-th column of U
519: *
520: IF( K.GT.1 ) THEN
521: *
522: * Perform a rank-1 update of A(1:k-1,1:k-1) and
523: * store U(k) in column k
524: *
525: IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
526: *
527: * Perform a rank-1 update of A(1:k-1,1:k-1) as
528: * A := A - U(k)*D(k)*U(k)**T
529: * = A - W(k)*1/D(k)*W(k)**T
530: *
531: D11 = CONE / A( K, K )
532: CALL ZSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
533: *
534: * Store U(k) in column k
535: *
536: CALL ZSCAL( K-1, D11, A( 1, K ), 1 )
537: ELSE
538: *
539: * Store L(k) in column K
540: *
541: D11 = A( K, K )
542: DO 16 II = 1, K - 1
543: A( II, K ) = A( II, K ) / D11
544: 16 CONTINUE
545: *
546: * Perform a rank-1 update of A(k+1:n,k+1:n) as
547: * A := A - U(k)*D(k)*U(k)**T
548: * = A - W(k)*(1/D(k))*W(k)**T
549: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
550: *
551: CALL ZSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
552: END IF
553: *
554: * Store the superdiagonal element of D in array E
555: *
556: E( K ) = CZERO
557: *
558: END IF
559: *
560: ELSE
561: *
562: * 2-by-2 pivot block D(k): columns k and k-1 now hold
563: *
564: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
565: *
566: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
567: * of U
568: *
569: * Perform a rank-2 update of A(1:k-2,1:k-2) as
570: *
571: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
572: * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
573: *
574: * and store L(k) and L(k+1) in columns k and k+1
575: *
576: IF( K.GT.2 ) THEN
577: *
578: D12 = A( K-1, K )
579: D22 = A( K-1, K-1 ) / D12
580: D11 = A( K, K ) / D12
581: T = CONE / ( D11*D22-CONE )
582: *
583: DO 30 J = K - 2, 1, -1
584: *
585: WKM1 = T*( D11*A( J, K-1 )-A( J, K ) )
586: WK = T*( D22*A( J, K )-A( J, K-1 ) )
587: *
588: DO 20 I = J, 1, -1
589: A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK -
590: $ ( A( I, K-1 ) / D12 )*WKM1
591: 20 CONTINUE
592: *
593: * Store U(k) and U(k-1) in cols k and k-1 for row J
594: *
595: A( J, K ) = WK / D12
596: A( J, K-1 ) = WKM1 / D12
597: *
598: 30 CONTINUE
599: *
600: END IF
601: *
602: * Copy superdiagonal elements of D(K) to E(K) and
603: * ZERO out superdiagonal entry of A
604: *
605: E( K ) = A( K-1, K )
606: E( K-1 ) = CZERO
607: A( K-1, K ) = CZERO
608: *
609: END IF
610: *
611: * End column K is nonsingular
612: *
613: END IF
614: *
615: * Store details of the interchanges in IPIV
616: *
617: IF( KSTEP.EQ.1 ) THEN
618: IPIV( K ) = KP
619: ELSE
620: IPIV( K ) = -P
621: IPIV( K-1 ) = -KP
622: END IF
623: *
624: * Decrease K and return to the start of the main loop
625: *
626: K = K - KSTEP
627: GO TO 10
628: *
629: 34 CONTINUE
630: *
631: ELSE
632: *
633: * Factorize A as L*D*L**T using the lower triangle of A
634: *
635: * Initilize the unused last entry of the subdiagonal array E.
636: *
637: E( N ) = CZERO
638: *
639: * K is the main loop index, increasing from 1 to N in steps of
640: * 1 or 2
641: *
642: K = 1
643: 40 CONTINUE
644: *
645: * If K > N, exit from loop
646: *
647: IF( K.GT.N )
648: $ GO TO 64
649: KSTEP = 1
650: P = K
651: *
652: * Determine rows and columns to be interchanged and whether
653: * a 1-by-1 or 2-by-2 pivot block will be used
654: *
655: ABSAKK = CABS1( A( K, K ) )
656: *
657: * IMAX is the row-index of the largest off-diagonal element in
658: * column K, and COLMAX is its absolute value.
659: * Determine both COLMAX and IMAX.
660: *
661: IF( K.LT.N ) THEN
662: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
663: COLMAX = CABS1( A( IMAX, K ) )
664: ELSE
665: COLMAX = ZERO
666: END IF
667: *
668: IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN
669: *
670: * Column K is zero or underflow: set INFO and continue
671: *
672: IF( INFO.EQ.0 )
673: $ INFO = K
674: KP = K
675: *
676: * Set E( K ) to zero
677: *
678: IF( K.LT.N )
679: $ E( K ) = CZERO
680: *
681: ELSE
682: *
683: * Test for interchange
684: *
685: * Equivalent to testing for (used to handle NaN and Inf)
686: * ABSAKK.GE.ALPHA*COLMAX
687: *
688: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
689: *
690: * no interchange, use 1-by-1 pivot block
691: *
692: KP = K
693: *
694: ELSE
695: *
696: DONE = .FALSE.
697: *
698: * Loop until pivot found
699: *
700: 42 CONTINUE
701: *
702: * Begin pivot search loop body
703: *
704: * JMAX is the column-index of the largest off-diagonal
705: * element in row IMAX, and ROWMAX is its absolute value.
706: * Determine both ROWMAX and JMAX.
707: *
708: IF( IMAX.NE.K ) THEN
709: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
710: ROWMAX = CABS1( A( IMAX, JMAX ) )
711: ELSE
712: ROWMAX = ZERO
713: END IF
714: *
715: IF( IMAX.LT.N ) THEN
716: ITEMP = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ),
717: $ 1 )
718: DTEMP = CABS1( A( ITEMP, IMAX ) )
719: IF( DTEMP.GT.ROWMAX ) THEN
720: ROWMAX = DTEMP
721: JMAX = ITEMP
722: END IF
723: END IF
724: *
725: * Equivalent to testing for (used to handle NaN and Inf)
726: * ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
727: *
728: IF( .NOT.( CABS1( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ))
729: $ THEN
730: *
731: * interchange rows and columns K and IMAX,
732: * use 1-by-1 pivot block
733: *
734: KP = IMAX
735: DONE = .TRUE.
736: *
737: * Equivalent to testing for ROWMAX .EQ. COLMAX,
738: * used to handle NaN and Inf
739: *
740: ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
741: *
742: * interchange rows and columns K+1 and IMAX,
743: * use 2-by-2 pivot block
744: *
745: KP = IMAX
746: KSTEP = 2
747: DONE = .TRUE.
748: ELSE
749: *
750: * Pivot NOT found, set variables and repeat
751: *
752: P = IMAX
753: COLMAX = ROWMAX
754: IMAX = JMAX
755: END IF
756: *
757: * End pivot search loop body
758: *
759: IF( .NOT. DONE ) GOTO 42
760: *
761: END IF
762: *
763: * Swap TWO rows and TWO columns
764: *
765: * First swap
766: *
767: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
768: *
769: * Interchange rows and column K and P in the trailing
770: * submatrix A(k:n,k:n) if we have a 2-by-2 pivot
771: *
772: IF( P.LT.N )
773: $ CALL ZSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
774: IF( P.GT.(K+1) )
775: $ CALL ZSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA )
776: T = A( K, K )
777: A( K, K ) = A( P, P )
778: A( P, P ) = T
779: *
780: * Convert lower triangle of A into L form by applying
781: * the interchanges in columns 1:k-1.
782: *
783: IF ( K.GT.1 )
784: $ CALL ZSWAP( K-1, A( K, 1 ), LDA, A( P, 1 ), LDA )
785: *
786: END IF
787: *
788: * Second swap
789: *
790: KK = K + KSTEP - 1
791: IF( KP.NE.KK ) THEN
792: *
793: * Interchange rows and columns KK and KP in the trailing
794: * submatrix A(k:n,k:n)
795: *
796: IF( KP.LT.N )
797: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
798: IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) )
799: $ CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
800: $ LDA )
801: T = A( KK, KK )
802: A( KK, KK ) = A( KP, KP )
803: A( KP, KP ) = T
804: IF( KSTEP.EQ.2 ) THEN
805: T = A( K+1, K )
806: A( K+1, K ) = A( KP, K )
807: A( KP, K ) = T
808: END IF
809: *
810: * Convert lower triangle of A into L form by applying
811: * the interchanges in columns 1:k-1.
812: *
813: IF ( K.GT.1 )
814: $ CALL ZSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
815: *
816: END IF
817: *
818: * Update the trailing submatrix
819: *
820: IF( KSTEP.EQ.1 ) THEN
821: *
822: * 1-by-1 pivot block D(k): column k now holds
823: *
824: * W(k) = L(k)*D(k)
825: *
826: * where L(k) is the k-th column of L
827: *
828: IF( K.LT.N ) THEN
829: *
830: * Perform a rank-1 update of A(k+1:n,k+1:n) and
831: * store L(k) in column k
832: *
833: IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
834: *
835: * Perform a rank-1 update of A(k+1:n,k+1:n) as
836: * A := A - L(k)*D(k)*L(k)**T
837: * = A - W(k)*(1/D(k))*W(k)**T
838: *
839: D11 = CONE / A( K, K )
840: CALL ZSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
841: $ A( K+1, K+1 ), LDA )
842: *
843: * Store L(k) in column k
844: *
845: CALL ZSCAL( N-K, D11, A( K+1, K ), 1 )
846: ELSE
847: *
848: * Store L(k) in column k
849: *
850: D11 = A( K, K )
851: DO 46 II = K + 1, N
852: A( II, K ) = A( II, K ) / D11
853: 46 CONTINUE
854: *
855: * Perform a rank-1 update of A(k+1:n,k+1:n) as
856: * A := A - L(k)*D(k)*L(k)**T
857: * = A - W(k)*(1/D(k))*W(k)**T
858: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
859: *
860: CALL ZSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
861: $ A( K+1, K+1 ), LDA )
862: END IF
863: *
864: * Store the subdiagonal element of D in array E
865: *
866: E( K ) = CZERO
867: *
868: END IF
869: *
870: ELSE
871: *
872: * 2-by-2 pivot block D(k): columns k and k+1 now hold
873: *
874: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
875: *
876: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
877: * of L
878: *
879: *
880: * Perform a rank-2 update of A(k+2:n,k+2:n) as
881: *
882: * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
883: * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
884: *
885: * and store L(k) and L(k+1) in columns k and k+1
886: *
887: IF( K.LT.N-1 ) THEN
888: *
889: D21 = A( K+1, K )
890: D11 = A( K+1, K+1 ) / D21
891: D22 = A( K, K ) / D21
892: T = CONE / ( D11*D22-CONE )
893: *
894: DO 60 J = K + 2, N
895: *
896: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
897: *
898: WK = T*( D11*A( J, K )-A( J, K+1 ) )
899: WKP1 = T*( D22*A( J, K+1 )-A( J, K ) )
900: *
901: * Perform a rank-2 update of A(k+2:n,k+2:n)
902: *
903: DO 50 I = J, N
904: A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK -
905: $ ( A( I, K+1 ) / D21 )*WKP1
906: 50 CONTINUE
907: *
908: * Store L(k) and L(k+1) in cols k and k+1 for row J
909: *
910: A( J, K ) = WK / D21
911: A( J, K+1 ) = WKP1 / D21
912: *
913: 60 CONTINUE
914: *
915: END IF
916: *
917: * Copy subdiagonal elements of D(K) to E(K) and
918: * ZERO out subdiagonal entry of A
919: *
920: E( K ) = A( K+1, K )
921: E( K+1 ) = CZERO
922: A( K+1, K ) = CZERO
923: *
924: END IF
925: *
926: * End column K is nonsingular
927: *
928: END IF
929: *
930: * Store details of the interchanges in IPIV
931: *
932: IF( KSTEP.EQ.1 ) THEN
933: IPIV( K ) = KP
934: ELSE
935: IPIV( K ) = -P
936: IPIV( K+1 ) = -KP
937: END IF
938: *
939: * Increase K and return to the start of the main loop
940: *
941: K = K + KSTEP
942: GO TO 40
943: *
944: 64 CONTINUE
945: *
946: END IF
947: *
948: RETURN
949: *
950: * End of ZSYTF2_RK
951: *
952: END