LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlasyf_rook()

subroutine zlasyf_rook ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

Download ZLASYF_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZLASYF_ROOK computes a partial factorization of a complex symmetric
 matrix A using the bounded Bunch-Kaufman ("rook") diagonal
 pivoting method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )

 A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
       ( L21  I ) (  0  A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.

 ZLASYF_ROOK is an auxiliary routine called by ZSYTRF_ROOK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
          n-by-n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n-by-n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.
          On exit, A contains details of the partial factorization.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.

          If UPLO = 'U':
             Only the last KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
             interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k-1 and -IPIV(k-1) were inerchaged,
             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
             Only the first KB elements of IPIV are set.

             If IPIV(k) > 0, then rows and columns k and IPIV(k)
             were interchanged and D(k,k) is a 1-by-1 diagonal block.

             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
             columns k and -IPIV(k) were interchanged and rows and
             columns k+1 and -IPIV(k+1) were inerchaged,
             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]W
          W is COMPLEX*16 array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
  November 2013,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 182 of file zlasyf_rook.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
Here is the call graph for this function:
Here is the caller graph for this function: