001:       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            INFO, LDA, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IPIV( * )
014:       COMPLEX*16         A( LDA, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZHETF2 computes the factorization of a complex Hermitian matrix A
021: *  using the Bunch-Kaufman diagonal pivoting method:
022: *
023: *     A = U*D*U'  or  A = L*D*L'
024: *
025: *  where U (or L) is a product of permutation and unit upper (lower)
026: *  triangular matrices, U' is the conjugate transpose of U, and D is
027: *  Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
028: *
029: *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  UPLO    (input) CHARACTER*1
035: *          Specifies whether the upper or lower triangular part of the
036: *          Hermitian matrix A is stored:
037: *          = 'U':  Upper triangular
038: *          = 'L':  Lower triangular
039: *
040: *  N       (input) INTEGER
041: *          The order of the matrix A.  N >= 0.
042: *
043: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
044: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
045: *          n-by-n upper triangular part of A contains the upper
046: *          triangular part of the matrix A, and the strictly lower
047: *          triangular part of A is not referenced.  If UPLO = 'L', the
048: *          leading n-by-n lower triangular part of A contains the lower
049: *          triangular part of the matrix A, and the strictly upper
050: *          triangular part of A is not referenced.
051: *
052: *          On exit, the block diagonal matrix D and the multipliers used
053: *          to obtain the factor U or L (see below for further details).
054: *
055: *  LDA     (input) INTEGER
056: *          The leading dimension of the array A.  LDA >= max(1,N).
057: *
058: *  IPIV    (output) INTEGER array, dimension (N)
059: *          Details of the interchanges and the block structure of D.
060: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
061: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
062: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
063: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
064: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
065: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
066: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
067: *
068: *  INFO    (output) INTEGER
069: *          = 0: successful exit
070: *          < 0: if INFO = -k, the k-th argument had an illegal value
071: *          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
072: *               has been completed, but the block diagonal matrix D is
073: *               exactly singular, and division by zero will occur if it
074: *               is used to solve a system of equations.
075: *
076: *  Further Details
077: *  ===============
078: *
079: *  09-29-06 - patch from
080: *    Bobby Cheng, MathWorks
081: *
082: *    Replace l.210 and l.393
083: *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
084: *    by
085: *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
086: *
087: *  01-01-96 - Based on modifications by
088: *    J. Lewis, Boeing Computer Services Company
089: *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
090: *
091: *  If UPLO = 'U', then A = U*D*U', where
092: *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
093: *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
094: *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
095: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
096: *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
097: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
098: *
099: *             (   I    v    0   )   k-s
100: *     U(k) =  (   0    I    0   )   s
101: *             (   0    0    I   )   n-k
102: *                k-s   s   n-k
103: *
104: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
105: *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
106: *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
107: *
108: *  If UPLO = 'L', then A = L*D*L', where
109: *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
110: *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
111: *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
112: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
113: *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
114: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
115: *
116: *             (   I    0     0   )  k-1
117: *     L(k) =  (   0    I     0   )  s
118: *             (   0    v     I   )  n-k-s+1
119: *                k-1   s  n-k-s+1
120: *
121: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
122: *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
123: *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
124: *
125: *  =====================================================================
126: *
127: *     .. Parameters ..
128:       DOUBLE PRECISION   ZERO, ONE
129:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
130:       DOUBLE PRECISION   EIGHT, SEVTEN
131:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
132: *     ..
133: *     .. Local Scalars ..
134:       LOGICAL            UPPER
135:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
136:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
137:      $                   TT
138:       COMPLEX*16         D12, D21, T, WK, WKM1, WKP1, ZDUM
139: *     ..
140: *     .. External Functions ..
141:       LOGICAL            LSAME, DISNAN
142:       INTEGER            IZAMAX
143:       DOUBLE PRECISION   DLAPY2
144:       EXTERNAL           LSAME, IZAMAX, DLAPY2, DISNAN
145: *     ..
146: *     .. External Subroutines ..
147:       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZSWAP
148: *     ..
149: *     .. Intrinsic Functions ..
150:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
151: *     ..
152: *     .. Statement Functions ..
153:       DOUBLE PRECISION   CABS1
154: *     ..
155: *     .. Statement Function definitions ..
156:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
157: *     ..
158: *     .. Executable Statements ..
159: *
160: *     Test the input parameters.
161: *
162:       INFO = 0
163:       UPPER = LSAME( UPLO, 'U' )
164:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
165:          INFO = -1
166:       ELSE IF( N.LT.0 ) THEN
167:          INFO = -2
168:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
169:          INFO = -4
170:       END IF
171:       IF( INFO.NE.0 ) THEN
172:          CALL XERBLA( 'ZHETF2', -INFO )
173:          RETURN
174:       END IF
175: *
176: *     Initialize ALPHA for use in choosing pivot block size.
177: *
178:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
179: *
180:       IF( UPPER ) THEN
181: *
182: *        Factorize A as U*D*U' using the upper triangle of A
183: *
184: *        K is the main loop index, decreasing from N to 1 in steps of
185: *        1 or 2
186: *
187:          K = N
188:    10    CONTINUE
189: *
190: *        If K < 1, exit from loop
191: *
192:          IF( K.LT.1 )
193:      $      GO TO 90
194:          KSTEP = 1
195: *
196: *        Determine rows and columns to be interchanged and whether
197: *        a 1-by-1 or 2-by-2 pivot block will be used
198: *
199:          ABSAKK = ABS( DBLE( A( K, K ) ) )
200: *
201: *        IMAX is the row-index of the largest off-diagonal element in
202: *        column K, and COLMAX is its absolute value
203: *
204:          IF( K.GT.1 ) THEN
205:             IMAX = IZAMAX( K-1, A( 1, K ), 1 )
206:             COLMAX = CABS1( A( IMAX, K ) )
207:          ELSE
208:             COLMAX = ZERO
209:          END IF
210: *
211:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
212: *
213: *           Column K is zero or contains a NaN: set INFO and continue
214: *
215:             IF( INFO.EQ.0 )
216:      $         INFO = K
217:             KP = K
218:             A( K, K ) = DBLE( A( K, K ) )
219:          ELSE
220:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
221: *
222: *              no interchange, use 1-by-1 pivot block
223: *
224:                KP = K
225:             ELSE
226: *
227: *              JMAX is the column-index of the largest off-diagonal
228: *              element in row IMAX, and ROWMAX is its absolute value
229: *
230:                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
231:                ROWMAX = CABS1( A( IMAX, JMAX ) )
232:                IF( IMAX.GT.1 ) THEN
233:                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
234:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
235:                END IF
236: *
237:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
238: *
239: *                 no interchange, use 1-by-1 pivot block
240: *
241:                   KP = K
242:                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
243:      $                   THEN
244: *
245: *                 interchange rows and columns K and IMAX, use 1-by-1
246: *                 pivot block
247: *
248:                   KP = IMAX
249:                ELSE
250: *
251: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
252: *                 pivot block
253: *
254:                   KP = IMAX
255:                   KSTEP = 2
256:                END IF
257:             END IF
258: *
259:             KK = K - KSTEP + 1
260:             IF( KP.NE.KK ) THEN
261: *
262: *              Interchange rows and columns KK and KP in the leading
263: *              submatrix A(1:k,1:k)
264: *
265:                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
266:                DO 20 J = KP + 1, KK - 1
267:                   T = DCONJG( A( J, KK ) )
268:                   A( J, KK ) = DCONJG( A( KP, J ) )
269:                   A( KP, J ) = T
270:    20          CONTINUE
271:                A( KP, KK ) = DCONJG( A( KP, KK ) )
272:                R1 = DBLE( A( KK, KK ) )
273:                A( KK, KK ) = DBLE( A( KP, KP ) )
274:                A( KP, KP ) = R1
275:                IF( KSTEP.EQ.2 ) THEN
276:                   A( K, K ) = DBLE( A( K, K ) )
277:                   T = A( K-1, K )
278:                   A( K-1, K ) = A( KP, K )
279:                   A( KP, K ) = T
280:                END IF
281:             ELSE
282:                A( K, K ) = DBLE( A( K, K ) )
283:                IF( KSTEP.EQ.2 )
284:      $            A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
285:             END IF
286: *
287: *           Update the leading submatrix
288: *
289:             IF( KSTEP.EQ.1 ) THEN
290: *
291: *              1-by-1 pivot block D(k): column k now holds
292: *
293: *              W(k) = U(k)*D(k)
294: *
295: *              where U(k) is the k-th column of U
296: *
297: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
298: *
299: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
300: *
301:                R1 = ONE / DBLE( A( K, K ) )
302:                CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
303: *
304: *              Store U(k) in column k
305: *
306:                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
307:             ELSE
308: *
309: *              2-by-2 pivot block D(k): columns k and k-1 now hold
310: *
311: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
312: *
313: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
314: *              of U
315: *
316: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
317: *
318: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
319: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
320: *
321:                IF( K.GT.2 ) THEN
322: *
323:                   D = DLAPY2( DBLE( A( K-1, K ) ),
324:      $                DIMAG( A( K-1, K ) ) )
325:                   D22 = DBLE( A( K-1, K-1 ) ) / D
326:                   D11 = DBLE( A( K, K ) ) / D
327:                   TT = ONE / ( D11*D22-ONE )
328:                   D12 = A( K-1, K ) / D
329:                   D = TT / D
330: *
331:                   DO 40 J = K - 2, 1, -1
332:                      WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
333:      $                      A( J, K ) )
334:                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
335:                      DO 30 I = J, 1, -1
336:                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
337:      $                              A( I, K-1 )*DCONJG( WKM1 )
338:    30                CONTINUE
339:                      A( J, K ) = WK
340:                      A( J, K-1 ) = WKM1
341:                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
342:    40             CONTINUE
343: *
344:                END IF
345: *
346:             END IF
347:          END IF
348: *
349: *        Store details of the interchanges in IPIV
350: *
351:          IF( KSTEP.EQ.1 ) THEN
352:             IPIV( K ) = KP
353:          ELSE
354:             IPIV( K ) = -KP
355:             IPIV( K-1 ) = -KP
356:          END IF
357: *
358: *        Decrease K and return to the start of the main loop
359: *
360:          K = K - KSTEP
361:          GO TO 10
362: *
363:       ELSE
364: *
365: *        Factorize A as L*D*L' using the lower triangle of A
366: *
367: *        K is the main loop index, increasing from 1 to N in steps of
368: *        1 or 2
369: *
370:          K = 1
371:    50    CONTINUE
372: *
373: *        If K > N, exit from loop
374: *
375:          IF( K.GT.N )
376:      $      GO TO 90
377:          KSTEP = 1
378: *
379: *        Determine rows and columns to be interchanged and whether
380: *        a 1-by-1 or 2-by-2 pivot block will be used
381: *
382:          ABSAKK = ABS( DBLE( A( K, K ) ) )
383: *
384: *        IMAX is the row-index of the largest off-diagonal element in
385: *        column K, and COLMAX is its absolute value
386: *
387:          IF( K.LT.N ) THEN
388:             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
389:             COLMAX = CABS1( A( IMAX, K ) )
390:          ELSE
391:             COLMAX = ZERO
392:          END IF
393: *
394:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
395: *
396: *           Column K is zero or contains a NaN: set INFO and continue
397: *
398:             IF( INFO.EQ.0 )
399:      $         INFO = K
400:             KP = K
401:             A( K, K ) = DBLE( A( K, K ) )
402:          ELSE
403:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
404: *
405: *              no interchange, use 1-by-1 pivot block
406: *
407:                KP = K
408:             ELSE
409: *
410: *              JMAX is the column-index of the largest off-diagonal
411: *              element in row IMAX, and ROWMAX is its absolute value
412: *
413:                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
414:                ROWMAX = CABS1( A( IMAX, JMAX ) )
415:                IF( IMAX.LT.N ) THEN
416:                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
417:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
418:                END IF
419: *
420:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
421: *
422: *                 no interchange, use 1-by-1 pivot block
423: *
424:                   KP = K
425:                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
426:      $                   THEN
427: *
428: *                 interchange rows and columns K and IMAX, use 1-by-1
429: *                 pivot block
430: *
431:                   KP = IMAX
432:                ELSE
433: *
434: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
435: *                 pivot block
436: *
437:                   KP = IMAX
438:                   KSTEP = 2
439:                END IF
440:             END IF
441: *
442:             KK = K + KSTEP - 1
443:             IF( KP.NE.KK ) THEN
444: *
445: *              Interchange rows and columns KK and KP in the trailing
446: *              submatrix A(k:n,k:n)
447: *
448:                IF( KP.LT.N )
449:      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
450:                DO 60 J = KK + 1, KP - 1
451:                   T = DCONJG( A( J, KK ) )
452:                   A( J, KK ) = DCONJG( A( KP, J ) )
453:                   A( KP, J ) = T
454:    60          CONTINUE
455:                A( KP, KK ) = DCONJG( A( KP, KK ) )
456:                R1 = DBLE( A( KK, KK ) )
457:                A( KK, KK ) = DBLE( A( KP, KP ) )
458:                A( KP, KP ) = R1
459:                IF( KSTEP.EQ.2 ) THEN
460:                   A( K, K ) = DBLE( A( K, K ) )
461:                   T = A( K+1, K )
462:                   A( K+1, K ) = A( KP, K )
463:                   A( KP, K ) = T
464:                END IF
465:             ELSE
466:                A( K, K ) = DBLE( A( K, K ) )
467:                IF( KSTEP.EQ.2 )
468:      $            A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
469:             END IF
470: *
471: *           Update the trailing submatrix
472: *
473:             IF( KSTEP.EQ.1 ) THEN
474: *
475: *              1-by-1 pivot block D(k): column k now holds
476: *
477: *              W(k) = L(k)*D(k)
478: *
479: *              where L(k) is the k-th column of L
480: *
481:                IF( K.LT.N ) THEN
482: *
483: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
484: *
485: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
486: *
487:                   R1 = ONE / DBLE( A( K, K ) )
488:                   CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
489:      $                       A( K+1, K+1 ), LDA )
490: *
491: *                 Store L(k) in column K
492: *
493:                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
494:                END IF
495:             ELSE
496: *
497: *              2-by-2 pivot block D(k)
498: *
499:                IF( K.LT.N-1 ) THEN
500: *
501: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
502: *
503: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
504: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
505: *
506: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
507: *                 columns of L
508: *
509:                   D = DLAPY2( DBLE( A( K+1, K ) ),
510:      $                DIMAG( A( K+1, K ) ) )
511:                   D11 = DBLE( A( K+1, K+1 ) ) / D
512:                   D22 = DBLE( A( K, K ) ) / D
513:                   TT = ONE / ( D11*D22-ONE )
514:                   D21 = A( K+1, K ) / D
515:                   D = TT / D
516: *
517:                   DO 80 J = K + 2, N
518:                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
519:                      WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
520:      $                      A( J, K ) )
521:                      DO 70 I = J, N
522:                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
523:      $                              A( I, K+1 )*DCONJG( WKP1 )
524:    70                CONTINUE
525:                      A( J, K ) = WK
526:                      A( J, K+1 ) = WKP1
527:                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
528:    80             CONTINUE
529:                END IF
530:             END IF
531:          END IF
532: *
533: *        Store details of the interchanges in IPIV
534: *
535:          IF( KSTEP.EQ.1 ) THEN
536:             IPIV( K ) = KP
537:          ELSE
538:             IPIV( K ) = -KP
539:             IPIV( K+1 ) = -KP
540:          END IF
541: *
542: *        Increase K and return to the start of the main loop
543: *
544:          K = K + KSTEP
545:          GO TO 50
546: *
547:       END IF
548: *
549:    90 CONTINUE
550:       RETURN
551: *
552: *     End of ZHETF2
553: *
554:       END
555: