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