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