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