Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / clasyf.f
1 *> \brief \b CLASYF computes a partial factorization of a complex symmetric matrix using the Bunch-Kaufman diagonal pivoting method.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLASYF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
22 *
23 *       .. Scalar Arguments ..
24 *       CHARACTER          UPLO
25 *       INTEGER            INFO, KB, LDA, LDW, N, NB
26 *       ..
27 *       .. Array Arguments ..
28 *       INTEGER            IPIV( * )
29 *       COMPLEX            A( LDA, * ), W( LDW, * )
30 *       ..
31 *
32 *
33 *> \par Purpose:
34 *  =============
35 *>
36 *> \verbatim
37 *>
38 *> CLASYF computes a partial factorization of a complex symmetric matrix
39 *> A using the Bunch-Kaufman diagonal pivoting method. The partial
40 *> factorization has the form:
41 *>
42 *> A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
43 *>       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
44 *>
45 *> A  =  ( L11  0 ) ( D    0  ) ( L11**T L21**T )  if UPLO = 'L'
46 *>       ( L21  I ) ( 0   A22 ) (  0       I    )
47 *>
48 *> where the order of D is at most NB. The actual order is returned in
49 *> the argument KB, and is either NB or NB-1, or N if N <= NB.
50 *> Note that U**T denotes the transpose of U.
51 *>
52 *> CLASYF is an auxiliary routine called by CSYTRF. It uses blocked code
53 *> (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
54 *> A22 (if UPLO = 'L').
55 *> \endverbatim
56 *
57 *  Arguments:
58 *  ==========
59 *
60 *> \param[in] UPLO
61 *> \verbatim
62 *>          UPLO is CHARACTER*1
63 *>          Specifies whether the upper or lower triangular part of the
64 *>          symmetric matrix A is stored:
65 *>          = 'U':  Upper triangular
66 *>          = 'L':  Lower triangular
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *>          N is INTEGER
72 *>          The order of the matrix A.  N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in] NB
76 *> \verbatim
77 *>          NB is INTEGER
78 *>          The maximum number of columns of the matrix A that should be
79 *>          factored.  NB should be at least 2 to allow for 2-by-2 pivot
80 *>          blocks.
81 *> \endverbatim
82 *>
83 *> \param[out] KB
84 *> \verbatim
85 *>          KB is INTEGER
86 *>          The number of columns of A that were actually factored.
87 *>          KB is either NB-1 or NB, or N if N <= NB.
88 *> \endverbatim
89 *>
90 *> \param[in,out] A
91 *> \verbatim
92 *>          A is COMPLEX array, dimension (LDA,N)
93 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
94 *>          n-by-n upper triangular part of A contains the upper
95 *>          triangular part of the matrix A, and the strictly lower
96 *>          triangular part of A is not referenced.  If UPLO = 'L', the
97 *>          leading n-by-n lower triangular part of A contains the lower
98 *>          triangular part of the matrix A, and the strictly upper
99 *>          triangular part of A is not referenced.
100 *>          On exit, A contains details of the partial factorization.
101 *> \endverbatim
102 *>
103 *> \param[in] LDA
104 *> \verbatim
105 *>          LDA is INTEGER
106 *>          The leading dimension of the array A.  LDA >= max(1,N).
107 *> \endverbatim
108 *>
109 *> \param[out] IPIV
110 *> \verbatim
111 *>          IPIV is INTEGER array, dimension (N)
112 *>          Details of the interchanges and the block structure of D.
113 *>
114 *>          If UPLO = 'U':
115 *>             Only the last KB elements of IPIV are set.
116 *>
117 *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
118 *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
119 *>
120 *>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
121 *>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
122 *>             is a 2-by-2 diagonal block.
123 *>
124 *>          If UPLO = 'L':
125 *>             Only the first KB elements of IPIV are set.
126 *>
127 *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
128 *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
129 *>
130 *>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
131 *>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
132 *>             is a 2-by-2 diagonal block.
133 *> \endverbatim
134 *>
135 *> \param[out] W
136 *> \verbatim
137 *>          W is COMPLEX array, dimension (LDW,NB)
138 *> \endverbatim
139 *>
140 *> \param[in] LDW
141 *> \verbatim
142 *>          LDW is INTEGER
143 *>          The leading dimension of the array W.  LDW >= max(1,N).
144 *> \endverbatim
145 *>
146 *> \param[out] INFO
147 *> \verbatim
148 *>          INFO is INTEGER
149 *>          = 0: successful exit
150 *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
151 *>               has been completed, but the block diagonal matrix D is
152 *>               exactly singular.
153 *> \endverbatim
154 *
155 *  Authors:
156 *  ========
157 *
158 *> \author Univ. of Tennessee
159 *> \author Univ. of California Berkeley
160 *> \author Univ. of Colorado Denver
161 *> \author NAG Ltd.
162 *
163 *> \date November 2013
164 *
165 *> \ingroup complexSYcomputational
166 *
167 *> \par Contributors:
168 *  ==================
169 *>
170 *> \verbatim
171 *>
172 *>  November 2013,  Igor Kozachenko,
173 *>                  Computer Science Division,
174 *>                  University of California, Berkeley
175 *> \endverbatim
176 *
177 *  =====================================================================
178       SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
179 *
180 *  -- LAPACK computational routine (version 3.5.0) --
181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 *     November 2013
184 *
185 *     .. Scalar Arguments ..
186       CHARACTER          UPLO
187       INTEGER            INFO, KB, LDA, LDW, N, NB
188 *     ..
189 *     .. Array Arguments ..
190       INTEGER            IPIV( * )
191       COMPLEX            A( LDA, * ), W( LDW, * )
192 *     ..
193 *
194 *  =====================================================================
195 *
196 *     .. Parameters ..
197       REAL               ZERO, ONE
198       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
199       REAL               EIGHT, SEVTEN
200       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
201       COMPLEX            CONE
202       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
203 *     ..
204 *     .. Local Scalars ..
205       INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
206      $                   KSTEP, KW
207       REAL               ABSAKK, ALPHA, COLMAX, ROWMAX
208       COMPLEX            D11, D21, D22, R1, T, Z
209 *     ..
210 *     .. External Functions ..
211       LOGICAL            LSAME
212       INTEGER            ICAMAX
213       EXTERNAL           LSAME, ICAMAX
214 *     ..
215 *     .. External Subroutines ..
216       EXTERNAL           CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
217 *     ..
218 *     .. Intrinsic Functions ..
219       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL, SQRT
220 *     ..
221 *     .. Statement Functions ..
222       REAL               CABS1
223 *     ..
224 *     .. Statement Function definitions ..
225       CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
226 *     ..
227 *     .. Executable Statements ..
228 *
229       INFO = 0
230 *
231 *     Initialize ALPHA for use in choosing pivot block size.
232 *
233       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
234 *
235       IF( LSAME( UPLO, 'U' ) ) THEN
236 *
237 *        Factorize the trailing columns of A using the upper triangle
238 *        of A and working backwards, and compute the matrix W = U12*D
239 *        for use in updating A11
240 *
241 *        K is the main loop index, decreasing from N in steps of 1 or 2
242 *
243 *        KW is the column of W which corresponds to column K of A
244 *
245          K = N
246    10    CONTINUE
247          KW = NB + K - N
248 *
249 *        Exit from loop
250 *
251          IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
252      $      GO TO 30
253 *
254 *        Copy column K of A to column KW of W and update it
255 *
256          CALL CCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
257          IF( K.LT.N )
258      $      CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
259      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
260 *
261          KSTEP = 1
262 *
263 *        Determine rows and columns to be interchanged and whether
264 *        a 1-by-1 or 2-by-2 pivot block will be used
265 *
266          ABSAKK = CABS1( W( K, KW ) )
267 *
268 *        IMAX is the row-index of the largest off-diagonal element in
269 *        column K, and COLMAX is its absolute value.
270 *        Determine both COLMAX and IMAX.
271 *
272          IF( K.GT.1 ) THEN
273             IMAX = ICAMAX( K-1, W( 1, KW ), 1 )
274             COLMAX = CABS1( W( IMAX, KW ) )
275          ELSE
276             COLMAX = ZERO
277          END IF
278 *
279          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
280 *
281 *           Column K is zero or underflow: set INFO and continue
282 *
283             IF( INFO.EQ.0 )
284      $         INFO = K
285             KP = K
286          ELSE
287             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
288 *
289 *              no interchange, use 1-by-1 pivot block
290 *
291                KP = K
292             ELSE
293 *
294 *              Copy column IMAX to column KW-1 of W and update it
295 *
296                CALL CCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
297                CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
298      $                     W( IMAX+1, KW-1 ), 1 )
299                IF( K.LT.N )
300      $            CALL CGEMV( 'No transpose', K, N-K, -CONE,
301      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
302      $                        CONE, W( 1, KW-1 ), 1 )
303 *
304 *              JMAX is the column-index of the largest off-diagonal
305 *              element in row IMAX, and ROWMAX is its absolute value
306 *
307                JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
308                ROWMAX = CABS1( W( JMAX, KW-1 ) )
309                IF( IMAX.GT.1 ) THEN
310                   JMAX = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 )
311                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
312                END IF
313 *
314                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
315 *
316 *                 no interchange, use 1-by-1 pivot block
317 *
318                   KP = K
319                ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
320 *
321 *                 interchange rows and columns K and IMAX, use 1-by-1
322 *                 pivot block
323 *
324                   KP = IMAX
325 *
326 *                 copy column KW-1 of W to column KW of W
327 *
328                   CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
329                ELSE
330 *
331 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
332 *                 pivot block
333 *
334                   KP = IMAX
335                   KSTEP = 2
336                END IF
337             END IF
338 *
339 *           ============================================================
340 *
341 *           KK is the column of A where pivoting step stopped
342 *
343             KK = K - KSTEP + 1
344 *
345 *           KKW is the column of W which corresponds to column KK of A
346 *
347             KKW = NB + KK - N
348 *
349 *           Interchange rows and columns KP and KK.
350 *           Updated column KP is already stored in column KKW of W.
351 *
352             IF( KP.NE.KK ) THEN
353 *
354 *              Copy non-updated column KK to column KP of submatrix A
355 *              at step K. No need to copy element into column K
356 *              (or K and K-1 for 2-by-2 pivot) of A, since these columns
357 *              will be later overwritten.
358 *
359                A( KP, KP ) = A( KK, KK )
360                CALL CCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
361      $                     LDA )
362                IF( KP.GT.1 )
363      $            CALL CCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
364 *
365 *              Interchange rows KK and KP in last K+1 to N columns of A
366 *              (columns K (or K and K-1 for 2-by-2 pivot) of A will be
367 *              later overwritten). Interchange rows KK and KP
368 *              in last KKW to NB columns of W.
369 *
370                IF( K.LT.N )
371      $            CALL CSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
372      $                        LDA )
373                CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
374      $                     LDW )
375             END IF
376 *
377             IF( KSTEP.EQ.1 ) THEN
378 *
379 *              1-by-1 pivot block D(k): column kw of W now holds
380 *
381 *              W(kw) = U(k)*D(k),
382 *
383 *              where U(k) is the k-th column of U
384 *
385 *              Store subdiag. elements of column U(k)
386 *              and 1-by-1 block D(k) in column k of A.
387 *              NOTE: Diagonal element U(k,k) is a UNIT element
388 *              and not stored.
389 *                 A(k,k) := D(k,k) = W(k,kw)
390 *                 A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
391 *
392                CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
393                R1 = CONE / A( K, K )
394                CALL CSCAL( K-1, R1, A( 1, K ), 1 )
395 *
396             ELSE
397 *
398 *              2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
399 *
400 *              ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
401 *
402 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
403 *              of U
404 *
405 *              Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
406 *              block D(k-1:k,k-1:k) in columns k-1 and k of A.
407 *              NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
408 *              block and not stored.
409 *                 A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
410 *                 A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
411 *                 = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
412 *
413                IF( K.GT.2 ) THEN
414 *
415 *                 Compose the columns of the inverse of 2-by-2 pivot
416 *                 block D in the following way to reduce the number
417 *                 of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
418 *                 this inverse
419 *
420 *                 D**(-1) = ( d11 d21 )**(-1) =
421 *                           ( d21 d22 )
422 *
423 *                 = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
424 *                                        ( (-d21 ) ( d11 ) )
425 *
426 *                 = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
427 *
428 *                   * ( ( d22/d21 ) (      -1 ) ) =
429 *                     ( (      -1 ) ( d11/d21 ) )
430 *
431 *                 = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) (  -1 ) ) =
432 *                                           ( ( -1  ) ( D22 ) )
433 *
434 *                 = 1/d21 * T * ( ( D11 ) (  -1 ) )
435 *                               ( (  -1 ) ( D22 ) )
436 *
437 *                 = D21 * ( ( D11 ) (  -1 ) )
438 *                         ( (  -1 ) ( D22 ) )
439 *
440                   D21 = W( K-1, KW )
441                   D11 = W( K, KW ) / D21
442                   D22 = W( K-1, KW-1 ) / D21
443                   T = CONE / ( D11*D22-CONE )
444 *
445 *                 Update elements in columns A(k-1) and A(k) as
446 *                 dot products of rows of ( W(kw-1) W(kw) ) and columns
447 *                 of D**(-1)
448 *
449                   D21 = T / D21
450                   DO 20 J = 1, K - 2
451                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
452                      A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
453    20             CONTINUE
454                END IF
455 *
456 *              Copy D(k) to A
457 *
458                A( K-1, K-1 ) = W( K-1, KW-1 )
459                A( K-1, K ) = W( K-1, KW )
460                A( K, K ) = W( K, KW )
461 *
462             END IF
463 *
464          END IF
465 *
466 *        Store details of the interchanges in IPIV
467 *
468          IF( KSTEP.EQ.1 ) THEN
469             IPIV( K ) = KP
470          ELSE
471             IPIV( K ) = -KP
472             IPIV( K-1 ) = -KP
473          END IF
474 *
475 *        Decrease K and return to the start of the main loop
476 *
477          K = K - KSTEP
478          GO TO 10
479 *
480    30    CONTINUE
481 *
482 *        Update the upper triangle of A11 (= A(1:k,1:k)) as
483 *
484 *        A11 := A11 - U12*D*U12**T = A11 - U12*W**T
485 *
486 *        computing blocks of NB columns at a time
487 *
488          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
489             JB = MIN( NB, K-J+1 )
490 *
491 *           Update the upper triangle of the diagonal block
492 *
493             DO 40 JJ = J, J + JB - 1
494                CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
495      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
496      $                     A( J, JJ ), 1 )
497    40       CONTINUE
498 *
499 *           Update the rectangular superdiagonal block
500 *
501             CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
502      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
503      $                  CONE, A( 1, J ), LDA )
504    50    CONTINUE
505 *
506 *        Put U12 in standard form by partially undoing the interchanges
507 *        in columns k+1:n looping backwards from k+1 to n
508 *
509          J = K + 1
510    60    CONTINUE
511 *
512 *           Undo the interchanges (if any) of rows JJ and JP at each
513 *           step J
514 *
515 *           (Here, J is a diagonal index)
516             JJ = J
517             JP = IPIV( J )
518             IF( JP.LT.0 ) THEN
519                JP = -JP
520 *              (Here, J is a diagonal index)
521                J = J + 1
522             END IF
523 *           (NOTE: Here, J is used to determine row length. Length N-J+1
524 *           of the rows to swap back doesn't include diagonal element)
525             J = J + 1
526             IF( JP.NE.JJ .AND. J.LE.N )
527      $         CALL CSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
528          IF( J.LT.N )
529      $      GO TO 60
530 *
531 *        Set KB to the number of columns factorized
532 *
533          KB = N - K
534 *
535       ELSE
536 *
537 *        Factorize the leading columns of A using the lower triangle
538 *        of A and working forwards, and compute the matrix W = L21*D
539 *        for use in updating A22
540 *
541 *        K is the main loop index, increasing from 1 in steps of 1 or 2
542 *
543          K = 1
544    70    CONTINUE
545 *
546 *        Exit from loop
547 *
548          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
549      $      GO TO 90
550 *
551 *        Copy column K of A to column K of W and update it
552 *
553          CALL CCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
554          CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
555      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
556 *
557          KSTEP = 1
558 *
559 *        Determine rows and columns to be interchanged and whether
560 *        a 1-by-1 or 2-by-2 pivot block will be used
561 *
562          ABSAKK = CABS1( W( K, K ) )
563 *
564 *        IMAX is the row-index of the largest off-diagonal element in
565 *        column K, and COLMAX is its absolute value.
566 *        Determine both COLMAX and IMAX.
567 *
568          IF( K.LT.N ) THEN
569             IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 )
570             COLMAX = CABS1( W( IMAX, K ) )
571          ELSE
572             COLMAX = ZERO
573          END IF
574 *
575          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
576 *
577 *           Column K is zero or underflow: set INFO and continue
578 *
579             IF( INFO.EQ.0 )
580      $         INFO = K
581             KP = K
582          ELSE
583             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
584 *
585 *              no interchange, use 1-by-1 pivot block
586 *
587                KP = K
588             ELSE
589 *
590 *              Copy column IMAX to column K+1 of W and update it
591 *
592                CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
593                CALL CCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
594      $                     1 )
595                CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
596      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
597      $                     1 )
598 *
599 *              JMAX is the column-index of the largest off-diagonal
600 *              element in row IMAX, and ROWMAX is its absolute value
601 *
602                JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 )
603                ROWMAX = CABS1( W( JMAX, K+1 ) )
604                IF( IMAX.LT.N ) THEN
605                   JMAX = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
606                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
607                END IF
608 *
609                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
610 *
611 *                 no interchange, use 1-by-1 pivot block
612 *
613                   KP = K
614                ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
615 *
616 *                 interchange rows and columns K and IMAX, use 1-by-1
617 *                 pivot block
618 *
619                   KP = IMAX
620 *
621 *                 copy column K+1 of W to column K of W
622 *
623                   CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
624                ELSE
625 *
626 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
627 *                 pivot block
628 *
629                   KP = IMAX
630                   KSTEP = 2
631                END IF
632             END IF
633 *
634 *           ============================================================
635 *
636 *           KK is the column of A where pivoting step stopped
637 *
638             KK = K + KSTEP - 1
639 *
640 *           Interchange rows and columns KP and KK.
641 *           Updated column KP is already stored in column KK of W.
642 *
643             IF( KP.NE.KK ) THEN
644 *
645 *              Copy non-updated column KK to column KP of submatrix A
646 *              at step K. No need to copy element into column K
647 *              (or K and K+1 for 2-by-2 pivot) of A, since these columns
648 *              will be later overwritten.
649 *
650                A( KP, KP ) = A( KK, KK )
651                CALL CCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
652      $                     LDA )
653                IF( KP.LT.N )
654      $            CALL CCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
655 *
656 *              Interchange rows KK and KP in first K-1 columns of A
657 *              (columns K (or K and K+1 for 2-by-2 pivot) of A will be
658 *              later overwritten). Interchange rows KK and KP
659 *              in first KK columns of W.
660 *
661                IF( K.GT.1 )
662      $            CALL CSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
663                CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
664             END IF
665 *
666             IF( KSTEP.EQ.1 ) THEN
667 *
668 *              1-by-1 pivot block D(k): column k of W now holds
669 *
670 *              W(k) = L(k)*D(k),
671 *
672 *              where L(k) is the k-th column of L
673 *
674 *              Store subdiag. elements of column L(k)
675 *              and 1-by-1 block D(k) in column k of A.
676 *              (NOTE: Diagonal element L(k,k) is a UNIT element
677 *              and not stored)
678 *                 A(k,k) := D(k,k) = W(k,k)
679 *                 A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
680 *
681                CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
682                IF( K.LT.N ) THEN
683                   R1 = CONE / A( K, K )
684                   CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
685                END IF
686 *
687             ELSE
688 *
689 *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
690 *
691 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
692 *
693 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
694 *              of L
695 *
696 *              Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
697 *              block D(k:k+1,k:k+1) in columns k and k+1 of A.
698 *              (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
699 *              block and not stored)
700 *                 A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
701 *                 A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
702 *                 = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
703 *
704                IF( K.LT.N-1 ) THEN
705 *
706 *                 Compose the columns of the inverse of 2-by-2 pivot
707 *                 block D in the following way to reduce the number
708 *                 of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
709 *                 this inverse
710 *
711 *                 D**(-1) = ( d11 d21 )**(-1) =
712 *                           ( d21 d22 )
713 *
714 *                 = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
715 *                                        ( (-d21 ) ( d11 ) )
716 *
717 *                 = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
718 *
719 *                   * ( ( d22/d21 ) (      -1 ) ) =
720 *                     ( (      -1 ) ( d11/d21 ) )
721 *
722 *                 = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) (  -1 ) ) =
723 *                                           ( ( -1  ) ( D22 ) )
724 *
725 *                 = 1/d21 * T * ( ( D11 ) (  -1 ) )
726 *                               ( (  -1 ) ( D22 ) )
727 *
728 *                 = D21 * ( ( D11 ) (  -1 ) )
729 *                         ( (  -1 ) ( D22 ) )
730 *
731                   D21 = W( K+1, K )
732                   D11 = W( K+1, K+1 ) / D21
733                   D22 = W( K, K ) / D21
734                   T = CONE / ( D11*D22-CONE )
735                   D21 = T / D21
736 *
737 *                 Update elements in columns A(k) and A(k+1) as
738 *                 dot products of rows of ( W(k) W(k+1) ) and columns
739 *                 of D**(-1)
740 *
741                   DO 80 J = K + 2, N
742                      A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
743                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
744    80             CONTINUE
745                END IF
746 *
747 *              Copy D(k) to A
748 *
749                A( K, K ) = W( K, K )
750                A( K+1, K ) = W( K+1, K )
751                A( K+1, K+1 ) = W( K+1, K+1 )
752 *
753             END IF
754 *
755          END IF
756 *
757 *        Store details of the interchanges in IPIV
758 *
759          IF( KSTEP.EQ.1 ) THEN
760             IPIV( K ) = KP
761          ELSE
762             IPIV( K ) = -KP
763             IPIV( K+1 ) = -KP
764          END IF
765 *
766 *        Increase K and return to the start of the main loop
767 *
768          K = K + KSTEP
769          GO TO 70
770 *
771    90    CONTINUE
772 *
773 *        Update the lower triangle of A22 (= A(k:n,k:n)) as
774 *
775 *        A22 := A22 - L21*D*L21**T = A22 - L21*W**T
776 *
777 *        computing blocks of NB columns at a time
778 *
779          DO 110 J = K, N, NB
780             JB = MIN( NB, N-J+1 )
781 *
782 *           Update the lower triangle of the diagonal block
783 *
784             DO 100 JJ = J, J + JB - 1
785                CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
786      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
787      $                     A( JJ, JJ ), 1 )
788   100       CONTINUE
789 *
790 *           Update the rectangular subdiagonal block
791 *
792             IF( J+JB.LE.N )
793      $         CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
794      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
795      $                     LDW, CONE, A( J+JB, J ), LDA )
796   110    CONTINUE
797 *
798 *        Put L21 in standard form by partially undoing the interchanges
799 *        of rows in columns 1:k-1 looping backwards from k-1 to 1
800 *
801          J = K - 1
802   120    CONTINUE
803 *
804 *           Undo the interchanges (if any) of rows JJ and JP at each
805 *           step J
806 *
807 *           (Here, J is a diagonal index)
808             JJ = J
809             JP = IPIV( J )
810             IF( JP.LT.0 ) THEN
811                JP = -JP
812 *              (Here, J is a diagonal index)
813                J = J - 1
814             END IF
815 *           (NOTE: Here, J is used to determine row length. Length J
816 *           of the rows to swap back doesn't include diagonal element)
817             J = J - 1
818             IF( JP.NE.JJ .AND. J.GE.1 )
819      $         CALL CSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
820          IF( J.GT.1 )
821      $      GO TO 120
822 *
823 *        Set KB to the number of columns factorized
824 *
825          KB = K - 1
826 *
827       END IF
828       RETURN
829 *
830 *     End of CLASYF
831 *
832       END