LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zlasyf_rk.f
Go to the documentation of this file.
1 *> \brief \b ZLASYF_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 ZLASYF_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLASYF_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*16 A( LDA, * ), E( * ), W( LDW, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *> ZLASYF_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 *> ZLASYF_RK is an auxiliary routine called by ZSYTRF_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*16 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*16 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*16 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 *> \ingroup complex16SYcomputational
243 *
244 *> \par Contributors:
245 * ==================
246 *>
247 *> \verbatim
248 *>
249 *> December 2016, Igor Kozachenko,
250 *> Computer Science Division,
251 *> University of California, Berkeley
252 *>
253 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
254 *> School of Mathematics,
255 *> University of Manchester
256 *>
257 *> \endverbatim
258 *
259 * =====================================================================
260  SUBROUTINE zlasyf_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
261  $ INFO )
262 *
263 * -- LAPACK computational routine --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266 *
267 * .. Scalar Arguments ..
268  CHARACTER UPLO
269  INTEGER INFO, KB, LDA, LDW, N, NB
270 * ..
271 * .. Array Arguments ..
272  INTEGER IPIV( * )
273  COMPLEX*16 A( LDA, * ), E( * ), W( LDW, * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Parameters ..
279  DOUBLE PRECISION ZERO, ONE
280  parameter( zero = 0.0d+0, one = 1.0d+0 )
281  DOUBLE PRECISION EIGHT, SEVTEN
282  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
283  COMPLEX*16 CONE, CZERO
284  parameter( cone = ( 1.0d+0, 0.0d+0 ),
285  $ czero = ( 0.0d+0, 0.0d+0 ) )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL DONE
289  INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, K, KK, KW, KKW,
290  $ kp, kstep, p, ii
291  DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, SFMIN, DTEMP
292  COMPLEX*16 D11, D12, D21, D22, R1, T, Z
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME
296  INTEGER IZAMAX
297  DOUBLE PRECISION DLAMCH
298  EXTERNAL lsame, izamax, dlamch
299 * ..
300 * .. External Subroutines ..
301  EXTERNAL zcopy, zgemm, zgemv, zscal, zswap
302 * ..
303 * .. Intrinsic Functions ..
304  INTRINSIC abs, dble, dimag, max, min, sqrt
305 * ..
306 * .. Statement Functions ..
307  DOUBLE PRECISION CABS1
308 * ..
309 * .. Statement Function definitions ..
310  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
311 * ..
312 * .. Executable Statements ..
313 *
314  info = 0
315 *
316 * Initialize ALPHA for use in choosing pivot block size.
317 *
318  alpha = ( one+sqrt( sevten ) ) / eight
319 *
320 * Compute machine safe minimum
321 *
322  sfmin = dlamch( 'S' )
323 *
324  IF( lsame( uplo, 'U' ) ) THEN
325 *
326 * Factorize the trailing columns of A using the upper triangle
327 * of A and working backwards, and compute the matrix W = U12*D
328 * for use in updating A11
329 *
330 * Initialize the first entry of array E, where superdiagonal
331 * elements of D are stored
332 *
333  e( 1 ) = czero
334 *
335 * K is the main loop index, decreasing from N in steps of 1 or 2
336 *
337  k = n
338  10 CONTINUE
339 *
340 * KW is the column of W which corresponds to column K of A
341 *
342  kw = nb + k - n
343 *
344 * Exit from loop
345 *
346  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
347  $ GO TO 30
348 *
349  kstep = 1
350  p = k
351 *
352 * Copy column K of A to column KW of W and update it
353 *
354  CALL zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
355  IF( k.LT.n )
356  $ CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
357  $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
358 *
359 * Determine rows and columns to be interchanged and whether
360 * a 1-by-1 or 2-by-2 pivot block will be used
361 *
362  absakk = cabs1( w( k, kw ) )
363 *
364 * IMAX is the row-index of the largest off-diagonal element in
365 * column K, and COLMAX is its absolute value.
366 * Determine both COLMAX and IMAX.
367 *
368  IF( k.GT.1 ) THEN
369  imax = izamax( k-1, w( 1, kw ), 1 )
370  colmax = cabs1( w( imax, kw ) )
371  ELSE
372  colmax = zero
373  END IF
374 *
375  IF( max( absakk, colmax ).EQ.zero ) THEN
376 *
377 * Column K is zero or underflow: set INFO and continue
378 *
379  IF( info.EQ.0 )
380  $ info = k
381  kp = k
382  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
383 *
384 * Set E( K ) to zero
385 *
386  IF( k.GT.1 )
387  $ e( k ) = czero
388 *
389  ELSE
390 *
391 * ============================================================
392 *
393 * Test for interchange
394 *
395 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
396 * (used to handle NaN and Inf)
397 *
398  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
399 *
400 * no interchange, use 1-by-1 pivot block
401 *
402  kp = k
403 *
404  ELSE
405 *
406  done = .false.
407 *
408 * Loop until pivot found
409 *
410  12 CONTINUE
411 *
412 * Begin pivot search loop body
413 *
414 *
415 * Copy column IMAX to column KW-1 of W and update it
416 *
417  CALL zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
418  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
419  $ w( imax+1, kw-1 ), 1 )
420 *
421  IF( k.LT.n )
422  $ CALL zgemv( 'No transpose', k, n-k, -cone,
423  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
424  $ cone, w( 1, kw-1 ), 1 )
425 *
426 * JMAX is the column-index of the largest off-diagonal
427 * element in row IMAX, and ROWMAX is its absolute value.
428 * Determine both ROWMAX and JMAX.
429 *
430  IF( imax.NE.k ) THEN
431  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
432  $ 1 )
433  rowmax = cabs1( w( jmax, kw-1 ) )
434  ELSE
435  rowmax = zero
436  END IF
437 *
438  IF( imax.GT.1 ) THEN
439  itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
440  dtemp = cabs1( w( itemp, kw-1 ) )
441  IF( dtemp.GT.rowmax ) THEN
442  rowmax = dtemp
443  jmax = itemp
444  END IF
445  END IF
446 *
447 * Equivalent to testing for
448 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
449 * (used to handle NaN and Inf)
450 *
451  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
452  $ THEN
453 *
454 * interchange rows and columns K and IMAX,
455 * use 1-by-1 pivot block
456 *
457  kp = imax
458 *
459 * copy column KW-1 of W to column KW of W
460 *
461  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
462 *
463  done = .true.
464 *
465 * Equivalent to testing for ROWMAX.EQ.COLMAX,
466 * (used to handle NaN and Inf)
467 *
468  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
469  $ THEN
470 *
471 * interchange rows and columns K-1 and IMAX,
472 * use 2-by-2 pivot block
473 *
474  kp = imax
475  kstep = 2
476  done = .true.
477  ELSE
478 *
479 * Pivot not found: set params and repeat
480 *
481  p = imax
482  colmax = rowmax
483  imax = jmax
484 *
485 * Copy updated JMAXth (next IMAXth) column to Kth of W
486 *
487  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
488 *
489  END IF
490 *
491 * End pivot search loop body
492 *
493  IF( .NOT. done ) GOTO 12
494 *
495  END IF
496 *
497 * ============================================================
498 *
499  kk = k - kstep + 1
500 *
501 * KKW is the column of W which corresponds to column KK of A
502 *
503  kkw = nb + kk - n
504 *
505  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
506 *
507 * Copy non-updated column K to column P
508 *
509  CALL zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
510  CALL zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
511 *
512 * Interchange rows K and P in last N-K+1 columns of A
513 * and last N-K+2 columns of W
514 *
515  CALL zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
516  CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
517  END IF
518 *
519 * Updated column KP is already stored in column KKW of W
520 *
521  IF( kp.NE.kk ) THEN
522 *
523 * Copy non-updated column KK to column KP
524 *
525  a( kp, k ) = a( kk, k )
526  CALL zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
527  $ lda )
528  CALL zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
529 *
530 * Interchange rows KK and KP in last N-KK+1 columns
531 * of A and W
532 *
533  CALL zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
534  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
535  $ ldw )
536  END IF
537 *
538  IF( kstep.EQ.1 ) THEN
539 *
540 * 1-by-1 pivot block D(k): column KW of W now holds
541 *
542 * W(k) = U(k)*D(k)
543 *
544 * where U(k) is the k-th column of U
545 *
546 * Store U(k) in column k of A
547 *
548  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
549  IF( k.GT.1 ) THEN
550  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
551  r1 = cone / a( k, k )
552  CALL zscal( k-1, r1, a( 1, k ), 1 )
553  ELSE IF( a( k, k ).NE.czero ) THEN
554  DO 14 ii = 1, k - 1
555  a( ii, k ) = a( ii, k ) / a( k, k )
556  14 CONTINUE
557  END IF
558 *
559 * Store the superdiagonal element of D in array E
560 *
561  e( k ) = czero
562 *
563  END IF
564 *
565  ELSE
566 *
567 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
568 * hold
569 *
570 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
571 *
572 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
573 * of U
574 *
575  IF( k.GT.2 ) THEN
576 *
577 * Store U(k) and U(k-1) in columns k and k-1 of A
578 *
579  d12 = w( k-1, kw )
580  d11 = w( k, kw ) / d12
581  d22 = w( k-1, kw-1 ) / d12
582  t = cone / ( d11*d22-cone )
583  DO 20 j = 1, k - 2
584  a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
585  $ d12 )
586  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
587  $ d12 )
588  20 CONTINUE
589  END IF
590 *
591 * Copy diagonal elements of D(K) to A,
592 * copy superdiagonal element of D(K) to E(K) and
593 * ZERO out superdiagonal entry of A
594 *
595  a( k-1, k-1 ) = w( k-1, kw-1 )
596  a( k-1, k ) = czero
597  a( k, k ) = w( k, kw )
598  e( k ) = w( k-1, kw )
599  e( k-1 ) = czero
600 *
601  END IF
602 *
603 * End column K is nonsingular
604 *
605  END IF
606 *
607 * Store details of the interchanges in IPIV
608 *
609  IF( kstep.EQ.1 ) THEN
610  ipiv( k ) = kp
611  ELSE
612  ipiv( k ) = -p
613  ipiv( k-1 ) = -kp
614  END IF
615 *
616 * Decrease K and return to the start of the main loop
617 *
618  k = k - kstep
619  GO TO 10
620 *
621  30 CONTINUE
622 *
623 * Update the upper triangle of A11 (= A(1:k,1:k)) as
624 *
625 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
626 *
627 * computing blocks of NB columns at a time
628 *
629  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
630  jb = min( nb, k-j+1 )
631 *
632 * Update the upper triangle of the diagonal block
633 *
634  DO 40 jj = j, j + jb - 1
635  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
636  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
637  $ a( j, jj ), 1 )
638  40 CONTINUE
639 *
640 * Update the rectangular superdiagonal block
641 *
642  IF( j.GE.2 )
643  $ CALL zgemm( 'No transpose', 'Transpose', j-1, jb,
644  $ n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ),
645  $ ldw, cone, a( 1, j ), lda )
646  50 CONTINUE
647 *
648 * Set KB to the number of columns factorized
649 *
650  kb = n - k
651 *
652  ELSE
653 *
654 * Factorize the leading columns of A using the lower triangle
655 * of A and working forwards, and compute the matrix W = L21*D
656 * for use in updating A22
657 *
658 * Initialize the unused last entry of the subdiagonal array E.
659 *
660  e( n ) = czero
661 *
662 * K is the main loop index, increasing from 1 in steps of 1 or 2
663 *
664  k = 1
665  70 CONTINUE
666 *
667 * Exit from loop
668 *
669  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
670  $ GO TO 90
671 *
672  kstep = 1
673  p = k
674 *
675 * Copy column K of A to column K of W and update it
676 *
677  CALL zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
678  IF( k.GT.1 )
679  $ CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
680  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
681 *
682 * Determine rows and columns to be interchanged and whether
683 * a 1-by-1 or 2-by-2 pivot block will be used
684 *
685  absakk = cabs1( w( k, k ) )
686 *
687 * IMAX is the row-index of the largest off-diagonal element in
688 * column K, and COLMAX is its absolute value.
689 * Determine both COLMAX and IMAX.
690 *
691  IF( k.LT.n ) THEN
692  imax = k + izamax( n-k, w( k+1, k ), 1 )
693  colmax = cabs1( w( imax, k ) )
694  ELSE
695  colmax = zero
696  END IF
697 *
698  IF( max( absakk, colmax ).EQ.zero ) THEN
699 *
700 * Column K is zero or underflow: set INFO and continue
701 *
702  IF( info.EQ.0 )
703  $ info = k
704  kp = k
705  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
706 *
707 * Set E( K ) to zero
708 *
709  IF( k.LT.n )
710  $ e( k ) = czero
711 *
712  ELSE
713 *
714 * ============================================================
715 *
716 * Test for interchange
717 *
718 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
719 * (used to handle NaN and Inf)
720 *
721  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
722 *
723 * no interchange, use 1-by-1 pivot block
724 *
725  kp = k
726 *
727  ELSE
728 *
729  done = .false.
730 *
731 * Loop until pivot found
732 *
733  72 CONTINUE
734 *
735 * Begin pivot search loop body
736 *
737 *
738 * Copy column IMAX to column K+1 of W and update it
739 *
740  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
741  CALL zcopy( n-imax+1, a( imax, imax ), 1,
742  $ w( imax, k+1 ), 1 )
743  IF( k.GT.1 )
744  $ CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
745  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
746  $ cone, w( k, k+1 ), 1 )
747 *
748 * JMAX is the column-index of the largest off-diagonal
749 * element in row IMAX, and ROWMAX is its absolute value.
750 * Determine both ROWMAX and JMAX.
751 *
752  IF( imax.NE.k ) THEN
753  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
754  rowmax = cabs1( w( jmax, k+1 ) )
755  ELSE
756  rowmax = zero
757  END IF
758 *
759  IF( imax.LT.n ) THEN
760  itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
761  dtemp = cabs1( w( itemp, k+1 ) )
762  IF( dtemp.GT.rowmax ) THEN
763  rowmax = dtemp
764  jmax = itemp
765  END IF
766  END IF
767 *
768 * Equivalent to testing for
769 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
770 * (used to handle NaN and Inf)
771 *
772  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
773  $ THEN
774 *
775 * interchange rows and columns K and IMAX,
776 * use 1-by-1 pivot block
777 *
778  kp = imax
779 *
780 * copy column K+1 of W to column K of W
781 *
782  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
783 *
784  done = .true.
785 *
786 * Equivalent to testing for ROWMAX.EQ.COLMAX,
787 * (used to handle NaN and Inf)
788 *
789  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
790  $ THEN
791 *
792 * interchange rows and columns K+1 and IMAX,
793 * use 2-by-2 pivot block
794 *
795  kp = imax
796  kstep = 2
797  done = .true.
798  ELSE
799 *
800 * Pivot not found: set params and repeat
801 *
802  p = imax
803  colmax = rowmax
804  imax = jmax
805 *
806 * Copy updated JMAXth (next IMAXth) column to Kth of W
807 *
808  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
809 *
810  END IF
811 *
812 * End pivot search loop body
813 *
814  IF( .NOT. done ) GOTO 72
815 *
816  END IF
817 *
818 * ============================================================
819 *
820  kk = k + kstep - 1
821 *
822  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
823 *
824 * Copy non-updated column K to column P
825 *
826  CALL zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
827  CALL zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
828 *
829 * Interchange rows K and P in first K columns of A
830 * and first K+1 columns of W
831 *
832  CALL zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
833  CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
834  END IF
835 *
836 * Updated column KP is already stored in column KK of W
837 *
838  IF( kp.NE.kk ) THEN
839 *
840 * Copy non-updated column KK to column KP
841 *
842  a( kp, k ) = a( kk, k )
843  CALL zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
844  CALL zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
845 *
846 * Interchange rows KK and KP in first KK columns of A and W
847 *
848  CALL zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
849  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
850  END IF
851 *
852  IF( kstep.EQ.1 ) THEN
853 *
854 * 1-by-1 pivot block D(k): column k of W now holds
855 *
856 * W(k) = L(k)*D(k)
857 *
858 * where L(k) is the k-th column of L
859 *
860 * Store L(k) in column k of A
861 *
862  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
863  IF( k.LT.n ) THEN
864  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
865  r1 = cone / a( k, k )
866  CALL zscal( n-k, r1, a( k+1, k ), 1 )
867  ELSE IF( a( k, k ).NE.czero ) THEN
868  DO 74 ii = k + 1, n
869  a( ii, k ) = a( ii, k ) / a( k, k )
870  74 CONTINUE
871  END IF
872 *
873 * Store the subdiagonal element of D in array E
874 *
875  e( k ) = czero
876 *
877  END IF
878 *
879  ELSE
880 *
881 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
882 *
883 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
884 *
885 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
886 * of L
887 *
888  IF( k.LT.n-1 ) THEN
889 *
890 * Store L(k) and L(k+1) in columns k and k+1 of A
891 *
892  d21 = w( k+1, k )
893  d11 = w( k+1, k+1 ) / d21
894  d22 = w( k, k ) / d21
895  t = cone / ( d11*d22-cone )
896  DO 80 j = k + 2, n
897  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
898  $ d21 )
899  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
900  $ d21 )
901  80 CONTINUE
902  END IF
903 *
904 * Copy diagonal elements of D(K) to A,
905 * copy subdiagonal element of D(K) to E(K) and
906 * ZERO out subdiagonal entry of A
907 *
908  a( k, k ) = w( k, k )
909  a( k+1, k ) = czero
910  a( k+1, k+1 ) = w( k+1, k+1 )
911  e( k ) = w( k+1, k )
912  e( k+1 ) = czero
913 *
914  END IF
915 *
916 * End column K is nonsingular
917 *
918  END IF
919 *
920 * Store details of the interchanges in IPIV
921 *
922  IF( kstep.EQ.1 ) THEN
923  ipiv( k ) = kp
924  ELSE
925  ipiv( k ) = -p
926  ipiv( k+1 ) = -kp
927  END IF
928 *
929 * Increase K and return to the start of the main loop
930 *
931  k = k + kstep
932  GO TO 70
933 *
934  90 CONTINUE
935 *
936 * Update the lower triangle of A22 (= A(k:n,k:n)) as
937 *
938 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
939 *
940 * computing blocks of NB columns at a time
941 *
942  DO 110 j = k, n, nb
943  jb = min( nb, n-j+1 )
944 *
945 * Update the lower triangle of the diagonal block
946 *
947  DO 100 jj = j, j + jb - 1
948  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
949  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
950  $ a( jj, jj ), 1 )
951  100 CONTINUE
952 *
953 * Update the rectangular subdiagonal block
954 *
955  IF( j+jb.LE.n )
956  $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
957  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
958  $ ldw, cone, a( j+jb, j ), lda )
959  110 CONTINUE
960 *
961 * Set KB to the number of columns factorized
962 *
963  kb = k - 1
964 *
965  END IF
966 *
967  RETURN
968 *
969 * End of ZLASYF_RK
970 *
971  END
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zlasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
ZLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bun...
Definition: zlasyf_rk.f:262