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