001:       SUBROUTINE CHETF2( 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            A( LDA, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  CHETF2 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 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.392
083: *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
084: *    by
085: *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(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:       REAL               ZERO, ONE
129:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
130:       REAL               EIGHT, SEVTEN
131:       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
132: *     ..
133: *     .. Local Scalars ..
134:       LOGICAL            UPPER
135:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
136:       REAL               ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
137:      $                   TT
138:       COMPLEX            D12, D21, T, WK, WKM1, WKP1, ZDUM
139: *     ..
140: *     .. External Functions ..
141:       LOGICAL            LSAME, SISNAN
142:       INTEGER            ICAMAX
143:       REAL               SLAPY2
144:       EXTERNAL           LSAME, ICAMAX, SLAPY2, SISNAN
145: *     ..
146: *     .. External Subroutines ..
147:       EXTERNAL           CHER, CSSCAL, CSWAP, XERBLA
148: *     ..
149: *     .. Intrinsic Functions ..
150:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
151: *     ..
152: *     .. Statement Functions ..
153:       REAL               CABS1
154: *     ..
155: *     .. Statement Function definitions ..
156:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( 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( 'CHETF2', -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( REAL( 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 = ICAMAX( 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. SISNAN(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 ) = REAL( 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 + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
231:                ROWMAX = CABS1( A( IMAX, JMAX ) )
232:                IF( IMAX.GT.1 ) THEN
233:                   JMAX = ICAMAX( 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( REAL( 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 CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
266:                DO 20 J = KP + 1, KK - 1
267:                   T = CONJG( A( J, KK ) )
268:                   A( J, KK ) = CONJG( A( KP, J ) )
269:                   A( KP, J ) = T
270:    20          CONTINUE
271:                A( KP, KK ) = CONJG( A( KP, KK ) )
272:                R1 = REAL( A( KK, KK ) )
273:                A( KK, KK ) = REAL( A( KP, KP ) )
274:                A( KP, KP ) = R1
275:                IF( KSTEP.EQ.2 ) THEN
276:                   A( K, K ) = REAL( 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 ) = REAL( A( K, K ) )
283:                IF( KSTEP.EQ.2 )
284:      $            A( K-1, K-1 ) = REAL( 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 / REAL( A( K, K ) )
302:                CALL CHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
303: *
304: *              Store U(k) in column k
305: *
306:                CALL CSSCAL( 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 = SLAPY2( REAL( A( K-1, K ) ),
324:      $                AIMAG( A( K-1, K ) ) )
325:                   D22 = REAL( A( K-1, K-1 ) ) / D
326:                   D11 = REAL( 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 )-CONJG( D12 )*A( J, K ) )
333:                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
334:                      DO 30 I = J, 1, -1
335:                         A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
336:      $                              A( I, K-1 )*CONJG( WKM1 )
337:    30                CONTINUE
338:                      A( J, K ) = WK
339:                      A( J, K-1 ) = WKM1
340:                      A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
341:    40             CONTINUE
342: *
343:                END IF
344: *
345:             END IF
346:          END IF
347: *
348: *        Store details of the interchanges in IPIV
349: *
350:          IF( KSTEP.EQ.1 ) THEN
351:             IPIV( K ) = KP
352:          ELSE
353:             IPIV( K ) = -KP
354:             IPIV( K-1 ) = -KP
355:          END IF
356: *
357: *        Decrease K and return to the start of the main loop
358: *
359:          K = K - KSTEP
360:          GO TO 10
361: *
362:       ELSE
363: *
364: *        Factorize A as L*D*L' using the lower triangle of A
365: *
366: *        K is the main loop index, increasing from 1 to N in steps of
367: *        1 or 2
368: *
369:          K = 1
370:    50    CONTINUE
371: *
372: *        If K > N, exit from loop
373: *
374:          IF( K.GT.N )
375:      $      GO TO 90
376:          KSTEP = 1
377: *
378: *        Determine rows and columns to be interchanged and whether
379: *        a 1-by-1 or 2-by-2 pivot block will be used
380: *
381:          ABSAKK = ABS( REAL( A( K, K ) ) )
382: *
383: *        IMAX is the row-index of the largest off-diagonal element in
384: *        column K, and COLMAX is its absolute value
385: *
386:          IF( K.LT.N ) THEN
387:             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
388:             COLMAX = CABS1( A( IMAX, K ) )
389:          ELSE
390:             COLMAX = ZERO
391:          END IF
392: *
393:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
394: *
395: *           Column K is zero or contains a NaN: set INFO and continue
396: *
397:             IF( INFO.EQ.0 )
398:      $         INFO = K
399:             KP = K
400:             A( K, K ) = REAL( A( K, K ) )
401:          ELSE
402:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
403: *
404: *              no interchange, use 1-by-1 pivot block
405: *
406:                KP = K
407:             ELSE
408: *
409: *              JMAX is the column-index of the largest off-diagonal
410: *              element in row IMAX, and ROWMAX is its absolute value
411: *
412:                JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
413:                ROWMAX = CABS1( A( IMAX, JMAX ) )
414:                IF( IMAX.LT.N ) THEN
415:                   JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
416:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
417:                END IF
418: *
419:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
420: *
421: *                 no interchange, use 1-by-1 pivot block
422: *
423:                   KP = K
424:                ELSE IF( ABS( REAL( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
425:      $                   THEN
426: *
427: *                 interchange rows and columns K and IMAX, use 1-by-1
428: *                 pivot block
429: *
430:                   KP = IMAX
431:                ELSE
432: *
433: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
434: *                 pivot block
435: *
436:                   KP = IMAX
437:                   KSTEP = 2
438:                END IF
439:             END IF
440: *
441:             KK = K + KSTEP - 1
442:             IF( KP.NE.KK ) THEN
443: *
444: *              Interchange rows and columns KK and KP in the trailing
445: *              submatrix A(k:n,k:n)
446: *
447:                IF( KP.LT.N )
448:      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
449:                DO 60 J = KK + 1, KP - 1
450:                   T = CONJG( A( J, KK ) )
451:                   A( J, KK ) = CONJG( A( KP, J ) )
452:                   A( KP, J ) = T
453:    60          CONTINUE
454:                A( KP, KK ) = CONJG( A( KP, KK ) )
455:                R1 = REAL( A( KK, KK ) )
456:                A( KK, KK ) = REAL( A( KP, KP ) )
457:                A( KP, KP ) = R1
458:                IF( KSTEP.EQ.2 ) THEN
459:                   A( K, K ) = REAL( A( K, K ) )
460:                   T = A( K+1, K )
461:                   A( K+1, K ) = A( KP, K )
462:                   A( KP, K ) = T
463:                END IF
464:             ELSE
465:                A( K, K ) = REAL( A( K, K ) )
466:                IF( KSTEP.EQ.2 )
467:      $            A( K+1, K+1 ) = REAL( A( K+1, K+1 ) )
468:             END IF
469: *
470: *           Update the trailing submatrix
471: *
472:             IF( KSTEP.EQ.1 ) THEN
473: *
474: *              1-by-1 pivot block D(k): column k now holds
475: *
476: *              W(k) = L(k)*D(k)
477: *
478: *              where L(k) is the k-th column of L
479: *
480:                IF( K.LT.N ) THEN
481: *
482: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
483: *
484: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
485: *
486:                   R1 = ONE / REAL( A( K, K ) )
487:                   CALL CHER( UPLO, N-K, -R1, A( K+1, K ), 1,
488:      $                       A( K+1, K+1 ), LDA )
489: *
490: *                 Store L(k) in column K
491: *
492:                   CALL CSSCAL( N-K, R1, A( K+1, K ), 1 )
493:                END IF
494:             ELSE
495: *
496: *              2-by-2 pivot block D(k)
497: *
498:                IF( K.LT.N-1 ) THEN
499: *
500: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
501: *
502: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
503: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
504: *
505: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
506: *                 columns of L
507: *
508:                   D = SLAPY2( REAL( A( K+1, K ) ),
509:      $                        AIMAG( A( K+1, K ) ) )
510:                   D11 = REAL( A( K+1, K+1 ) ) / D
511:                   D22 = REAL( A( K, K ) ) / D
512:                   TT = ONE / ( D11*D22-ONE )
513:                   D21 = A( K+1, K ) / D
514:                   D =  TT / D
515: *
516:                   DO 80 J = K + 2, N
517:                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
518:                      WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) )
519:                      DO 70 I = J, N
520:                         A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
521:      $                              A( I, K+1 )*CONJG( WKP1 )
522:    70                CONTINUE
523:                      A( J, K ) = WK
524:                      A( J, K+1 ) = WKP1
525:                      A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
526:    80             CONTINUE
527:                END IF
528:             END IF
529:          END IF
530: *
531: *        Store details of the interchanges in IPIV
532: *
533:          IF( KSTEP.EQ.1 ) THEN
534:             IPIV( K ) = KP
535:          ELSE
536:             IPIV( K ) = -KP
537:             IPIV( K+1 ) = -KP
538:          END IF
539: *
540: *        Increase K and return to the start of the main loop
541: *
542:          K = K + KSTEP
543:          GO TO 50
544: *
545:       END IF
546: *
547:    90 CONTINUE
548:       RETURN
549: *
550: *     End of CHETF2
551: *
552:       END
553: