1: *> \brief \b ZHETF2_ROOK computes the factorization of a complex Hermitian 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 ZHETF2_ROOK + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2_rook.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2_rook.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2_rook.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETF2_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: *> ZHETF2_ROOK computes the factorization of a complex Hermitian matrix A
39: *> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:
40: *>
41: *> A = U*D*U**H or A = L*D*L**H
42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, U**H is the conjugate transpose of U, and D is
45: *> Hermitian and 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: *> Hermitian 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 Hermitian 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) were
96: *> 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 complex16HEcomputational
135: *
136: *> \par Further Details:
137: * =====================
138: *>
139: *> \verbatim
140: *>
141: *> If UPLO = 'U', then A = U*D*U**H, 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**H, 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, USA
192: *> \endverbatim
193: *
194: * =====================================================================
195: SUBROUTINE ZHETF2_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: * ..
219: * .. Local Scalars ..
220: LOGICAL DONE, UPPER
221: INTEGER I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP,
222: $ P
223: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, DTEMP,
224: $ ROWMAX, TT, SFMIN
225: COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, Z
226: * ..
227: * .. External Functions ..
228: *
229: LOGICAL LSAME
230: INTEGER IZAMAX
231: DOUBLE PRECISION DLAMCH, DLAPY2
232: EXTERNAL LSAME, IZAMAX, DLAMCH, DLAPY2
233: * ..
234: * .. External Subroutines ..
235: EXTERNAL XERBLA, ZDSCAL, ZHER, ZSWAP
236: * ..
237: * .. Intrinsic Functions ..
238: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
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( 'ZHETF2_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**H 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 = ABS( DBLE( 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: A( K, K ) = DBLE( A( K, K ) )
313: ELSE
314: *
315: * ============================================================
316: *
317: * BEGIN pivot search
318: *
319: * Case(1)
320: * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
321: * (used to handle NaN and Inf)
322: *
323: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
324: *
325: * no interchange, use 1-by-1 pivot block
326: *
327: KP = K
328: *
329: ELSE
330: *
331: DONE = .FALSE.
332: *
333: * Loop until pivot found
334: *
335: 12 CONTINUE
336: *
337: * BEGIN pivot search loop body
338: *
339: *
340: * JMAX is the column-index of the largest off-diagonal
341: * element in row IMAX, and ROWMAX is its absolute value.
342: * Determine both ROWMAX and JMAX.
343: *
344: IF( IMAX.NE.K ) THEN
345: JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ),
346: $ LDA )
347: ROWMAX = CABS1( A( IMAX, JMAX ) )
348: ELSE
349: ROWMAX = ZERO
350: END IF
351: *
352: IF( IMAX.GT.1 ) THEN
353: ITEMP = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
354: DTEMP = CABS1( A( ITEMP, IMAX ) )
355: IF( DTEMP.GT.ROWMAX ) THEN
356: ROWMAX = DTEMP
357: JMAX = ITEMP
358: END IF
359: END IF
360: *
361: * Case(2)
362: * Equivalent to testing for
363: * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
364: * (used to handle NaN and Inf)
365: *
366: IF( .NOT.( ABS( DBLE( A( IMAX, IMAX ) ) )
367: $ .LT.ALPHA*ROWMAX ) ) THEN
368: *
369: * interchange rows and columns K and IMAX,
370: * use 1-by-1 pivot block
371: *
372: KP = IMAX
373: DONE = .TRUE.
374: *
375: * Case(3)
376: * Equivalent to testing for ROWMAX.EQ.COLMAX,
377: * (used to handle NaN and Inf)
378: *
379: ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
380: $ THEN
381: *
382: * interchange rows and columns K-1 and IMAX,
383: * use 2-by-2 pivot block
384: *
385: KP = IMAX
386: KSTEP = 2
387: DONE = .TRUE.
388: *
389: * Case(4)
390: ELSE
391: *
392: * Pivot not found: set params and repeat
393: *
394: P = IMAX
395: COLMAX = ROWMAX
396: IMAX = JMAX
397: END IF
398: *
399: * END pivot search loop body
400: *
401: IF( .NOT.DONE ) GOTO 12
402: *
403: END IF
404: *
405: * END pivot search
406: *
407: * ============================================================
408: *
409: * KK is the column of A where pivoting step stopped
410: *
411: KK = K - KSTEP + 1
412: *
413: * For only a 2x2 pivot, interchange rows and columns K and P
414: * in the leading submatrix A(1:k,1:k)
415: *
416: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
417: * (1) Swap columnar parts
418: IF( P.GT.1 )
419: $ CALL ZSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
420: * (2) Swap and conjugate middle parts
421: DO 14 J = P + 1, K - 1
422: T = DCONJG( A( J, K ) )
423: A( J, K ) = DCONJG( A( P, J ) )
424: A( P, J ) = T
425: 14 CONTINUE
426: * (3) Swap and conjugate corner elements at row-col interserction
427: A( P, K ) = DCONJG( A( P, K ) )
428: * (4) Swap diagonal elements at row-col intersection
429: R1 = DBLE( A( K, K ) )
430: A( K, K ) = DBLE( A( P, P ) )
431: A( P, P ) = R1
432: END IF
433: *
434: * For both 1x1 and 2x2 pivots, interchange rows and
435: * columns KK and KP in the leading submatrix A(1:k,1:k)
436: *
437: IF( KP.NE.KK ) THEN
438: * (1) Swap columnar parts
439: IF( KP.GT.1 )
440: $ CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
441: * (2) Swap and conjugate middle parts
442: DO 15 J = KP + 1, KK - 1
443: T = DCONJG( A( J, KK ) )
444: A( J, KK ) = DCONJG( A( KP, J ) )
445: A( KP, J ) = T
446: 15 CONTINUE
447: * (3) Swap and conjugate corner elements at row-col interserction
448: A( KP, KK ) = DCONJG( A( KP, KK ) )
449: * (4) Swap diagonal elements at row-col intersection
450: R1 = DBLE( A( KK, KK ) )
451: A( KK, KK ) = DBLE( A( KP, KP ) )
452: A( KP, KP ) = R1
453: *
454: IF( KSTEP.EQ.2 ) THEN
455: * (*) Make sure that diagonal element of pivot is real
456: A( K, K ) = DBLE( A( K, K ) )
457: * (5) Swap row elements
458: T = A( K-1, K )
459: A( K-1, K ) = A( KP, K )
460: A( KP, K ) = T
461: END IF
462: ELSE
463: * (*) Make sure that diagonal element of pivot is real
464: A( K, K ) = DBLE( A( K, K ) )
465: IF( KSTEP.EQ.2 )
466: $ A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
467: END IF
468: *
469: * Update the leading submatrix
470: *
471: IF( KSTEP.EQ.1 ) THEN
472: *
473: * 1-by-1 pivot block D(k): column k now holds
474: *
475: * W(k) = U(k)*D(k)
476: *
477: * where U(k) is the k-th column of U
478: *
479: IF( K.GT.1 ) THEN
480: *
481: * Perform a rank-1 update of A(1:k-1,1:k-1) and
482: * store U(k) in column k
483: *
484: IF( ABS( DBLE( A( K, K ) ) ).GE.SFMIN ) THEN
485: *
486: * Perform a rank-1 update of A(1:k-1,1:k-1) as
487: * A := A - U(k)*D(k)*U(k)**T
488: * = A - W(k)*1/D(k)*W(k)**T
489: *
490: D11 = ONE / DBLE( A( K, K ) )
491: CALL ZHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
492: *
493: * Store U(k) in column k
494: *
495: CALL ZDSCAL( K-1, D11, A( 1, K ), 1 )
496: ELSE
497: *
498: * Store L(k) in column K
499: *
500: D11 = DBLE( A( K, K ) )
501: DO 16 II = 1, K - 1
502: A( II, K ) = A( II, K ) / D11
503: 16 CONTINUE
504: *
505: * Perform a rank-1 update of A(k+1:n,k+1:n) as
506: * A := A - U(k)*D(k)*U(k)**T
507: * = A - W(k)*(1/D(k))*W(k)**T
508: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
509: *
510: CALL ZHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
511: END IF
512: END IF
513: *
514: ELSE
515: *
516: * 2-by-2 pivot block D(k): columns k and k-1 now hold
517: *
518: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
519: *
520: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
521: * of U
522: *
523: * Perform a rank-2 update of A(1:k-2,1:k-2) as
524: *
525: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
526: * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
527: *
528: * and store L(k) and L(k+1) in columns k and k+1
529: *
530: IF( K.GT.2 ) THEN
531: * D = |A12|
532: D = DLAPY2( DBLE( A( K-1, K ) ),
533: $ DIMAG( A( K-1, K ) ) )
534: D11 = A( K, K ) / D
535: D22 = A( K-1, K-1 ) / D
536: D12 = A( K-1, K ) / D
537: TT = ONE / ( D11*D22-ONE )
538: *
539: DO 30 J = K - 2, 1, -1
540: *
541: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
542: *
543: WKM1 = TT*( D11*A( J, K-1 )-DCONJG( D12 )*
544: $ A( J, K ) )
545: WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) )
546: *
547: * Perform a rank-2 update of A(1:k-2,1:k-2)
548: *
549: DO 20 I = J, 1, -1
550: A( I, J ) = A( I, J ) -
551: $ ( A( I, K ) / D )*DCONJG( WK ) -
552: $ ( A( I, K-1 ) / D )*DCONJG( WKM1 )
553: 20 CONTINUE
554: *
555: * Store U(k) and U(k-1) in cols k and k-1 for row J
556: *
557: A( J, K ) = WK / D
558: A( J, K-1 ) = WKM1 / D
559: * (*) Make sure that diagonal element of pivot is real
560: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO )
561: *
562: 30 CONTINUE
563: *
564: END IF
565: *
566: END IF
567: *
568: END IF
569: *
570: * Store details of the interchanges in IPIV
571: *
572: IF( KSTEP.EQ.1 ) THEN
573: IPIV( K ) = KP
574: ELSE
575: IPIV( K ) = -P
576: IPIV( K-1 ) = -KP
577: END IF
578: *
579: * Decrease K and return to the start of the main loop
580: *
581: K = K - KSTEP
582: GO TO 10
583: *
584: ELSE
585: *
586: * Factorize A as L*D*L**H using the lower triangle of A
587: *
588: * K is the main loop index, increasing from 1 to N in steps of
589: * 1 or 2
590: *
591: K = 1
592: 40 CONTINUE
593: *
594: * If K > N, exit from loop
595: *
596: IF( K.GT.N )
597: $ GO TO 70
598: KSTEP = 1
599: P = K
600: *
601: * Determine rows and columns to be interchanged and whether
602: * a 1-by-1 or 2-by-2 pivot block will be used
603: *
604: ABSAKK = ABS( DBLE( A( K, K ) ) )
605: *
606: * IMAX is the row-index of the largest off-diagonal element in
607: * column K, and COLMAX is its absolute value.
608: * Determine both COLMAX and IMAX.
609: *
610: IF( K.LT.N ) THEN
611: IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
612: COLMAX = CABS1( A( IMAX, K ) )
613: ELSE
614: COLMAX = ZERO
615: END IF
616: *
617: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
618: *
619: * Column K is zero or underflow: set INFO and continue
620: *
621: IF( INFO.EQ.0 )
622: $ INFO = K
623: KP = K
624: A( K, K ) = DBLE( A( K, K ) )
625: ELSE
626: *
627: * ============================================================
628: *
629: * BEGIN pivot search
630: *
631: * Case(1)
632: * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
633: * (used to handle NaN and Inf)
634: *
635: IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
636: *
637: * no interchange, use 1-by-1 pivot block
638: *
639: KP = K
640: *
641: ELSE
642: *
643: DONE = .FALSE.
644: *
645: * Loop until pivot found
646: *
647: 42 CONTINUE
648: *
649: * BEGIN pivot search loop body
650: *
651: *
652: * JMAX is the column-index of the largest off-diagonal
653: * element in row IMAX, and ROWMAX is its absolute value.
654: * Determine both ROWMAX and JMAX.
655: *
656: IF( IMAX.NE.K ) THEN
657: JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
658: ROWMAX = CABS1( A( IMAX, JMAX ) )
659: ELSE
660: ROWMAX = ZERO
661: END IF
662: *
663: IF( IMAX.LT.N ) THEN
664: ITEMP = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ),
665: $ 1 )
666: DTEMP = CABS1( A( ITEMP, IMAX ) )
667: IF( DTEMP.GT.ROWMAX ) THEN
668: ROWMAX = DTEMP
669: JMAX = ITEMP
670: END IF
671: END IF
672: *
673: * Case(2)
674: * Equivalent to testing for
675: * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
676: * (used to handle NaN and Inf)
677: *
678: IF( .NOT.( ABS( DBLE( A( IMAX, IMAX ) ) )
679: $ .LT.ALPHA*ROWMAX ) ) THEN
680: *
681: * interchange rows and columns K and IMAX,
682: * use 1-by-1 pivot block
683: *
684: KP = IMAX
685: DONE = .TRUE.
686: *
687: * Case(3)
688: * Equivalent to testing for ROWMAX.EQ.COLMAX,
689: * (used to handle NaN and Inf)
690: *
691: ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
692: $ THEN
693: *
694: * interchange rows and columns K+1 and IMAX,
695: * use 2-by-2 pivot block
696: *
697: KP = IMAX
698: KSTEP = 2
699: DONE = .TRUE.
700: *
701: * Case(4)
702: ELSE
703: *
704: * Pivot not found: set params and repeat
705: *
706: P = IMAX
707: COLMAX = ROWMAX
708: IMAX = JMAX
709: END IF
710: *
711: *
712: * END pivot search loop body
713: *
714: IF( .NOT.DONE ) GOTO 42
715: *
716: END IF
717: *
718: * END pivot search
719: *
720: * ============================================================
721: *
722: * KK is the column of A where pivoting step stopped
723: *
724: KK = K + KSTEP - 1
725: *
726: * For only a 2x2 pivot, interchange rows and columns K and P
727: * in the trailing submatrix A(k:n,k:n)
728: *
729: IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
730: * (1) Swap columnar parts
731: IF( P.LT.N )
732: $ CALL ZSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
733: * (2) Swap and conjugate middle parts
734: DO 44 J = K + 1, P - 1
735: T = DCONJG( A( J, K ) )
736: A( J, K ) = DCONJG( A( P, J ) )
737: A( P, J ) = T
738: 44 CONTINUE
739: * (3) Swap and conjugate corner elements at row-col interserction
740: A( P, K ) = DCONJG( A( P, K ) )
741: * (4) Swap diagonal elements at row-col intersection
742: R1 = DBLE( A( K, K ) )
743: A( K, K ) = DBLE( A( P, P ) )
744: A( P, P ) = R1
745: END IF
746: *
747: * For both 1x1 and 2x2 pivots, interchange rows and
748: * columns KK and KP in the trailing submatrix A(k:n,k:n)
749: *
750: IF( KP.NE.KK ) THEN
751: * (1) Swap columnar parts
752: IF( KP.LT.N )
753: $ CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
754: * (2) Swap and conjugate middle parts
755: DO 45 J = KK + 1, KP - 1
756: T = DCONJG( A( J, KK ) )
757: A( J, KK ) = DCONJG( A( KP, J ) )
758: A( KP, J ) = T
759: 45 CONTINUE
760: * (3) Swap and conjugate corner elements at row-col interserction
761: A( KP, KK ) = DCONJG( A( KP, KK ) )
762: * (4) Swap diagonal elements at row-col intersection
763: R1 = DBLE( A( KK, KK ) )
764: A( KK, KK ) = DBLE( A( KP, KP ) )
765: A( KP, KP ) = R1
766: *
767: IF( KSTEP.EQ.2 ) THEN
768: * (*) Make sure that diagonal element of pivot is real
769: A( K, K ) = DBLE( A( K, K ) )
770: * (5) Swap row elements
771: T = A( K+1, K )
772: A( K+1, K ) = A( KP, K )
773: A( KP, K ) = T
774: END IF
775: ELSE
776: * (*) Make sure that diagonal element of pivot is real
777: A( K, K ) = DBLE( A( K, K ) )
778: IF( KSTEP.EQ.2 )
779: $ A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
780: END IF
781: *
782: * Update the trailing submatrix
783: *
784: IF( KSTEP.EQ.1 ) THEN
785: *
786: * 1-by-1 pivot block D(k): column k of A now holds
787: *
788: * W(k) = L(k)*D(k),
789: *
790: * where L(k) is the k-th column of L
791: *
792: IF( K.LT.N ) THEN
793: *
794: * Perform a rank-1 update of A(k+1:n,k+1:n) and
795: * store L(k) in column k
796: *
797: * Handle division by a small number
798: *
799: IF( ABS( DBLE( A( K, K ) ) ).GE.SFMIN ) THEN
800: *
801: * Perform a rank-1 update of A(k+1:n,k+1:n) as
802: * A := A - L(k)*D(k)*L(k)**T
803: * = A - W(k)*(1/D(k))*W(k)**T
804: *
805: D11 = ONE / DBLE( A( K, K ) )
806: CALL ZHER( UPLO, N-K, -D11, A( K+1, K ), 1,
807: $ A( K+1, K+1 ), LDA )
808: *
809: * Store L(k) in column k
810: *
811: CALL ZDSCAL( N-K, D11, A( K+1, K ), 1 )
812: ELSE
813: *
814: * Store L(k) in column k
815: *
816: D11 = DBLE( A( K, K ) )
817: DO 46 II = K + 1, N
818: A( II, K ) = A( II, K ) / D11
819: 46 CONTINUE
820: *
821: * Perform a rank-1 update of A(k+1:n,k+1:n) as
822: * A := A - L(k)*D(k)*L(k)**T
823: * = A - W(k)*(1/D(k))*W(k)**T
824: * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
825: *
826: CALL ZHER( UPLO, N-K, -D11, A( K+1, K ), 1,
827: $ A( K+1, K+1 ), LDA )
828: END IF
829: END IF
830: *
831: ELSE
832: *
833: * 2-by-2 pivot block D(k): columns k and k+1 now hold
834: *
835: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
836: *
837: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
838: * of L
839: *
840: *
841: * Perform a rank-2 update of A(k+2:n,k+2:n) as
842: *
843: * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
844: * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
845: *
846: * and store L(k) and L(k+1) in columns k and k+1
847: *
848: IF( K.LT.N-1 ) THEN
849: * D = |A21|
850: D = DLAPY2( DBLE( A( K+1, K ) ),
851: $ DIMAG( A( K+1, K ) ) )
852: D11 = DBLE( A( K+1, K+1 ) ) / D
853: D22 = DBLE( A( K, K ) ) / D
854: D21 = A( K+1, K ) / D
855: TT = ONE / ( D11*D22-ONE )
856: *
857: DO 60 J = K + 2, N
858: *
859: * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
860: *
861: WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) )
862: WKP1 = TT*( D22*A( J, K+1 )-DCONJG( D21 )*
863: $ A( J, K ) )
864: *
865: * Perform a rank-2 update of A(k+2:n,k+2:n)
866: *
867: DO 50 I = J, N
868: A( I, J ) = A( I, J ) -
869: $ ( A( I, K ) / D )*DCONJG( WK ) -
870: $ ( A( I, K+1 ) / D )*DCONJG( WKP1 )
871: 50 CONTINUE
872: *
873: * Store L(k) and L(k+1) in cols k and k+1 for row J
874: *
875: A( J, K ) = WK / D
876: A( J, K+1 ) = WKP1 / D
877: * (*) Make sure that diagonal element of pivot is real
878: A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO )
879: *
880: 60 CONTINUE
881: *
882: END IF
883: *
884: END IF
885: *
886: END IF
887: *
888: * Store details of the interchanges in IPIV
889: *
890: IF( KSTEP.EQ.1 ) THEN
891: IPIV( K ) = KP
892: ELSE
893: IPIV( K ) = -P
894: IPIV( K+1 ) = -KP
895: END IF
896: *
897: * Increase K and return to the start of the main loop
898: *
899: K = K + KSTEP
900: GO TO 40
901: *
902: END IF
903: *
904: 70 CONTINUE
905: *
906: RETURN
907: *
908: * End of ZHETF2_ROOK
909: *
910: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetf2_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack