001:       SUBROUTINE CLAHEF( 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: *  CLAHEF computes a partial factorization of a complex Hermitian
021: *  matrix A using the Bunch-Kaufman diagonal pivoting method. The
022: *  partial 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 conjugate transpose of U.
033: *
034: *  CLAHEF is an auxiliary routine called by CHETRF. 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: *          Hermitian 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 Hermitian 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:       COMPLEX            CONE
102:       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
103:       REAL               EIGHT, SEVTEN
104:       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.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, R1, ROWMAX, T
110:       COMPLEX            D11, D21, D22, Z
111: *     ..
112: *     .. External Functions ..
113:       LOGICAL            LSAME
114:       INTEGER            ICAMAX
115:       EXTERNAL           LSAME, ICAMAX
116: *     ..
117: *     .. External Subroutines ..
118:       EXTERNAL           CCOPY, CGEMM, CGEMV, CLACGV, CSSCAL, CSWAP
119: *     ..
120: *     .. Intrinsic Functions ..
121:       INTRINSIC          ABS, AIMAG, CONJG, 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 (note that conjg(W) is actually stored)
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-1, A( 1, K ), 1, W( 1, KW ), 1 )
159:          W( K, KW ) = REAL( A( K, K ) )
160:          IF( K.LT.N ) THEN
161:             CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
162:      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
163:             W( K, KW ) = REAL( W( K, KW ) )
164:          END IF
165: *
166:          KSTEP = 1
167: *
168: *        Determine rows and columns to be interchanged and whether
169: *        a 1-by-1 or 2-by-2 pivot block will be used
170: *
171:          ABSAKK = ABS( REAL( W( K, KW ) ) )
172: *
173: *        IMAX is the row-index of the largest off-diagonal element in
174: *        column K, and COLMAX is its absolute value
175: *
176:          IF( K.GT.1 ) THEN
177:             IMAX = ICAMAX( K-1, W( 1, KW ), 1 )
178:             COLMAX = CABS1( W( IMAX, KW ) )
179:          ELSE
180:             COLMAX = ZERO
181:          END IF
182: *
183:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
184: *
185: *           Column K is zero: set INFO and continue
186: *
187:             IF( INFO.EQ.0 )
188:      $         INFO = K
189:             KP = K
190:             A( K, K ) = REAL( A( K, K ) )
191:          ELSE
192:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
193: *
194: *              no interchange, use 1-by-1 pivot block
195: *
196:                KP = K
197:             ELSE
198: *
199: *              Copy column IMAX to column KW-1 of W and update it
200: *
201:                CALL CCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
202:                W( IMAX, KW-1 ) = REAL( A( IMAX, IMAX ) )
203:                CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
204:      $                     W( IMAX+1, KW-1 ), 1 )
205:                CALL CLACGV( K-IMAX, W( IMAX+1, KW-1 ), 1 )
206:                IF( K.LT.N ) THEN
207:                   CALL CGEMV( 'No transpose', K, N-K, -CONE,
208:      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
209:      $                        CONE, W( 1, KW-1 ), 1 )
210:                   W( IMAX, KW-1 ) = REAL( W( IMAX, KW-1 ) )
211:                END IF
212: *
213: *              JMAX is the column-index of the largest off-diagonal
214: *              element in row IMAX, and ROWMAX is its absolute value
215: *
216:                JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
217:                ROWMAX = CABS1( W( JMAX, KW-1 ) )
218:                IF( IMAX.GT.1 ) THEN
219:                   JMAX = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 )
220:                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
221:                END IF
222: *
223:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
224: *
225: *                 no interchange, use 1-by-1 pivot block
226: *
227:                   KP = K
228:                ELSE IF( ABS( REAL( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )
229:      $                   THEN
230: *
231: *                 interchange rows and columns K and IMAX, use 1-by-1
232: *                 pivot block
233: *
234:                   KP = IMAX
235: *
236: *                 copy column KW-1 of W to column KW
237: *
238:                   CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
239:                ELSE
240: *
241: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
242: *                 pivot block
243: *
244:                   KP = IMAX
245:                   KSTEP = 2
246:                END IF
247:             END IF
248: *
249:             KK = K - KSTEP + 1
250:             KKW = NB + KK - N
251: *
252: *           Updated column KP is already stored in column KKW of W
253: *
254:             IF( KP.NE.KK ) THEN
255: *
256: *              Copy non-updated column KK to column KP
257: *
258:                A( KP, KP ) = REAL( A( KK, KK ) )
259:                CALL CCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
260:      $                     LDA )
261:                CALL CLACGV( KK-1-KP, A( KP, KP+1 ), LDA )
262:                CALL CCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
263: *
264: *              Interchange rows KK and KP in last KK columns of A and W
265: *
266:                IF( KK.LT.N )
267:      $            CALL CSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ),
268:      $                        LDA )
269:                CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
270:      $                     LDW )
271:             END IF
272: *
273:             IF( KSTEP.EQ.1 ) THEN
274: *
275: *              1-by-1 pivot block D(k): column KW of W now holds
276: *
277: *              W(k) = U(k)*D(k)
278: *
279: *              where U(k) is the k-th column of U
280: *
281: *              Store U(k) in column k of A
282: *
283:                CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
284:                R1 = ONE / REAL( A( K, K ) )
285:                CALL CSSCAL( K-1, R1, A( 1, K ), 1 )
286: *
287: *              Conjugate W(k)
288: *
289:                CALL CLACGV( K-1, W( 1, KW ), 1 )
290:             ELSE
291: *
292: *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
293: *              hold
294: *
295: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
296: *
297: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
298: *              of U
299: *
300:                IF( K.GT.2 ) THEN
301: *
302: *                 Store U(k) and U(k-1) in columns k and k-1 of A
303: *
304:                   D21 = W( K-1, KW )
305:                   D11 = W( K, KW ) / CONJG( D21 )
306:                   D22 = W( K-1, KW-1 ) / D21
307:                   T = ONE / ( REAL( D11*D22 )-ONE )
308:                   D21 = T / D21
309:                   DO 20 J = 1, K - 2
310:                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
311:                      A( J, K ) = CONJG( D21 )*
312:      $                           ( D22*W( J, KW )-W( J, KW-1 ) )
313:    20             CONTINUE
314:                END IF
315: *
316: *              Copy D(k) to A
317: *
318:                A( K-1, K-1 ) = W( K-1, KW-1 )
319:                A( K-1, K ) = W( K-1, KW )
320:                A( K, K ) = W( K, KW )
321: *
322: *              Conjugate W(k) and W(k-1)
323: *
324:                CALL CLACGV( K-1, W( 1, KW ), 1 )
325:                CALL CLACGV( K-2, W( 1, KW-1 ), 1 )
326:             END IF
327:          END IF
328: *
329: *        Store details of the interchanges in IPIV
330: *
331:          IF( KSTEP.EQ.1 ) THEN
332:             IPIV( K ) = KP
333:          ELSE
334:             IPIV( K ) = -KP
335:             IPIV( K-1 ) = -KP
336:          END IF
337: *
338: *        Decrease K and return to the start of the main loop
339: *
340:          K = K - KSTEP
341:          GO TO 10
342: *
343:    30    CONTINUE
344: *
345: *        Update the upper triangle of A11 (= A(1:k,1:k)) as
346: *
347: *        A11 := A11 - U12*D*U12' = A11 - U12*W'
348: *
349: *        computing blocks of NB columns at a time (note that conjg(W) is
350: *        actually stored)
351: *
352:          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
353:             JB = MIN( NB, K-J+1 )
354: *
355: *           Update the upper triangle of the diagonal block
356: *
357:             DO 40 JJ = J, J + JB - 1
358:                A( JJ, JJ ) = REAL( A( JJ, JJ ) )
359:                CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
360:      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
361:      $                     A( J, JJ ), 1 )
362:                A( JJ, JJ ) = REAL( A( JJ, JJ ) )
363:    40       CONTINUE
364: *
365: *           Update the rectangular superdiagonal block
366: *
367:             CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
368:      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
369:      $                  CONE, A( 1, J ), LDA )
370:    50    CONTINUE
371: *
372: *        Put U12 in standard form by partially undoing the interchanges
373: *        in columns k+1:n
374: *
375:          J = K + 1
376:    60    CONTINUE
377:          JJ = J
378:          JP = IPIV( J )
379:          IF( JP.LT.0 ) THEN
380:             JP = -JP
381:             J = J + 1
382:          END IF
383:          J = J + 1
384:          IF( JP.NE.JJ .AND. J.LE.N )
385:      $      CALL CSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
386:          IF( J.LE.N )
387:      $      GO TO 60
388: *
389: *        Set KB to the number of columns factorized
390: *
391:          KB = N - K
392: *
393:       ELSE
394: *
395: *        Factorize the leading columns of A using the lower triangle
396: *        of A and working forwards, and compute the matrix W = L21*D
397: *        for use in updating A22 (note that conjg(W) is actually stored)
398: *
399: *        K is the main loop index, increasing from 1 in steps of 1 or 2
400: *
401:          K = 1
402:    70    CONTINUE
403: *
404: *        Exit from loop
405: *
406:          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
407:      $      GO TO 90
408: *
409: *        Copy column K of A to column K of W and update it
410: *
411:          W( K, K ) = REAL( A( K, K ) )
412:          IF( K.LT.N )
413:      $      CALL CCOPY( N-K, A( K+1, K ), 1, W( K+1, K ), 1 )
414:          CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
415:      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
416:          W( K, K ) = REAL( W( K, K ) )
417: *
418:          KSTEP = 1
419: *
420: *        Determine rows and columns to be interchanged and whether
421: *        a 1-by-1 or 2-by-2 pivot block will be used
422: *
423:          ABSAKK = ABS( REAL( W( K, K ) ) )
424: *
425: *        IMAX is the row-index of the largest off-diagonal element in
426: *        column K, and COLMAX is its absolute value
427: *
428:          IF( K.LT.N ) THEN
429:             IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 )
430:             COLMAX = CABS1( W( IMAX, K ) )
431:          ELSE
432:             COLMAX = ZERO
433:          END IF
434: *
435:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
436: *
437: *           Column K is zero: set INFO and continue
438: *
439:             IF( INFO.EQ.0 )
440:      $         INFO = K
441:             KP = K
442:             A( K, K ) = REAL( A( K, K ) )
443:          ELSE
444:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
445: *
446: *              no interchange, use 1-by-1 pivot block
447: *
448:                KP = K
449:             ELSE
450: *
451: *              Copy column IMAX to column K+1 of W and update it
452: *
453:                CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
454:                CALL CLACGV( IMAX-K, W( K, K+1 ), 1 )
455:                W( IMAX, K+1 ) = REAL( A( IMAX, IMAX ) )
456:                IF( IMAX.LT.N )
457:      $            CALL CCOPY( N-IMAX, A( IMAX+1, IMAX ), 1,
458:      $                        W( IMAX+1, K+1 ), 1 )
459:                CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
460:      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
461:      $                     1 )
462:                W( IMAX, K+1 ) = REAL( W( IMAX, K+1 ) )
463: *
464: *              JMAX is the column-index of the largest off-diagonal
465: *              element in row IMAX, and ROWMAX is its absolute value
466: *
467:                JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 )
468:                ROWMAX = CABS1( W( JMAX, K+1 ) )
469:                IF( IMAX.LT.N ) THEN
470:                   JMAX = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
471:                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
472:                END IF
473: *
474:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
475: *
476: *                 no interchange, use 1-by-1 pivot block
477: *
478:                   KP = K
479:                ELSE IF( ABS( REAL( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )
480:      $                   THEN
481: *
482: *                 interchange rows and columns K and IMAX, use 1-by-1
483: *                 pivot block
484: *
485:                   KP = IMAX
486: *
487: *                 copy column K+1 of W to column K
488: *
489:                   CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
490:                ELSE
491: *
492: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
493: *                 pivot block
494: *
495:                   KP = IMAX
496:                   KSTEP = 2
497:                END IF
498:             END IF
499: *
500:             KK = K + KSTEP - 1
501: *
502: *           Updated column KP is already stored in column KK of W
503: *
504:             IF( KP.NE.KK ) THEN
505: *
506: *              Copy non-updated column KK to column KP
507: *
508:                A( KP, KP ) = REAL( A( KK, KK ) )
509:                CALL CCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
510:      $                     LDA )
511:                CALL CLACGV( KP-KK-1, A( KP, KK+1 ), LDA )
512:                IF( KP.LT.N )
513:      $            CALL CCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
514: *
515: *              Interchange rows KK and KP in first KK columns of A and W
516: *
517:                CALL CSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
518:                CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
519:             END IF
520: *
521:             IF( KSTEP.EQ.1 ) THEN
522: *
523: *              1-by-1 pivot block D(k): column k of W now holds
524: *
525: *              W(k) = L(k)*D(k)
526: *
527: *              where L(k) is the k-th column of L
528: *
529: *              Store L(k) in column k of A
530: *
531:                CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
532:                IF( K.LT.N ) THEN
533:                   R1 = ONE / REAL( A( K, K ) )
534:                   CALL CSSCAL( N-K, R1, A( K+1, K ), 1 )
535: *
536: *                 Conjugate W(k)
537: *
538:                   CALL CLACGV( N-K, W( K+1, K ), 1 )
539:                END IF
540:             ELSE
541: *
542: *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
543: *
544: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
545: *
546: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
547: *              of L
548: *
549:                IF( K.LT.N-1 ) THEN
550: *
551: *                 Store L(k) and L(k+1) in columns k and k+1 of A
552: *
553:                   D21 = W( K+1, K )
554:                   D11 = W( K+1, K+1 ) / D21
555:                   D22 = W( K, K ) / CONJG( D21 )
556:                   T = ONE / ( REAL( D11*D22 )-ONE )
557:                   D21 = T / D21
558:                   DO 80 J = K + 2, N
559:                      A( J, K ) = CONJG( D21 )*
560:      $                           ( D11*W( J, K )-W( J, K+1 ) )
561:                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
562:    80             CONTINUE
563:                END IF
564: *
565: *              Copy D(k) to A
566: *
567:                A( K, K ) = W( K, K )
568:                A( K+1, K ) = W( K+1, K )
569:                A( K+1, K+1 ) = W( K+1, K+1 )
570: *
571: *              Conjugate W(k) and W(k+1)
572: *
573:                CALL CLACGV( N-K, W( K+1, K ), 1 )
574:                CALL CLACGV( N-K-1, W( K+2, K+1 ), 1 )
575:             END IF
576:          END IF
577: *
578: *        Store details of the interchanges in IPIV
579: *
580:          IF( KSTEP.EQ.1 ) THEN
581:             IPIV( K ) = KP
582:          ELSE
583:             IPIV( K ) = -KP
584:             IPIV( K+1 ) = -KP
585:          END IF
586: *
587: *        Increase K and return to the start of the main loop
588: *
589:          K = K + KSTEP
590:          GO TO 70
591: *
592:    90    CONTINUE
593: *
594: *        Update the lower triangle of A22 (= A(k:n,k:n)) as
595: *
596: *        A22 := A22 - L21*D*L21' = A22 - L21*W'
597: *
598: *        computing blocks of NB columns at a time (note that conjg(W) is
599: *        actually stored)
600: *
601:          DO 110 J = K, N, NB
602:             JB = MIN( NB, N-J+1 )
603: *
604: *           Update the lower triangle of the diagonal block
605: *
606:             DO 100 JJ = J, J + JB - 1
607:                A( JJ, JJ ) = REAL( A( JJ, JJ ) )
608:                CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
609:      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
610:      $                     A( JJ, JJ ), 1 )
611:                A( JJ, JJ ) = REAL( A( JJ, JJ ) )
612:   100       CONTINUE
613: *
614: *           Update the rectangular subdiagonal block
615: *
616:             IF( J+JB.LE.N )
617:      $         CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
618:      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
619:      $                     LDW, CONE, A( J+JB, J ), LDA )
620:   110    CONTINUE
621: *
622: *        Put L21 in standard form by partially undoing the interchanges
623: *        in columns 1:k-1
624: *
625:          J = K - 1
626:   120    CONTINUE
627:          JJ = J
628:          JP = IPIV( J )
629:          IF( JP.LT.0 ) THEN
630:             JP = -JP
631:             J = J - 1
632:          END IF
633:          J = J - 1
634:          IF( JP.NE.JJ .AND. J.GE.1 )
635:      $      CALL CSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
636:          IF( J.GE.1 )
637:      $      GO TO 120
638: *
639: *        Set KB to the number of columns factorized
640: *
641:          KB = K - 1
642: *
643:       END IF
644:       RETURN
645: *
646: *     End of CLAHEF
647: *
648:       END
649: