ssytrf()

subroutine ssytrf	(	character	uplo,
		integer	n,
		real, dimension( lda, * )	a,
		integer	lda,
		integer, dimension( * )	ipiv,
		real, dimension( * )	work,
		integer	lwork,
		integer	info )

SSYTRF

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

Purpose:

!>
!> SSYTRF computes the factorization of a real symmetric matrix A using
!> the Bunch-Kaufman diagonal pivoting method.  The form of the
!> factorization is
!>
!>    A = U**T*D*U  or  A = L*D*L**T
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, and D is symmetric and block diagonal with
!> 1-by-1 and 2-by-2 diagonal blocks.
!>
!> This is the blocked version of the algorithm, calling Level 3 BLAS.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> = 'U': Upper triangle of A is stored; !> = 'L': Lower triangle of A is stored. !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[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, the block diagonal matrix D and the multipliers used !> to obtain the factor U or L (see below for further details). !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,N). !>
[out]	IPIV	!> IPIV is INTEGER array, dimension (N) !> Details of the interchanges and the block structure of D. !> If IPIV(k) > 0, then rows and columns k and IPIV(k) were !> interchanged and D(k,k) is a 1-by-1 diagonal block. !> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and !> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) !> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = !> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were !> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. !>
[out]	WORK	!> WORK is REAL array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !>
[in]	LWORK	!> LWORK is INTEGER !> The length of WORK. LWORK >=1. For best performance !> LWORK >= N*NB, where NB is the block size returned by ILAENV. !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !> > 0: if INFO = i, D(i,i) is exactly zero. 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. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  If UPLO = 'U', then A = U**T*D*U, where
!>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
!>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
!>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    v    0   )   k-s
!>     U(k) =  (   0    I    0   )   s
!>             (   0    0    I   )   n-k
!>                k-s   s   n-k
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
!>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
!>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
!>
!>  If UPLO = 'L', then A = L*D*L**T, where
!>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
!>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
!>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    0     0   )  k-1
!>     L(k) =  (   0    I     0   )  s
!>             (   0    v     I   )  n-k-s+1
!>                k-1   s  n-k-s+1
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
!>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
!>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
!> 

Definition at line 179 of file ssytrf.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, LWORK, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               A( LDA, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            LQUERY, UPPER
      INTEGER            IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      REAL               SROUNDUP_LWORK
      EXTERNAL           lsame, ilaenv, sroundup_lwork
*     ..
*     .. External Subroutines ..
      EXTERNAL           slasyf, ssytf2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      upper = lsame( uplo, 'U' )
      lquery = ( lwork.EQ.-1 )
      IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -4
      ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
         info = -7
      END IF
*
      IF( info.EQ.0 ) THEN
*
*        Determine the block size
*
         nb = ilaenv( 1, 'SSYTRF', uplo, n, -1, -1, -1 )
         lwkopt = max( 1, n*nb )
         work( 1 ) = sroundup_lwork( lwkopt )
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SSYTRF', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
      nbmin = 2
      ldwork = n
      IF( nb.GT.1 .AND. nb.LT.n ) THEN
         iws = ldwork*nb
         IF( lwork.LT.iws ) THEN
            nb = max( lwork / ldwork, 1 )
            nbmin = max( 2, ilaenv( 2, 'SSYTRF', uplo, n, -1, -1,
     $                   -1 ) )
         END IF
      ELSE
         iws = 1
      END IF
      IF( nb.LT.nbmin )
     $   nb = n
*
      IF( upper ) THEN
*
*        Factorize A as U**T*D*U using the upper triangle of A
*
*        K is the main loop index, decreasing from N to 1 in steps of
*        KB, where KB is the number of columns factorized by SLASYF;
*        KB is either NB or NB-1, or K for the last block
*
         k = n
   10    CONTINUE
*
*        If K < 1, exit from loop
*
         IF( k.LT.1 )
     $      GO TO 40
*
         IF( k.GT.nb ) THEN
*
*           Factorize columns k-kb+1:k of A and use blocked code to
*           update columns 1:k-kb
*
            CALL slasyf( uplo, k, nb, kb, a, lda, ipiv, work, ldwork,
     $                   iinfo )
         ELSE
*
*           Use unblocked code to factorize columns 1:k of A
*
            CALL ssytf2( uplo, k, a, lda, ipiv, iinfo )
            kb = k
         END IF
*
*        Set INFO on the first occurrence of a zero pivot
*
         IF( info.EQ.0 .AND. iinfo.GT.0 )
     $      info = iinfo
*
*        Decrease K and return to the start of the main loop
*
         k = k - kb
         GO TO 10
*
      ELSE
*
*        Factorize A as L*D*L**T using the lower triangle of A
*
*        K is the main loop index, increasing from 1 to N in steps of
*        KB, where KB is the number of columns factorized by SLASYF;
*        KB is either NB or NB-1, or N-K+1 for the last block
*
         k = 1
   20    CONTINUE
*
*        If K > N, exit from loop
*
         IF( k.GT.n )
     $      GO TO 40
*
         IF( k.LE.n-nb ) THEN
*
*           Factorize columns k:k+kb-1 of A and use blocked code to
*           update columns k+kb:n
*
            CALL slasyf( uplo, n-k+1, nb, kb, a( k, k ), lda,
     $                   ipiv( k ),
     $                   work, ldwork, iinfo )
         ELSE
*
*           Use unblocked code to factorize columns k:n of A
*
            CALL ssytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ),
     $                   iinfo )
            kb = n - k + 1
         END IF
*
*        Set INFO on the first occurrence of a zero pivot
*
         IF( info.EQ.0 .AND. iinfo.GT.0 )
     $      info = iinfo + k - 1
*
*        Adjust IPIV
*
         DO 30 j = k, k + kb - 1
            IF( ipiv( j ).GT.0 ) THEN
               ipiv( j ) = ipiv( j ) + k - 1
            ELSE
               ipiv( j ) = ipiv( j ) - k + 1
            END IF
   30    CONTINUE
*
*        Increase K and return to the start of the main loop
*
         k = k + kb
         GO TO 20
*
      END IF
*
   40 CONTINUE
*
      work( 1 ) = sroundup_lwork( lwkopt )
      RETURN
*
*     End of SSYTRF
*

Here is the call graph for this function:

Here is the caller graph for this function: