001:       SUBROUTINE ZSYTF2( 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: *  ZSYTF2 computes the factorization of a complex symmetric 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 transpose of U, and D is symmetric and
027: *  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: *          symmetric 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 symmetric 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.209 and l.377
083: *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
084: *    by
085: *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
086: *
087: *  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
088: *         Company
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:       DOUBLE PRECISION   ZERO, ONE
128:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
129:       DOUBLE PRECISION   EIGHT, SEVTEN
130:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
131:       COMPLEX*16         CONE
132:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
133: *     ..
134: *     .. Local Scalars ..
135:       LOGICAL            UPPER
136:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
137:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, ROWMAX
138:       COMPLEX*16         D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
139: *     ..
140: *     .. External Functions ..
141:       LOGICAL            DISNAN, LSAME
142:       INTEGER            IZAMAX
143:       EXTERNAL           DISNAN, LSAME, IZAMAX
144: *     ..
145: *     .. External Subroutines ..
146:       EXTERNAL           XERBLA, ZSCAL, ZSWAP, ZSYR
147: *     ..
148: *     .. Intrinsic Functions ..
149:       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
150: *     ..
151: *     .. Statement Functions ..
152:       DOUBLE PRECISION   CABS1
153: *     ..
154: *     .. Statement Function definitions ..
155:       CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
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( 'ZSYTF2', -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 70
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 = CABS1( 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 = IZAMAX( 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. DISNAN(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:          ELSE
218:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
219: *
220: *              no interchange, use 1-by-1 pivot block
221: *
222:                KP = K
223:             ELSE
224: *
225: *              JMAX is the column-index of the largest off-diagonal
226: *              element in row IMAX, and ROWMAX is its absolute value
227: *
228:                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
229:                ROWMAX = CABS1( A( IMAX, JMAX ) )
230:                IF( IMAX.GT.1 ) THEN
231:                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
232:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
233:                END IF
234: *
235:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
236: *
237: *                 no interchange, use 1-by-1 pivot block
238: *
239:                   KP = K
240:                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
241: *
242: *                 interchange rows and columns K and IMAX, use 1-by-1
243: *                 pivot block
244: *
245:                   KP = IMAX
246:                ELSE
247: *
248: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
249: *                 pivot block
250: *
251:                   KP = IMAX
252:                   KSTEP = 2
253:                END IF
254:             END IF
255: *
256:             KK = K - KSTEP + 1
257:             IF( KP.NE.KK ) THEN
258: *
259: *              Interchange rows and columns KK and KP in the leading
260: *              submatrix A(1:k,1:k)
261: *
262:                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
263:                CALL ZSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
264:      $                     LDA )
265:                T = A( KK, KK )
266:                A( KK, KK ) = A( KP, KP )
267:                A( KP, KP ) = T
268:                IF( KSTEP.EQ.2 ) THEN
269:                   T = A( K-1, K )
270:                   A( K-1, K ) = A( KP, K )
271:                   A( KP, K ) = T
272:                END IF
273:             END IF
274: *
275: *           Update the leading submatrix
276: *
277:             IF( KSTEP.EQ.1 ) THEN
278: *
279: *              1-by-1 pivot block D(k): column k now holds
280: *
281: *              W(k) = U(k)*D(k)
282: *
283: *              where U(k) is the k-th column of U
284: *
285: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
286: *
287: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
288: *
289:                R1 = CONE / A( K, K )
290:                CALL ZSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
291: *
292: *              Store U(k) in column k
293: *
294:                CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
295:             ELSE
296: *
297: *              2-by-2 pivot block D(k): columns k and k-1 now hold
298: *
299: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
300: *
301: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
302: *              of U
303: *
304: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
305: *
306: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
307: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
308: *
309:                IF( K.GT.2 ) THEN
310: *
311:                   D12 = A( K-1, K )
312:                   D22 = A( K-1, K-1 ) / D12
313:                   D11 = A( K, K ) / D12
314:                   T = CONE / ( D11*D22-CONE )
315:                   D12 = T / D12
316: *
317:                   DO 30 J = K - 2, 1, -1
318:                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
319:                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
320:                      DO 20 I = J, 1, -1
321:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
322:      $                              A( I, K-1 )*WKM1
323:    20                CONTINUE
324:                      A( J, K ) = WK
325:                      A( J, K-1 ) = WKM1
326:    30             CONTINUE
327: *
328:                END IF
329: *
330:             END IF
331:          END IF
332: *
333: *        Store details of the interchanges in IPIV
334: *
335:          IF( KSTEP.EQ.1 ) THEN
336:             IPIV( K ) = KP
337:          ELSE
338:             IPIV( K ) = -KP
339:             IPIV( K-1 ) = -KP
340:          END IF
341: *
342: *        Decrease K and return to the start of the main loop
343: *
344:          K = K - KSTEP
345:          GO TO 10
346: *
347:       ELSE
348: *
349: *        Factorize A as L*D*L' using the lower triangle of A
350: *
351: *        K is the main loop index, increasing from 1 to N in steps of
352: *        1 or 2
353: *
354:          K = 1
355:    40    CONTINUE
356: *
357: *        If K > N, exit from loop
358: *
359:          IF( K.GT.N )
360:      $      GO TO 70
361:          KSTEP = 1
362: *
363: *        Determine rows and columns to be interchanged and whether
364: *        a 1-by-1 or 2-by-2 pivot block will be used
365: *
366:          ABSAKK = CABS1( A( K, K ) )
367: *
368: *        IMAX is the row-index of the largest off-diagonal element in
369: *        column K, and COLMAX is its absolute value
370: *
371:          IF( K.LT.N ) THEN
372:             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
373:             COLMAX = CABS1( A( IMAX, K ) )
374:          ELSE
375:             COLMAX = ZERO
376:          END IF
377: *
378:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. DISNAN(ABSAKK) ) THEN
379: *
380: *           Column K is zero or contains a NaN: set INFO and continue
381: *
382:             IF( INFO.EQ.0 )
383:      $         INFO = K
384:             KP = K
385:          ELSE
386:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
387: *
388: *              no interchange, use 1-by-1 pivot block
389: *
390:                KP = K
391:             ELSE
392: *
393: *              JMAX is the column-index of the largest off-diagonal
394: *              element in row IMAX, and ROWMAX is its absolute value
395: *
396:                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
397:                ROWMAX = CABS1( A( IMAX, JMAX ) )
398:                IF( IMAX.LT.N ) THEN
399:                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
400:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
401:                END IF
402: *
403:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
404: *
405: *                 no interchange, use 1-by-1 pivot block
406: *
407:                   KP = K
408:                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
409: *
410: *                 interchange rows and columns K and IMAX, use 1-by-1
411: *                 pivot block
412: *
413:                   KP = IMAX
414:                ELSE
415: *
416: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
417: *                 pivot block
418: *
419:                   KP = IMAX
420:                   KSTEP = 2
421:                END IF
422:             END IF
423: *
424:             KK = K + KSTEP - 1
425:             IF( KP.NE.KK ) THEN
426: *
427: *              Interchange rows and columns KK and KP in the trailing
428: *              submatrix A(k:n,k:n)
429: *
430:                IF( KP.LT.N )
431:      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
432:                CALL ZSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
433:      $                     LDA )
434:                T = A( KK, KK )
435:                A( KK, KK ) = A( KP, KP )
436:                A( KP, KP ) = T
437:                IF( KSTEP.EQ.2 ) THEN
438:                   T = A( K+1, K )
439:                   A( K+1, K ) = A( KP, K )
440:                   A( KP, K ) = T
441:                END IF
442:             END IF
443: *
444: *           Update the trailing submatrix
445: *
446:             IF( KSTEP.EQ.1 ) THEN
447: *
448: *              1-by-1 pivot block D(k): column k now holds
449: *
450: *              W(k) = L(k)*D(k)
451: *
452: *              where L(k) is the k-th column of L
453: *
454:                IF( K.LT.N ) THEN
455: *
456: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
457: *
458: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
459: *
460:                   R1 = CONE / A( K, K )
461:                   CALL ZSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
462:      $                       A( K+1, K+1 ), LDA )
463: *
464: *                 Store L(k) in column K
465: *
466:                   CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
467:                END IF
468:             ELSE
469: *
470: *              2-by-2 pivot block D(k)
471: *
472:                IF( K.LT.N-1 ) THEN
473: *
474: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
475: *
476: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
477: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
478: *
479: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
480: *                 columns of L
481: *
482:                   D21 = A( K+1, K )
483:                   D11 = A( K+1, K+1 ) / D21
484:                   D22 = A( K, K ) / D21
485:                   T = CONE / ( D11*D22-CONE )
486:                   D21 = T / D21
487: *
488:                   DO 60 J = K + 2, N
489:                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
490:                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
491:                      DO 50 I = J, N
492:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
493:      $                              A( I, K+1 )*WKP1
494:    50                CONTINUE
495:                      A( J, K ) = WK
496:                      A( J, K+1 ) = WKP1
497:    60             CONTINUE
498:                END IF
499:             END IF
500:          END IF
501: *
502: *        Store details of the interchanges in IPIV
503: *
504:          IF( KSTEP.EQ.1 ) THEN
505:             IPIV( K ) = KP
506:          ELSE
507:             IPIV( K ) = -KP
508:             IPIV( K+1 ) = -KP
509:          END IF
510: *
511: *        Increase K and return to the start of the main loop
512: *
513:          K = K + KSTEP
514:          GO TO 40
515: *
516:       END IF
517: *
518:    70 CONTINUE
519:       RETURN
520: *
521: *     End of ZSYTF2
522: *
523:       END
524: