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