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