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