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