001:       SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, 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, KB, LDA, LDW, N, NB
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       DOUBLE PRECISION   A( LDA, * ), W( LDW, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DLASYF computes a partial factorization of a real symmetric matrix A
020: *  using the Bunch-Kaufman diagonal pivoting method. The partial
021: *  factorization has the form:
022: *
023: *  A  =  ( I  U12 ) ( A11  0  ) (  I    0   )  if UPLO = 'U', or:
024: *        ( 0  U22 ) (  0   D  ) ( U12' U22' )
025: *
026: *  A  =  ( L11  0 ) (  D   0  ) ( L11' L21' )  if UPLO = 'L'
027: *        ( L21  I ) (  0  A22 ) (  0    I   )
028: *
029: *  where the order of D is at most NB. The actual order is returned in
030: *  the argument KB, and is either NB or NB-1, or N if N <= NB.
031: *
032: *  DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
033: *  (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
034: *  A22 (if UPLO = 'L').
035: *
036: *  Arguments
037: *  =========
038: *
039: *  UPLO    (input) CHARACTER*1
040: *          Specifies whether the upper or lower triangular part of the
041: *          symmetric matrix A is stored:
042: *          = 'U':  Upper triangular
043: *          = 'L':  Lower triangular
044: *
045: *  N       (input) INTEGER
046: *          The order of the matrix A.  N >= 0.
047: *
048: *  NB      (input) INTEGER
049: *          The maximum number of columns of the matrix A that should be
050: *          factored.  NB should be at least 2 to allow for 2-by-2 pivot
051: *          blocks.
052: *
053: *  KB      (output) INTEGER
054: *          The number of columns of A that were actually factored.
055: *          KB is either NB-1 or NB, or N if N <= NB.
056: *
057: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
058: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
059: *          n-by-n upper triangular part of A contains the upper
060: *          triangular part of the matrix A, and the strictly lower
061: *          triangular part of A is not referenced.  If UPLO = 'L', the
062: *          leading n-by-n lower triangular part of A contains the lower
063: *          triangular part of the matrix A, and the strictly upper
064: *          triangular part of A is not referenced.
065: *          On exit, A contains details of the partial factorization.
066: *
067: *  LDA     (input) INTEGER
068: *          The leading dimension of the array A.  LDA >= max(1,N).
069: *
070: *  IPIV    (output) INTEGER array, dimension (N)
071: *          Details of the interchanges and the block structure of D.
072: *          If UPLO = 'U', only the last KB elements of IPIV are set;
073: *          if UPLO = 'L', only the first KB elements are set.
074: *
075: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
076: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
077: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
078: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
079: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
080: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
081: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
082: *
083: *  W       (workspace) DOUBLE PRECISION array, dimension (LDW,NB)
084: *
085: *  LDW     (input) INTEGER
086: *          The leading dimension of the array W.  LDW >= max(1,N).
087: *
088: *  INFO    (output) INTEGER
089: *          = 0: successful exit
090: *          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
091: *               has been completed, but the block diagonal matrix D is
092: *               exactly singular.
093: *
094: *  =====================================================================
095: *
096: *     .. Parameters ..
097:       DOUBLE PRECISION   ZERO, ONE
098:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
099:       DOUBLE PRECISION   EIGHT, SEVTEN
100:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
101: *     ..
102: *     .. Local Scalars ..
103:       INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
104:      $                   KSTEP, KW
105:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
106:      $                   ROWMAX, T
107: *     ..
108: *     .. External Functions ..
109:       LOGICAL            LSAME
110:       INTEGER            IDAMAX
111:       EXTERNAL           LSAME, IDAMAX
112: *     ..
113: *     .. External Subroutines ..
114:       EXTERNAL           DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
115: *     ..
116: *     .. Intrinsic Functions ..
117:       INTRINSIC          ABS, MAX, MIN, SQRT
118: *     ..
119: *     .. Executable Statements ..
120: *
121:       INFO = 0
122: *
123: *     Initialize ALPHA for use in choosing pivot block size.
124: *
125:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
126: *
127:       IF( LSAME( UPLO, 'U' ) ) THEN
128: *
129: *        Factorize the trailing columns of A using the upper triangle
130: *        of A and working backwards, and compute the matrix W = U12*D
131: *        for use in updating A11
132: *
133: *        K is the main loop index, decreasing from N in steps of 1 or 2
134: *
135: *        KW is the column of W which corresponds to column K of A
136: *
137:          K = N
138:    10    CONTINUE
139:          KW = NB + K - N
140: *
141: *        Exit from loop
142: *
143:          IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
144:      $      GO TO 30
145: *
146: *        Copy column K of A to column KW of W and update it
147: *
148:          CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
149:          IF( K.LT.N )
150:      $      CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA,
151:      $                  W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
152: *
153:          KSTEP = 1
154: *
155: *        Determine rows and columns to be interchanged and whether
156: *        a 1-by-1 or 2-by-2 pivot block will be used
157: *
158:          ABSAKK = ABS( W( K, KW ) )
159: *
160: *        IMAX is the row-index of the largest off-diagonal element in
161: *        column K, and COLMAX is its absolute value
162: *
163:          IF( K.GT.1 ) THEN
164:             IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
165:             COLMAX = ABS( W( IMAX, KW ) )
166:          ELSE
167:             COLMAX = ZERO
168:          END IF
169: *
170:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
171: *
172: *           Column K is zero: set INFO and continue
173: *
174:             IF( INFO.EQ.0 )
175:      $         INFO = K
176:             KP = K
177:          ELSE
178:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
179: *
180: *              no interchange, use 1-by-1 pivot block
181: *
182:                KP = K
183:             ELSE
184: *
185: *              Copy column IMAX to column KW-1 of W and update it
186: *
187:                CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
188:                CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
189:      $                     W( IMAX+1, KW-1 ), 1 )
190:                IF( K.LT.N )
191:      $            CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ),
192:      $                        LDA, W( IMAX, KW+1 ), LDW, ONE,
193:      $                        W( 1, KW-1 ), 1 )
194: *
195: *              JMAX is the column-index of the largest off-diagonal
196: *              element in row IMAX, and ROWMAX is its absolute value
197: *
198:                JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
199:                ROWMAX = ABS( W( JMAX, KW-1 ) )
200:                IF( IMAX.GT.1 ) THEN
201:                   JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 )
202:                   ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
203:                END IF
204: *
205:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
206: *
207: *                 no interchange, use 1-by-1 pivot block
208: *
209:                   KP = K
210:                ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
211: *
212: *                 interchange rows and columns K and IMAX, use 1-by-1
213: *                 pivot block
214: *
215:                   KP = IMAX
216: *
217: *                 copy column KW-1 of W to column KW
218: *
219:                   CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
220:                ELSE
221: *
222: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
223: *                 pivot block
224: *
225:                   KP = IMAX
226:                   KSTEP = 2
227:                END IF
228:             END IF
229: *
230:             KK = K - KSTEP + 1
231:             KKW = NB + KK - N
232: *
233: *           Updated column KP is already stored in column KKW of W
234: *
235:             IF( KP.NE.KK ) THEN
236: *
237: *              Copy non-updated column KK to column KP
238: *
239:                A( KP, K ) = A( KK, K )
240:                CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
241:      $                     LDA )
242:                CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
243: *
244: *              Interchange rows KK and KP in last KK columns of A and W
245: *
246:                CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
247:                CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
248:      $                     LDW )
249:             END IF
250: *
251:             IF( KSTEP.EQ.1 ) THEN
252: *
253: *              1-by-1 pivot block D(k): column KW of W now holds
254: *
255: *              W(k) = U(k)*D(k)
256: *
257: *              where U(k) is the k-th column of U
258: *
259: *              Store U(k) in column k of A
260: *
261:                CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
262:                R1 = ONE / A( K, K )
263:                CALL DSCAL( K-1, R1, A( 1, K ), 1 )
264:             ELSE
265: *
266: *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
267: *              hold
268: *
269: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
270: *
271: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
272: *              of U
273: *
274:                IF( K.GT.2 ) THEN
275: *
276: *                 Store U(k) and U(k-1) in columns k and k-1 of A
277: *
278:                   D21 = W( K-1, KW )
279:                   D11 = W( K, KW ) / D21
280:                   D22 = W( K-1, KW-1 ) / D21
281:                   T = ONE / ( D11*D22-ONE )
282:                   D21 = T / D21
283:                   DO 20 J = 1, K - 2
284:                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
285:                      A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
286:    20             CONTINUE
287:                END IF
288: *
289: *              Copy D(k) to A
290: *
291:                A( K-1, K-1 ) = W( K-1, KW-1 )
292:                A( K-1, K ) = W( K-1, KW )
293:                A( K, K ) = W( K, KW )
294:             END IF
295:          END IF
296: *
297: *        Store details of the interchanges in IPIV
298: *
299:          IF( KSTEP.EQ.1 ) THEN
300:             IPIV( K ) = KP
301:          ELSE
302:             IPIV( K ) = -KP
303:             IPIV( K-1 ) = -KP
304:          END IF
305: *
306: *        Decrease K and return to the start of the main loop
307: *
308:          K = K - KSTEP
309:          GO TO 10
310: *
311:    30    CONTINUE
312: *
313: *        Update the upper triangle of A11 (= A(1:k,1:k)) as
314: *
315: *        A11 := A11 - U12*D*U12' = A11 - U12*W'
316: *
317: *        computing blocks of NB columns at a time
318: *
319:          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
320:             JB = MIN( NB, K-J+1 )
321: *
322: *           Update the upper triangle of the diagonal block
323: *
324:             DO 40 JJ = J, J + JB - 1
325:                CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
326:      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
327:      $                     A( J, JJ ), 1 )
328:    40       CONTINUE
329: *
330: *           Update the rectangular superdiagonal block
331: *
332:             CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE,
333:      $                  A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
334:      $                  A( 1, J ), LDA )
335:    50    CONTINUE
336: *
337: *        Put U12 in standard form by partially undoing the interchanges
338: *        in columns k+1:n
339: *
340:          J = K + 1
341:    60    CONTINUE
342:          JJ = J
343:          JP = IPIV( J )
344:          IF( JP.LT.0 ) THEN
345:             JP = -JP
346:             J = J + 1
347:          END IF
348:          J = J + 1
349:          IF( JP.NE.JJ .AND. J.LE.N )
350:      $      CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
351:          IF( J.LE.N )
352:      $      GO TO 60
353: *
354: *        Set KB to the number of columns factorized
355: *
356:          KB = N - K
357: *
358:       ELSE
359: *
360: *        Factorize the leading columns of A using the lower triangle
361: *        of A and working forwards, and compute the matrix W = L21*D
362: *        for use in updating A22
363: *
364: *        K is the main loop index, increasing from 1 in steps of 1 or 2
365: *
366:          K = 1
367:    70    CONTINUE
368: *
369: *        Exit from loop
370: *
371:          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
372:      $      GO TO 90
373: *
374: *        Copy column K of A to column K of W and update it
375: *
376:          CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
377:          CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA,
378:      $               W( K, 1 ), LDW, ONE, W( K, K ), 1 )
379: *
380:          KSTEP = 1
381: *
382: *        Determine rows and columns to be interchanged and whether
383: *        a 1-by-1 or 2-by-2 pivot block will be used
384: *
385:          ABSAKK = ABS( W( K, K ) )
386: *
387: *        IMAX is the row-index of the largest off-diagonal element in
388: *        column K, and COLMAX is its absolute value
389: *
390:          IF( K.LT.N ) THEN
391:             IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
392:             COLMAX = ABS( W( IMAX, K ) )
393:          ELSE
394:             COLMAX = ZERO
395:          END IF
396: *
397:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
398: *
399: *           Column K is zero: set INFO and continue
400: *
401:             IF( INFO.EQ.0 )
402:      $         INFO = K
403:             KP = K
404:          ELSE
405:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
406: *
407: *              no interchange, use 1-by-1 pivot block
408: *
409:                KP = K
410:             ELSE
411: *
412: *              Copy column IMAX to column K+1 of W and update it
413: *
414:                CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
415:                CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
416:      $                     1 )
417:                CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ),
418:      $                     LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
419: *
420: *              JMAX is the column-index of the largest off-diagonal
421: *              element in row IMAX, and ROWMAX is its absolute value
422: *
423:                JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 )
424:                ROWMAX = ABS( W( JMAX, K+1 ) )
425:                IF( IMAX.LT.N ) THEN
426:                   JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
427:                   ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
428:                END IF
429: *
430:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
431: *
432: *                 no interchange, use 1-by-1 pivot block
433: *
434:                   KP = K
435:                ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
436: *
437: *                 interchange rows and columns K and IMAX, use 1-by-1
438: *                 pivot block
439: *
440:                   KP = IMAX
441: *
442: *                 copy column K+1 of W to column K
443: *
444:                   CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
445:                ELSE
446: *
447: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
448: *                 pivot block
449: *
450:                   KP = IMAX
451:                   KSTEP = 2
452:                END IF
453:             END IF
454: *
455:             KK = K + KSTEP - 1
456: *
457: *           Updated column KP is already stored in column KK of W
458: *
459:             IF( KP.NE.KK ) THEN
460: *
461: *              Copy non-updated column KK to column KP
462: *
463:                A( KP, K ) = A( KK, K )
464:                CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
465:                CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
466: *
467: *              Interchange rows KK and KP in first KK columns of A and W
468: *
469:                CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
470:                CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
471:             END IF
472: *
473:             IF( KSTEP.EQ.1 ) THEN
474: *
475: *              1-by-1 pivot block D(k): column k of W now holds
476: *
477: *              W(k) = L(k)*D(k)
478: *
479: *              where L(k) is the k-th column of L
480: *
481: *              Store L(k) in column k of A
482: *
483:                CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
484:                IF( K.LT.N ) THEN
485:                   R1 = ONE / A( K, K )
486:                   CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
487:                END IF
488:             ELSE
489: *
490: *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
491: *
492: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
493: *
494: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
495: *              of L
496: *
497:                IF( K.LT.N-1 ) THEN
498: *
499: *                 Store L(k) and L(k+1) in columns k and k+1 of A
500: *
501:                   D21 = W( K+1, K )
502:                   D11 = W( K+1, K+1 ) / D21
503:                   D22 = W( K, K ) / D21
504:                   T = ONE / ( D11*D22-ONE )
505:                   D21 = T / D21
506:                   DO 80 J = K + 2, N
507:                      A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
508:                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
509:    80             CONTINUE
510:                END IF
511: *
512: *              Copy D(k) to A
513: *
514:                A( K, K ) = W( K, K )
515:                A( K+1, K ) = W( K+1, K )
516:                A( K+1, K+1 ) = W( K+1, K+1 )
517:             END IF
518:          END IF
519: *
520: *        Store details of the interchanges in IPIV
521: *
522:          IF( KSTEP.EQ.1 ) THEN
523:             IPIV( K ) = KP
524:          ELSE
525:             IPIV( K ) = -KP
526:             IPIV( K+1 ) = -KP
527:          END IF
528: *
529: *        Increase K and return to the start of the main loop
530: *
531:          K = K + KSTEP
532:          GO TO 70
533: *
534:    90    CONTINUE
535: *
536: *        Update the lower triangle of A22 (= A(k:n,k:n)) as
537: *
538: *        A22 := A22 - L21*D*L21' = A22 - L21*W'
539: *
540: *        computing blocks of NB columns at a time
541: *
542:          DO 110 J = K, N, NB
543:             JB = MIN( NB, N-J+1 )
544: *
545: *           Update the lower triangle of the diagonal block
546: *
547:             DO 100 JJ = J, J + JB - 1
548:                CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
549:      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
550:      $                     A( JJ, JJ ), 1 )
551:   100       CONTINUE
552: *
553: *           Update the rectangular subdiagonal block
554: *
555:             IF( J+JB.LE.N )
556:      $         CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
557:      $                     K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
558:      $                     ONE, A( J+JB, J ), LDA )
559:   110    CONTINUE
560: *
561: *        Put L21 in standard form by partially undoing the interchanges
562: *        in columns 1:k-1
563: *
564:          J = K - 1
565:   120    CONTINUE
566:          JJ = J
567:          JP = IPIV( J )
568:          IF( JP.LT.0 ) THEN
569:             JP = -JP
570:             J = J - 1
571:          END IF
572:          J = J - 1
573:          IF( JP.NE.JJ .AND. J.GE.1 )
574:      $      CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
575:          IF( J.GE.1 )
576:      $      GO TO 120
577: *
578: *        Set KB to the number of columns factorized
579: *
580:          KB = K - 1
581: *
582:       END IF
583:       RETURN
584: *
585: *     End of DLASYF
586: *
587:       END
588: