001:       SUBROUTINE DSYTF2( 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:       DOUBLE PRECISION   A( LDA, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DSYTF2 computes the factorization of a real symmetric matrix A using
021: *  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) DOUBLE PRECISION 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.204 and l.372
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: *  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
091: *         Company
092: *
093: *  If UPLO = 'U', then A = U*D*U', where
094: *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
095: *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
096: *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
097: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
098: *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
099: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
100: *
101: *             (   I    v    0   )   k-s
102: *     U(k) =  (   0    I    0   )   s
103: *             (   0    0    I   )   n-k
104: *                k-s   s   n-k
105: *
106: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
107: *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
108: *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
109: *
110: *  If UPLO = 'L', then A = L*D*L', where
111: *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
112: *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
113: *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
114: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
115: *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
116: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
117: *
118: *             (   I    0     0   )  k-1
119: *     L(k) =  (   0    I     0   )  s
120: *             (   0    v     I   )  n-k-s+1
121: *                k-1   s  n-k-s+1
122: *
123: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
124: *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
125: *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
126: *
127: *  =====================================================================
128: *
129: *     .. Parameters ..
130:       DOUBLE PRECISION   ZERO, ONE
131:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
132:       DOUBLE PRECISION   EIGHT, SEVTEN
133:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
134: *     ..
135: *     .. Local Scalars ..
136:       LOGICAL            UPPER
137:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
138:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
139:      $                   ROWMAX, T, WK, WKM1, WKP1
140: *     ..
141: *     .. External Functions ..
142:       LOGICAL            LSAME, DISNAN
143:       INTEGER            IDAMAX
144:       EXTERNAL           LSAME, IDAMAX, DISNAN
145: *     ..
146: *     .. External Subroutines ..
147:       EXTERNAL           DSCAL, DSWAP, DSYR, XERBLA
148: *     ..
149: *     .. Intrinsic Functions ..
150:       INTRINSIC          ABS, MAX, SQRT
151: *     ..
152: *     .. Executable Statements ..
153: *
154: *     Test the input parameters.
155: *
156:       INFO = 0
157:       UPPER = LSAME( UPLO, 'U' )
158:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
159:          INFO = -1
160:       ELSE IF( N.LT.0 ) THEN
161:          INFO = -2
162:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
163:          INFO = -4
164:       END IF
165:       IF( INFO.NE.0 ) THEN
166:          CALL XERBLA( 'DSYTF2', -INFO )
167:          RETURN
168:       END IF
169: *
170: *     Initialize ALPHA for use in choosing pivot block size.
171: *
172:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
173: *
174:       IF( UPPER ) THEN
175: *
176: *        Factorize A as U*D*U' using the upper triangle of A
177: *
178: *        K is the main loop index, decreasing from N to 1 in steps of
179: *        1 or 2
180: *
181:          K = N
182:    10    CONTINUE
183: *
184: *        If K < 1, exit from loop
185: *
186:          IF( K.LT.1 )
187:      $      GO TO 70
188:          KSTEP = 1
189: *
190: *        Determine rows and columns to be interchanged and whether
191: *        a 1-by-1 or 2-by-2 pivot block will be used
192: *
193:          ABSAKK = ABS( A( K, K ) )
194: *
195: *        IMAX is the row-index of the largest off-diagonal element in
196: *        column K, and COLMAX is its absolute value
197: *
198:          IF( K.GT.1 ) THEN
199:             IMAX = IDAMAX( K-1, A( 1, K ), 1 )
200:             COLMAX = ABS( A( IMAX, K ) )
201:          ELSE
202:             COLMAX = ZERO
203:          END IF
204: *
205:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
206: *
207: *           Column K is zero or contains a NaN: set INFO and continue
208: *
209:             IF( INFO.EQ.0 )
210:      $         INFO = K
211:             KP = K
212:          ELSE
213:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
214: *
215: *              no interchange, use 1-by-1 pivot block
216: *
217:                KP = K
218:             ELSE
219: *
220: *              JMAX is the column-index of the largest off-diagonal
221: *              element in row IMAX, and ROWMAX is its absolute value
222: *
223:                JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
224:                ROWMAX = ABS( A( IMAX, JMAX ) )
225:                IF( IMAX.GT.1 ) THEN
226:                   JMAX = IDAMAX( IMAX-1, A( 1, IMAX ), 1 )
227:                   ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
228:                END IF
229: *
230:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
231: *
232: *                 no interchange, use 1-by-1 pivot block
233: *
234:                   KP = K
235:                ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
236: *
237: *                 interchange rows and columns K and IMAX, use 1-by-1
238: *                 pivot block
239: *
240:                   KP = IMAX
241:                ELSE
242: *
243: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
244: *                 pivot block
245: *
246:                   KP = IMAX
247:                   KSTEP = 2
248:                END IF
249:             END IF
250: *
251:             KK = K - KSTEP + 1
252:             IF( KP.NE.KK ) THEN
253: *
254: *              Interchange rows and columns KK and KP in the leading
255: *              submatrix A(1:k,1:k)
256: *
257:                CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
258:                CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
259:      $                     LDA )
260:                T = A( KK, KK )
261:                A( KK, KK ) = A( KP, KP )
262:                A( KP, KP ) = T
263:                IF( KSTEP.EQ.2 ) THEN
264:                   T = A( K-1, K )
265:                   A( K-1, K ) = A( KP, K )
266:                   A( KP, K ) = T
267:                END IF
268:             END IF
269: *
270: *           Update the leading submatrix
271: *
272:             IF( KSTEP.EQ.1 ) THEN
273: *
274: *              1-by-1 pivot block D(k): column k now holds
275: *
276: *              W(k) = U(k)*D(k)
277: *
278: *              where U(k) is the k-th column of U
279: *
280: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
281: *
282: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
283: *
284:                R1 = ONE / A( K, K )
285:                CALL DSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
286: *
287: *              Store U(k) in column k
288: *
289:                CALL DSCAL( K-1, R1, A( 1, K ), 1 )
290:             ELSE
291: *
292: *              2-by-2 pivot block D(k): columns k and k-1 now hold
293: *
294: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
295: *
296: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
297: *              of U
298: *
299: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
300: *
301: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
302: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
303: *
304:                IF( K.GT.2 ) THEN
305: *
306:                   D12 = A( K-1, K )
307:                   D22 = A( K-1, K-1 ) / D12
308:                   D11 = A( K, K ) / D12
309:                   T = ONE / ( D11*D22-ONE )
310:                   D12 = T / D12
311: *
312:                   DO 30 J = K - 2, 1, -1
313:                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
314:                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
315:                      DO 20 I = J, 1, -1
316:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
317:      $                              A( I, K-1 )*WKM1
318:    20                CONTINUE
319:                      A( J, K ) = WK
320:                      A( J, K-1 ) = WKM1
321:    30             CONTINUE
322: *
323:                END IF
324: *
325:             END IF
326:          END IF
327: *
328: *        Store details of the interchanges in IPIV
329: *
330:          IF( KSTEP.EQ.1 ) THEN
331:             IPIV( K ) = KP
332:          ELSE
333:             IPIV( K ) = -KP
334:             IPIV( K-1 ) = -KP
335:          END IF
336: *
337: *        Decrease K and return to the start of the main loop
338: *
339:          K = K - KSTEP
340:          GO TO 10
341: *
342:       ELSE
343: *
344: *        Factorize A as L*D*L' using the lower triangle of A
345: *
346: *        K is the main loop index, increasing from 1 to N in steps of
347: *        1 or 2
348: *
349:          K = 1
350:    40    CONTINUE
351: *
352: *        If K > N, exit from loop
353: *
354:          IF( K.GT.N )
355:      $      GO TO 70
356:          KSTEP = 1
357: *
358: *        Determine rows and columns to be interchanged and whether
359: *        a 1-by-1 or 2-by-2 pivot block will be used
360: *
361:          ABSAKK = ABS( A( K, K ) )
362: *
363: *        IMAX is the row-index of the largest off-diagonal element in
364: *        column K, and COLMAX is its absolute value
365: *
366:          IF( K.LT.N ) THEN
367:             IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 )
368:             COLMAX = ABS( A( IMAX, K ) )
369:          ELSE
370:             COLMAX = ZERO
371:          END IF
372: *
373:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
374: *
375: *           Column K is zero or contains a NaN: set INFO and continue
376: *
377:             IF( INFO.EQ.0 )
378:      $         INFO = K
379:             KP = K
380:          ELSE
381:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
382: *
383: *              no interchange, use 1-by-1 pivot block
384: *
385:                KP = K
386:             ELSE
387: *
388: *              JMAX is the column-index of the largest off-diagonal
389: *              element in row IMAX, and ROWMAX is its absolute value
390: *
391:                JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA )
392:                ROWMAX = ABS( A( IMAX, JMAX ) )
393:                IF( IMAX.LT.N ) THEN
394:                   JMAX = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
395:                   ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
396:                END IF
397: *
398:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
399: *
400: *                 no interchange, use 1-by-1 pivot block
401: *
402:                   KP = K
403:                ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
404: *
405: *                 interchange rows and columns K and IMAX, use 1-by-1
406: *                 pivot block
407: *
408:                   KP = IMAX
409:                ELSE
410: *
411: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
412: *                 pivot block
413: *
414:                   KP = IMAX
415:                   KSTEP = 2
416:                END IF
417:             END IF
418: *
419:             KK = K + KSTEP - 1
420:             IF( KP.NE.KK ) THEN
421: *
422: *              Interchange rows and columns KK and KP in the trailing
423: *              submatrix A(k:n,k:n)
424: *
425:                IF( KP.LT.N )
426:      $            CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
427:                CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
428:      $                     LDA )
429:                T = A( KK, KK )
430:                A( KK, KK ) = A( KP, KP )
431:                A( KP, KP ) = T
432:                IF( KSTEP.EQ.2 ) THEN
433:                   T = A( K+1, K )
434:                   A( K+1, K ) = A( KP, K )
435:                   A( KP, K ) = T
436:                END IF
437:             END IF
438: *
439: *           Update the trailing submatrix
440: *
441:             IF( KSTEP.EQ.1 ) THEN
442: *
443: *              1-by-1 pivot block D(k): column k now holds
444: *
445: *              W(k) = L(k)*D(k)
446: *
447: *              where L(k) is the k-th column of L
448: *
449:                IF( K.LT.N ) THEN
450: *
451: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
452: *
453: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
454: *
455:                   D11 = ONE / A( K, K )
456:                   CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
457:      $                       A( K+1, K+1 ), LDA )
458: *
459: *                 Store L(k) in column K
460: *
461:                   CALL DSCAL( N-K, D11, A( K+1, K ), 1 )
462:                END IF
463:             ELSE
464: *
465: *              2-by-2 pivot block D(k)
466: *
467:                IF( K.LT.N-1 ) THEN
468: *
469: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
470: *
471: *                 A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))'
472: *
473: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
474: *                 columns of L
475: *
476:                   D21 = A( K+1, K )
477:                   D11 = A( K+1, K+1 ) / D21
478:                   D22 = A( K, K ) / D21
479:                   T = ONE / ( D11*D22-ONE )
480:                   D21 = T / D21
481: *
482:                   DO 60 J = K + 2, N
483: *
484:                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
485:                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
486: *
487:                      DO 50 I = J, N
488:                         A( I, J ) = A( I, J ) - A( I, K )*WK -
489:      $                              A( I, K+1 )*WKP1
490:    50                CONTINUE
491: *
492:                      A( J, K ) = WK
493:                      A( J, K+1 ) = WKP1
494: *
495:    60             CONTINUE
496:                END IF
497:             END IF
498:          END IF
499: *
500: *        Store details of the interchanges in IPIV
501: *
502:          IF( KSTEP.EQ.1 ) THEN
503:             IPIV( K ) = KP
504:          ELSE
505:             IPIV( K ) = -KP
506:             IPIV( K+1 ) = -KP
507:          END IF
508: *
509: *        Increase K and return to the start of the main loop
510: *
511:          K = K + KSTEP
512:          GO TO 40
513: *
514:       END IF
515: *
516:    70 CONTINUE
517: *
518:       RETURN
519: *
520: *     End of DSYTF2
521: *
522:       END
523: