LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
clasyf_rk.f
Go to the documentation of this file.
1 *> \brief \b CLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bunch-Kaufman (rook) 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_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLASYF_RK( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, KB, LDA, LDW, N, NB
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * )
30 * COMPLEX A( LDA, * ), E( * ), W( LDW, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *> CLASYF_RK computes a partial factorization of a complex symmetric
39 *> matrix A using the bounded Bunch-Kaufman (rook) diagonal
40 *> pivoting method. The partial 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 *>
51 *> CLASYF_RK is an auxiliary routine called by CSYTRF_RK. It uses
52 *> blocked code (calling Level 3 BLAS) to update the submatrix
53 *> A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] UPLO
60 *> \verbatim
61 *> UPLO is CHARACTER*1
62 *> Specifies whether the upper or lower triangular part of the
63 *> symmetric matrix A is stored:
64 *> = 'U': Upper triangular
65 *> = 'L': Lower triangular
66 *> \endverbatim
67 *>
68 *> \param[in] N
69 *> \verbatim
70 *> N is INTEGER
71 *> The order of the matrix A. N >= 0.
72 *> \endverbatim
73 *>
74 *> \param[in] NB
75 *> \verbatim
76 *> NB is INTEGER
77 *> The maximum number of columns of the matrix A that should be
78 *> factored. NB should be at least 2 to allow for 2-by-2 pivot
79 *> blocks.
80 *> \endverbatim
81 *>
82 *> \param[out] KB
83 *> \verbatim
84 *> KB is INTEGER
85 *> The number of columns of A that were actually factored.
86 *> KB is either NB-1 or NB, or N if N <= NB.
87 *> \endverbatim
88 *>
89 *> \param[in,out] A
90 *> \verbatim
91 *> A is COMPLEX array, dimension (LDA,N)
92 *> On entry, the symmetric matrix A.
93 *> If UPLO = 'U': the leading N-by-N upper triangular part
94 *> of A contains the upper triangular part of the matrix A,
95 *> and the strictly lower triangular part of A is not
96 *> referenced.
97 *>
98 *> If UPLO = 'L': the leading N-by-N lower triangular part
99 *> of A contains the lower triangular part of the matrix A,
100 *> and the strictly upper triangular part of A is not
101 *> referenced.
102 *>
103 *> On exit, contains:
104 *> a) ONLY diagonal elements of the symmetric block diagonal
105 *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
106 *> (superdiagonal (or subdiagonal) elements of D
107 *> are stored on exit in array E), and
108 *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
109 *> If UPLO = 'L': factor L in the subdiagonal part of A.
110 *> \endverbatim
111 *>
112 *> \param[in] LDA
113 *> \verbatim
114 *> LDA is INTEGER
115 *> The leading dimension of the array A. LDA >= max(1,N).
116 *> \endverbatim
117 *>
118 *> \param[out] E
119 *> \verbatim
120 *> E is COMPLEX array, dimension (N)
121 *> On exit, contains the superdiagonal (or subdiagonal)
122 *> elements of the symmetric block diagonal matrix D
123 *> with 1-by-1 or 2-by-2 diagonal blocks, where
124 *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
125 *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
126 *>
127 *> NOTE: For 1-by-1 diagonal block D(k), where
128 *> 1 <= k <= N, the element E(k) is set to 0 in both
129 *> UPLO = 'U' or UPLO = 'L' cases.
130 *> \endverbatim
131 *>
132 *> \param[out] IPIV
133 *> \verbatim
134 *> IPIV is INTEGER array, dimension (N)
135 *> IPIV describes the permutation matrix P in the factorization
136 *> of matrix A as follows. The absolute value of IPIV(k)
137 *> represents the index of row and column that were
138 *> interchanged with the k-th row and column. The value of UPLO
139 *> describes the order in which the interchanges were applied.
140 *> Also, the sign of IPIV represents the block structure of
141 *> the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
142 *> diagonal blocks which correspond to 1 or 2 interchanges
143 *> at each factorization step.
144 *>
145 *> If UPLO = 'U',
146 *> ( in factorization order, k decreases from N to 1 ):
147 *> a) A single positive entry IPIV(k) > 0 means:
148 *> D(k,k) is a 1-by-1 diagonal block.
149 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
150 *> interchanged in the submatrix A(1:N,N-KB+1:N);
151 *> If IPIV(k) = k, no interchange occurred.
152 *>
153 *>
154 *> b) A pair of consecutive negative entries
155 *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
156 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
157 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
158 *> 1) If -IPIV(k) != k, rows and columns
159 *> k and -IPIV(k) were interchanged
160 *> in the matrix A(1:N,N-KB+1:N).
161 *> If -IPIV(k) = k, no interchange occurred.
162 *> 2) If -IPIV(k-1) != k-1, rows and columns
163 *> k-1 and -IPIV(k-1) were interchanged
164 *> in the submatrix A(1:N,N-KB+1:N).
165 *> If -IPIV(k-1) = k-1, no interchange occurred.
166 *>
167 *> c) In both cases a) and b) is always ABS( IPIV(k) ) <= k.
168 *>
169 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
170 *>
171 *> If UPLO = 'L',
172 *> ( in factorization order, k increases from 1 to N ):
173 *> a) A single positive entry IPIV(k) > 0 means:
174 *> D(k,k) is a 1-by-1 diagonal block.
175 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
176 *> interchanged in the submatrix A(1:N,1:KB).
177 *> If IPIV(k) = k, no interchange occurred.
178 *>
179 *> b) A pair of consecutive negative entries
180 *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
181 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
182 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
183 *> 1) If -IPIV(k) != k, rows and columns
184 *> k and -IPIV(k) were interchanged
185 *> in the submatrix A(1:N,1:KB).
186 *> If -IPIV(k) = k, no interchange occurred.
187 *> 2) If -IPIV(k+1) != k+1, rows and columns
188 *> k-1 and -IPIV(k-1) were interchanged
189 *> in the submatrix A(1:N,1:KB).
190 *> If -IPIV(k+1) = k+1, no interchange occurred.
191 *>
192 *> c) In both cases a) and b) is always ABS( IPIV(k) ) >= k.
193 *>
194 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
195 *> \endverbatim
196 *>
197 *> \param[out] W
198 *> \verbatim
199 *> W is COMPLEX array, dimension (LDW,NB)
200 *> \endverbatim
201 *>
202 *> \param[in] LDW
203 *> \verbatim
204 *> LDW is INTEGER
205 *> The leading dimension of the array W. LDW >= max(1,N).
206 *> \endverbatim
207 *>
208 *> \param[out] INFO
209 *> \verbatim
210 *> INFO is INTEGER
211 *> = 0: successful exit
212 *>
213 *> < 0: If INFO = -k, the k-th argument had an illegal value
214 *>
215 *> > 0: If INFO = k, the matrix A is singular, because:
216 *> If UPLO = 'U': column k in the upper
217 *> triangular part of A contains all zeros.
218 *> If UPLO = 'L': column k in the lower
219 *> triangular part of A contains all zeros.
220 *>
221 *> Therefore D(k,k) is exactly zero, and superdiagonal
222 *> elements of column k of U (or subdiagonal elements of
223 *> column k of L ) are all zeros. The factorization has
224 *> been completed, but the block diagonal matrix D is
225 *> exactly singular, and division by zero will occur if
226 *> it is used to solve a system of equations.
227 *>
228 *> NOTE: INFO only stores the first occurrence of
229 *> a singularity, any subsequent occurrence of singularity
230 *> is not stored in INFO even though the factorization
231 *> always completes.
232 *> \endverbatim
233 *
234 * Authors:
235 * ========
236 *
237 *> \author Univ. of Tennessee
238 *> \author Univ. of California Berkeley
239 *> \author Univ. of Colorado Denver
240 *> \author NAG Ltd.
241 *
242 *> \date December 2016
243 *
244 *> \ingroup complexSYcomputational
245 *
246 *> \par Contributors:
247 * ==================
248 *>
249 *> \verbatim
250 *>
251 *> December 2016, Igor Kozachenko,
252 *> Computer Science Division,
253 *> University of California, Berkeley
254 *>
255 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
256 *> School of Mathematics,
257 *> University of Manchester
258 *>
259 *> \endverbatim
260 *
261 * =====================================================================
262  SUBROUTINE clasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
263  $ INFO )
264 *
265 * -- LAPACK computational routine (version 3.7.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * December 2016
269 *
270 * .. Scalar Arguments ..
271  CHARACTER UPLO
272  INTEGER INFO, KB, LDA, LDW, N, NB
273 * ..
274 * .. Array Arguments ..
275  INTEGER IPIV( * )
276  COMPLEX A( LDA, * ), E( * ), W( LDW, * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  REAL ZERO, ONE
283  parameter( zero = 0.0e+0, one = 1.0e+0 )
284  REAL EIGHT, SEVTEN
285  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
286  COMPLEX CONE, CZERO
287  parameter( cone = ( 1.0e+0, 0.0e+0 ),
288  $ czero = ( 0.0e+0, 0.0e+0 ) )
289 * ..
290 * .. Local Scalars ..
291  LOGICAL DONE
292  INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
293  $ kp, kstep, p, ii
294  REAL ABSAKK, ALPHA, COLMAX, ROWMAX, SFMIN, STEMP
295  COMPLEX D11, D12, D21, D22, R1, T, Z
296 * ..
297 * .. External Functions ..
298  LOGICAL LSAME
299  INTEGER ICAMAX
300  REAL SLAMCH
301  EXTERNAL lsame, icamax, slamch
302 * ..
303 * .. External Subroutines ..
304  EXTERNAL ccopy, cgemm, cgemv, cscal, cswap
305 * ..
306 * .. Intrinsic Functions ..
307  INTRINSIC abs, aimag, max, min, real, sqrt
308 * ..
309 * .. Statement Functions ..
310  REAL CABS1
311 * ..
312 * .. Statement Function definitions ..
313  cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
314 * ..
315 * .. Executable Statements ..
316 *
317  info = 0
318 *
319 * Initialize ALPHA for use in choosing pivot block size.
320 *
321  alpha = ( one+sqrt( sevten ) ) / eight
322 *
323 * Compute machine safe minimum
324 *
325  sfmin = slamch( 'S' )
326 *
327  IF( lsame( uplo, 'U' ) ) THEN
328 *
329 * Factorize the trailing columns of A using the upper triangle
330 * of A and working backwards, and compute the matrix W = U12*D
331 * for use in updating A11
332 *
333 * Initialize the first entry of array E, where superdiagonal
334 * elements of D are stored
335 *
336  e( 1 ) = czero
337 *
338 * K is the main loop index, decreasing from N in steps of 1 or 2
339 *
340  k = n
341  10 CONTINUE
342 *
343 * KW is the column of W which corresponds to column K of A
344 *
345  kw = nb + k - n
346 *
347 * Exit from loop
348 *
349  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
350  $ GO TO 30
351 *
352  kstep = 1
353  p = k
354 *
355 * Copy column K of A to column KW of W and update it
356 *
357  CALL ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
358  IF( k.LT.n )
359  $ CALL cgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
360  $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
361 *
362 * Determine rows and columns to be interchanged and whether
363 * a 1-by-1 or 2-by-2 pivot block will be used
364 *
365  absakk = cabs1( w( k, kw ) )
366 *
367 * IMAX is the row-index of the largest off-diagonal element in
368 * column K, and COLMAX is its absolute value.
369 * Determine both COLMAX and IMAX.
370 *
371  IF( k.GT.1 ) THEN
372  imax = icamax( k-1, w( 1, kw ), 1 )
373  colmax = cabs1( w( imax, kw ) )
374  ELSE
375  colmax = zero
376  END IF
377 *
378  IF( max( absakk, colmax ).EQ.zero ) THEN
379 *
380 * Column K is zero or underflow: set INFO and continue
381 *
382  IF( info.EQ.0 )
383  $ info = k
384  kp = k
385  CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
386 *
387 * Set E( K ) to zero
388 *
389  IF( k.GT.1 )
390  $ e( k ) = czero
391 *
392  ELSE
393 *
394 * ============================================================
395 *
396 * Test for interchange
397 *
398 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
399 * (used to handle NaN and Inf)
400 *
401  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
402 *
403 * no interchange, use 1-by-1 pivot block
404 *
405  kp = k
406 *
407  ELSE
408 *
409  done = .false.
410 *
411 * Loop until pivot found
412 *
413  12 CONTINUE
414 *
415 * Begin pivot search loop body
416 *
417 *
418 * Copy column IMAX to column KW-1 of W and update it
419 *
420  CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
421  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
422  $ w( imax+1, kw-1 ), 1 )
423 *
424  IF( k.LT.n )
425  $ CALL cgemv( 'No transpose', k, n-k, -cone,
426  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
427  $ cone, w( 1, kw-1 ), 1 )
428 *
429 * JMAX is the column-index of the largest off-diagonal
430 * element in row IMAX, and ROWMAX is its absolute value.
431 * Determine both ROWMAX and JMAX.
432 *
433  IF( imax.NE.k ) THEN
434  jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
435  $ 1 )
436  rowmax = cabs1( w( jmax, kw-1 ) )
437  ELSE
438  rowmax = zero
439  END IF
440 *
441  IF( imax.GT.1 ) THEN
442  itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
443  stemp = cabs1( w( itemp, kw-1 ) )
444  IF( stemp.GT.rowmax ) THEN
445  rowmax = stemp
446  jmax = itemp
447  END IF
448  END IF
449 *
450 * Equivalent to testing for
451 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
452 * (used to handle NaN and Inf)
453 *
454  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
455  $ THEN
456 *
457 * interchange rows and columns K and IMAX,
458 * use 1-by-1 pivot block
459 *
460  kp = imax
461 *
462 * copy column KW-1 of W to column KW of W
463 *
464  CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
465 *
466  done = .true.
467 *
468 * Equivalent to testing for ROWMAX.EQ.COLMAX,
469 * (used to handle NaN and Inf)
470 *
471  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
472  $ THEN
473 *
474 * interchange rows and columns K-1 and IMAX,
475 * use 2-by-2 pivot block
476 *
477  kp = imax
478  kstep = 2
479  done = .true.
480  ELSE
481 *
482 * Pivot not found: set params and repeat
483 *
484  p = imax
485  colmax = rowmax
486  imax = jmax
487 *
488 * Copy updated JMAXth (next IMAXth) column to Kth of W
489 *
490  CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
491 *
492  END IF
493 *
494 * End pivot search loop body
495 *
496  IF( .NOT. done ) GOTO 12
497 *
498  END IF
499 *
500 * ============================================================
501 *
502  kk = k - kstep + 1
503 *
504 * KKW is the column of W which corresponds to column KK of A
505 *
506  kkw = nb + kk - n
507 *
508  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
509 *
510 * Copy non-updated column K to column P
511 *
512  CALL ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
513  CALL ccopy( p, a( 1, k ), 1, a( 1, p ), 1 )
514 *
515 * Interchange rows K and P in last N-K+1 columns of A
516 * and last N-K+2 columns of W
517 *
518  CALL cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
519  CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
520  END IF
521 *
522 * Updated column KP is already stored in column KKW of W
523 *
524  IF( kp.NE.kk ) THEN
525 *
526 * Copy non-updated column KK to column KP
527 *
528  a( kp, k ) = a( kk, k )
529  CALL ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
530  $ lda )
531  CALL ccopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
532 *
533 * Interchange rows KK and KP in last N-KK+1 columns
534 * of A and W
535 *
536  CALL cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
537  CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
538  $ ldw )
539  END IF
540 *
541  IF( kstep.EQ.1 ) THEN
542 *
543 * 1-by-1 pivot block D(k): column KW of W now holds
544 *
545 * W(k) = U(k)*D(k)
546 *
547 * where U(k) is the k-th column of U
548 *
549 * Store U(k) in column k of A
550 *
551  CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
552  IF( k.GT.1 ) THEN
553  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
554  r1 = cone / a( k, k )
555  CALL cscal( k-1, r1, a( 1, k ), 1 )
556  ELSE IF( a( k, k ).NE.czero ) THEN
557  DO 14 ii = 1, k - 1
558  a( ii, k ) = a( ii, k ) / a( k, k )
559  14 CONTINUE
560  END IF
561 *
562 * Store the superdiagonal element of D in array E
563 *
564  e( k ) = czero
565 *
566  END IF
567 *
568  ELSE
569 *
570 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
571 * hold
572 *
573 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
574 *
575 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
576 * of U
577 *
578  IF( k.GT.2 ) THEN
579 *
580 * Store U(k) and U(k-1) in columns k and k-1 of A
581 *
582  d12 = w( k-1, kw )
583  d11 = w( k, kw ) / d12
584  d22 = w( k-1, kw-1 ) / d12
585  t = cone / ( d11*d22-cone )
586  DO 20 j = 1, k - 2
587  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
588  $ d12 )
589  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
590  $ d12 )
591  20 CONTINUE
592  END IF
593 *
594 * Copy diagonal elements of D(K) to A,
595 * copy superdiagonal element of D(K) to E(K) and
596 * ZERO out superdiagonal entry of A
597 *
598  a( k-1, k-1 ) = w( k-1, kw-1 )
599  a( k-1, k ) = czero
600  a( k, k ) = w( k, kw )
601  e( k ) = w( k-1, kw )
602  e( k-1 ) = czero
603 *
604  END IF
605 *
606 * End column K is nonsingular
607 *
608  END IF
609 *
610 * Store details of the interchanges in IPIV
611 *
612  IF( kstep.EQ.1 ) THEN
613  ipiv( k ) = kp
614  ELSE
615  ipiv( k ) = -p
616  ipiv( k-1 ) = -kp
617  END IF
618 *
619 * Decrease K and return to the start of the main loop
620 *
621  k = k - kstep
622  GO TO 10
623 *
624  30 CONTINUE
625 *
626 * Update the upper triangle of A11 (= A(1:k,1:k)) as
627 *
628 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
629 *
630 * computing blocks of NB columns at a time
631 *
632  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
633  jb = min( nb, k-j+1 )
634 *
635 * Update the upper triangle of the diagonal block
636 *
637  DO 40 jj = j, j + jb - 1
638  CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
639  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
640  $ a( j, jj ), 1 )
641  40 CONTINUE
642 *
643 * Update the rectangular superdiagonal block
644 *
645  IF( j.GE.2 )
646  $ CALL cgemm( 'No transpose', 'Transpose', j-1, jb,
647  $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ),
648  $ ldw, cone, a( 1, j ), lda )
649  50 CONTINUE
650 *
651 * Set KB to the number of columns factorized
652 *
653  kb = n - k
654 *
655  ELSE
656 *
657 * Factorize the leading columns of A using the lower triangle
658 * of A and working forwards, and compute the matrix W = L21*D
659 * for use in updating A22
660 *
661 * Initialize the unused last entry of the subdiagonal array E.
662 *
663  e( n ) = czero
664 *
665 * K is the main loop index, increasing from 1 in steps of 1 or 2
666 *
667  k = 1
668  70 CONTINUE
669 *
670 * Exit from loop
671 *
672  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
673  $ GO TO 90
674 *
675  kstep = 1
676  p = k
677 *
678 * Copy column K of A to column K of W and update it
679 *
680  CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
681  IF( k.GT.1 )
682  $ CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
683  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
684 *
685 * Determine rows and columns to be interchanged and whether
686 * a 1-by-1 or 2-by-2 pivot block will be used
687 *
688  absakk = cabs1( w( k, k ) )
689 *
690 * IMAX is the row-index of the largest off-diagonal element in
691 * column K, and COLMAX is its absolute value.
692 * Determine both COLMAX and IMAX.
693 *
694  IF( k.LT.n ) THEN
695  imax = k + icamax( n-k, w( k+1, k ), 1 )
696  colmax = cabs1( w( imax, k ) )
697  ELSE
698  colmax = zero
699  END IF
700 *
701  IF( max( absakk, colmax ).EQ.zero ) THEN
702 *
703 * Column K is zero or underflow: set INFO and continue
704 *
705  IF( info.EQ.0 )
706  $ info = k
707  kp = k
708  CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
709 *
710 * Set E( K ) to zero
711 *
712  IF( k.LT.n )
713  $ e( k ) = czero
714 *
715  ELSE
716 *
717 * ============================================================
718 *
719 * Test for interchange
720 *
721 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
722 * (used to handle NaN and Inf)
723 *
724  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
725 *
726 * no interchange, use 1-by-1 pivot block
727 *
728  kp = k
729 *
730  ELSE
731 *
732  done = .false.
733 *
734 * Loop until pivot found
735 *
736  72 CONTINUE
737 *
738 * Begin pivot search loop body
739 *
740 *
741 * Copy column IMAX to column K+1 of W and update it
742 *
743  CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
744  CALL ccopy( n-imax+1, a( imax, imax ), 1,
745  $ w( imax, k+1 ), 1 )
746  IF( k.GT.1 )
747  $ CALL cgemv( 'No transpose', n-k+1, k-1, -cone,
748  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
749  $ cone, w( k, k+1 ), 1 )
750 *
751 * JMAX is the column-index of the largest off-diagonal
752 * element in row IMAX, and ROWMAX is its absolute value.
753 * Determine both ROWMAX and JMAX.
754 *
755  IF( imax.NE.k ) THEN
756  jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
757  rowmax = cabs1( w( jmax, k+1 ) )
758  ELSE
759  rowmax = zero
760  END IF
761 *
762  IF( imax.LT.n ) THEN
763  itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
764  stemp = cabs1( w( itemp, k+1 ) )
765  IF( stemp.GT.rowmax ) THEN
766  rowmax = stemp
767  jmax = itemp
768  END IF
769  END IF
770 *
771 * Equivalent to testing for
772 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
773 * (used to handle NaN and Inf)
774 *
775  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
776  $ THEN
777 *
778 * interchange rows and columns K and IMAX,
779 * use 1-by-1 pivot block
780 *
781  kp = imax
782 *
783 * copy column K+1 of W to column K of W
784 *
785  CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
786 *
787  done = .true.
788 *
789 * Equivalent to testing for ROWMAX.EQ.COLMAX,
790 * (used to handle NaN and Inf)
791 *
792  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
793  $ THEN
794 *
795 * interchange rows and columns K+1 and IMAX,
796 * use 2-by-2 pivot block
797 *
798  kp = imax
799  kstep = 2
800  done = .true.
801  ELSE
802 *
803 * Pivot not found: set params and repeat
804 *
805  p = imax
806  colmax = rowmax
807  imax = jmax
808 *
809 * Copy updated JMAXth (next IMAXth) column to Kth of W
810 *
811  CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
812 *
813  END IF
814 *
815 * End pivot search loop body
816 *
817  IF( .NOT. done ) GOTO 72
818 *
819  END IF
820 *
821 * ============================================================
822 *
823  kk = k + kstep - 1
824 *
825  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
826 *
827 * Copy non-updated column K to column P
828 *
829  CALL ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
830  CALL ccopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
831 *
832 * Interchange rows K and P in first K columns of A
833 * and first K+1 columns of W
834 *
835  CALL cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
836  CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
837  END IF
838 *
839 * Updated column KP is already stored in column KK of W
840 *
841  IF( kp.NE.kk ) THEN
842 *
843 * Copy non-updated column KK to column KP
844 *
845  a( kp, k ) = a( kk, k )
846  CALL ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
847  CALL ccopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
848 *
849 * Interchange rows KK and KP in first KK columns of A and W
850 *
851  CALL cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
852  CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
853  END IF
854 *
855  IF( kstep.EQ.1 ) THEN
856 *
857 * 1-by-1 pivot block D(k): column k of W now holds
858 *
859 * W(k) = L(k)*D(k)
860 *
861 * where L(k) is the k-th column of L
862 *
863 * Store L(k) in column k of A
864 *
865  CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
866  IF( k.LT.n ) THEN
867  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
868  r1 = cone / a( k, k )
869  CALL cscal( n-k, r1, a( k+1, k ), 1 )
870  ELSE IF( a( k, k ).NE.czero ) THEN
871  DO 74 ii = k + 1, n
872  a( ii, k ) = a( ii, k ) / a( k, k )
873  74 CONTINUE
874  END IF
875 *
876 * Store the subdiagonal element of D in array E
877 *
878  e( k ) = czero
879 *
880  END IF
881 *
882  ELSE
883 *
884 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
885 *
886 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
887 *
888 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
889 * of L
890 *
891  IF( k.LT.n-1 ) THEN
892 *
893 * Store L(k) and L(k+1) in columns k and k+1 of A
894 *
895  d21 = w( k+1, k )
896  d11 = w( k+1, k+1 ) / d21
897  d22 = w( k, k ) / d21
898  t = cone / ( d11*d22-cone )
899  DO 80 j = k + 2, n
900  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
901  $ d21 )
902  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
903  $ d21 )
904  80 CONTINUE
905  END IF
906 *
907 * Copy diagonal elements of D(K) to A,
908 * copy subdiagonal element of D(K) to E(K) and
909 * ZERO out subdiagonal entry of A
910 *
911  a( k, k ) = w( k, k )
912  a( k+1, k ) = czero
913  a( k+1, k+1 ) = w( k+1, k+1 )
914  e( k ) = w( k+1, k )
915  e( k+1 ) = czero
916 *
917  END IF
918 *
919 * End column K is nonsingular
920 *
921  END IF
922 *
923 * Store details of the interchanges in IPIV
924 *
925  IF( kstep.EQ.1 ) THEN
926  ipiv( k ) = kp
927  ELSE
928  ipiv( k ) = -p
929  ipiv( k+1 ) = -kp
930  END IF
931 *
932 * Increase K and return to the start of the main loop
933 *
934  k = k + kstep
935  GO TO 70
936 *
937  90 CONTINUE
938 *
939 * Update the lower triangle of A22 (= A(k:n,k:n)) as
940 *
941 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
942 *
943 * computing blocks of NB columns at a time
944 *
945  DO 110 j = k, n, nb
946  jb = min( nb, n-j+1 )
947 *
948 * Update the lower triangle of the diagonal block
949 *
950  DO 100 jj = j, j + jb - 1
951  CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
952  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
953  $ a( jj, jj ), 1 )
954  100 CONTINUE
955 *
956 * Update the rectangular subdiagonal block
957 *
958  IF( j+jb.LE.n )
959  $ CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
960  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
961  $ ldw, cone, a( j+jb, j ), lda )
962  110 CONTINUE
963 *
964 * Set KB to the number of columns factorized
965 *
966  kb = k - 1
967 *
968  END IF
969 *
970  RETURN
971 *
972 * End of CLASYF_RK
973 *
974  END
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
clasyf_rk
subroutine clasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
CLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bun...
Definition: clasyf_rk.f:264
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83