LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlasyf_rk()

subroutine dlasyf_rk ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  E,
integer, dimension( * )  IPIV,
double precision, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

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

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

Purpose:
 DLASYF_RK computes a partial factorization of a real 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.

 DLASYF_RK is an auxiliary routine called by DSYTRF_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.
Date
December 2016
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 264 of file dlasyf_rk.f.

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