LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ slasyf_rk()

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

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

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

Purpose:
 SLASYF_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.

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