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