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