OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
slasyf_rk.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine slasyf_rk (uplo, n, nb, kb, a, lda, e, ipiv, w, ldw, info)
 SLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-Kaufman (rook) diagonal pivoting method.

Function/Subroutine Documentation

◆ 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*
alpha
#define alpha
Definition eval.h:35
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21