LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlasyf_rk()

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

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

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

Purpose:
 ZLASYF_RK 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_RK is an auxiliary routine called by ZSYTRF_RK. 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, contains:
            a) ONLY diagonal elements of the symmetric block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                are stored on exit in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]E
          E is COMPLEX*16 array, dimension (N)
          On exit, contains the superdiagonal (or subdiagonal)
          elements of the symmetric block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is set to 0 in both
          UPLO = 'U' or UPLO = 'L' cases.
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          IPIV describes the permutation matrix P in the factorization
          of matrix A as follows. The absolute value of IPIV(k)
          represents the index of row and column that were
          interchanged with the k-th row and column. The value of UPLO
          describes the order in which the interchanges were applied.
          Also, the sign of IPIV represents the block structure of
          the symmetric block diagonal matrix D with 1-by-1 or 2-by-2
          diagonal blocks which correspond to 1 or 2 interchanges
          at each factorization step.

          If UPLO = 'U',
          ( in factorization order, k decreases from N to 1 ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the submatrix A(1:N,N-KB+1:N);
               If IPIV(k) = k, no interchange occurred.


            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k-1) < 0 means:
               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the matrix A(1:N,N-KB+1:N).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k-1) != k-1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the submatrix A(1:N,N-KB+1:N).
                  If -IPIV(k-1) = k-1, no interchange occurred.

            c) In both cases a) and b) is always ABS( IPIV(k) ) <= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.

          If UPLO = 'L',
          ( in factorization order, k increases from 1 to N ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the submatrix A(1:N,1:KB).
               If IPIV(k) = k, no interchange occurred.

            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k+1) < 0 means:
               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the submatrix A(1:N,1:KB).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k+1) != k+1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the submatrix A(1:N,1:KB).
                  If -IPIV(k+1) = k+1, no interchange occurred.

            c) In both cases a) and b) is always ABS( IPIV(k) ) >= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
[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, the k-th argument had an illegal value

          > 0: If INFO = k, the matrix A is singular, because:
                 If UPLO = 'U': column k in the upper
                 triangular part of A contains all zeros.
                 If UPLO = 'L': column k in the lower
                 triangular part of A contains all zeros.

               Therefore D(k,k) is exactly zero, and superdiagonal
               elements of column k of U (or subdiagonal elements of
               column k of L ) are all zeros. The factorization has
               been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if
               it is used to solve a system of equations.

               NOTE: INFO only stores the first occurrence of
               a singularity, any subsequent occurrence of singularity
               is not stored in INFO even though the factorization
               always completes.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
  December 2016,  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 260 of file zlasyf_rk.f.

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 *
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: