real

Functions
subroutine	sgtcon (norm, n, dl, d, du, du2, ipiv, anorm, rcond, work, iwork, info)
	SGTCON
subroutine	sgtrfs (trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
	SGTRFS
subroutine	sgttrf (n, dl, d, du, du2, ipiv, info)
	SGTTRF
subroutine	sgttrs (trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
	SGTTRS
subroutine	sgtts2 (itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb)
	SGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.

Detailed Description

This is the group of real computational functions for GT matrices

Function Documentation

sgtcon()

subroutine sgtcon	(	character	norm,
		integer	n,
		real, dimension( * )	dl,
		real, dimension( * )	d,
		real, dimension( * )	du,
		real, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		real	anorm,
		real	rcond,
		real, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

SGTCON

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

Purpose:

!>
!> SGTCON estimates the reciprocal of the condition number of a real
!> tridiagonal matrix A using the LU factorization as computed by
!> SGTTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
!> 

Parameters

[in]	NORM	!> NORM is CHARACTER*1 !> Specifies whether the 1-norm condition number or the !> infinity-norm condition number is required: !> = '1' or 'O': 1-norm; !> = 'I': Infinity-norm. !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in]	DL	!> DL is REAL array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A as computed by SGTTRF. !>
[in]	D	!> D is REAL array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is REAL array, dimension (N-1) !> The (n-1) elements of the first superdiagonal of U. !>
[in]	DU2	!> DU2 is REAL array, dimension (N-2) !> The (n-2) elements of the second superdiagonal of U. !>
[in]	IPIV	!> IPIV is INTEGER array, dimension (N) !> The pivot indices; for 1 <= i <= n, row i of the matrix was !> interchanged with row IPIV(i). IPIV(i) will always be either !> i or i+1; IPIV(i) = i indicates a row interchange was not !> required. !>
[in]	ANORM	!> ANORM is REAL !> If NORM = '1' or 'O', the 1-norm of the original matrix A. !> If NORM = 'I', the infinity-norm of the original matrix A. !>
[out]	RCOND	!> RCOND is REAL !> The reciprocal of the condition number of the matrix A, !> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an !> estimate of the 1-norm of inv(A) computed in this routine. !>
[out]	WORK	!> WORK is REAL array, dimension (2*N) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (N) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 144 of file sgtcon.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          NORM
      INTEGER            INFO, N
      REAL               ANORM, RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IWORK( * )
      REAL               D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      parameter( one = 1.0e+0, zero = 0.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ONENRM
      INTEGER            I, KASE, KASE1
      REAL               AINVNM
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgttrs, slacn2, xerbla
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments.
*
      info = 0
      onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
      IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( anorm.LT.zero ) THEN
         info = -8
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SGTCON', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      rcond = zero
      IF( n.EQ.0 ) THEN
         rcond = one
         RETURN
      ELSE IF( anorm.EQ.zero ) THEN
         RETURN
      END IF
*
*     Check that D(1:N) is non-zero.
*
      DO 10 i = 1, n
         IF( d( i ).EQ.zero )
     $      RETURN
   10 CONTINUE
*
      ainvnm = zero
      IF( onenrm ) THEN
         kase1 = 1
      ELSE
         kase1 = 2
      END IF
      kase = 0
   20 CONTINUE
      CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
      IF( kase.NE.0 ) THEN
         IF( kase.EQ.kase1 ) THEN
*
*           Multiply by inv(U)*inv(L).
*
            CALL sgttrs( 'No transpose', n, 1, dl, d, du, du2, ipiv,
     $                   work, n, info )
         ELSE
*
*           Multiply by inv(L**T)*inv(U**T).
*
            CALL sgttrs( 'Transpose', n, 1, dl, d, du, du2, ipiv, work,
     $                   n, info )
         END IF
         GO TO 20
      END IF
*
*     Compute the estimate of the reciprocal condition number.
*
      IF( ainvnm.NE.zero )
     $   rcond = ( one / ainvnm ) / anorm
*
      RETURN
*
*     End of SGTCON
*

sgtrfs()

subroutine sgtrfs	(	character	trans,
		integer	n,
		integer	nrhs,
		real, dimension( * )	dl,
		real, dimension( * )	d,
		real, dimension( * )	du,
		real, dimension( * )	dlf,
		real, dimension( * )	df,
		real, dimension( * )	duf,
		real, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		real, dimension( ldb, * )	b,
		integer	ldb,
		real, dimension( ldx, * )	x,
		integer	ldx,
		real, dimension( * )	ferr,
		real, dimension( * )	berr,
		real, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

SGTRFS

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

Purpose:

!>
!> SGTRFS improves the computed solution to a system of linear
!> equations when the coefficient matrix is tridiagonal, and provides
!> error bounds and backward error estimates for the solution.
!> 

Parameters

[in]	TRANS	!> TRANS is CHARACTER*1 !> Specifies the form of the system of equations: !> = 'N': A * X = B (No transpose) !> = 'T': A**T * X = B (Transpose) !> = 'C': A**H * X = B (Conjugate transpose = Transpose) !>
[in]	N	!> N is INTEGER !> The order of the matrix A. N >= 0. !>
[in]	NRHS	!> NRHS is INTEGER !> The number of right hand sides, i.e., the number of columns !> of the matrix B. NRHS >= 0. !>
[in]	DL	!> DL is REAL array, dimension (N-1) !> The (n-1) subdiagonal elements of A. !>
[in]	D	!> D is REAL array, dimension (N) !> The diagonal elements of A. !>
[in]	DU	!> DU is REAL array, dimension (N-1) !> The (n-1) superdiagonal elements of A. !>
[in]	DLF	!> DLF is REAL array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A as computed by SGTTRF. !>
[in]	DF	!> DF is REAL array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DUF	!> DUF is REAL array, dimension (N-1) !> The (n-1) elements of the first superdiagonal of U. !>
[in]	DU2	!> DU2 is REAL array, dimension (N-2) !> The (n-2) elements of the second superdiagonal of U. !>
[in]	IPIV	!> IPIV is INTEGER array, dimension (N) !> The pivot indices; for 1 <= i <= n, row i of the matrix was !> interchanged with row IPIV(i). IPIV(i) will always be either !> i or i+1; IPIV(i) = i indicates a row interchange was not !> required. !>
[in]	B	!> B is REAL array, dimension (LDB,NRHS) !> The right hand side matrix B. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,N). !>
[in,out]	X	!> X is REAL array, dimension (LDX,NRHS) !> On entry, the solution matrix X, as computed by SGTTRS. !> On exit, the improved solution matrix X. !>
[in]	LDX	!> LDX is INTEGER !> The leading dimension of the array X. LDX >= max(1,N). !>
[out]	FERR	!> FERR is REAL array, dimension (NRHS) !> The estimated forward error bound for each solution vector !> X(j) (the j-th column of the solution matrix X). !> If XTRUE is the true solution corresponding to X(j), FERR(j) !> is an estimated upper bound for the magnitude of the largest !> element in (X(j) - XTRUE) divided by the magnitude of the !> largest element in X(j). The estimate is as reliable as !> the estimate for RCOND, and is almost always a slight !> overestimate of the true error. !>
[out]	BERR	!> BERR is REAL array, dimension (NRHS) !> The componentwise relative backward error of each solution !> vector X(j) (i.e., the smallest relative change in !> any element of A or B that makes X(j) an exact solution). !>
[out]	WORK	!> WORK is REAL array, dimension (3*N) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (N) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !>

Internal Parameters:

!>  ITMAX is the maximum number of steps of iterative refinement.
!> 

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 206 of file sgtrfs.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          TRANS
      INTEGER            INFO, LDB, LDX, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IWORK( * )
      REAL               B( LDB, * ), BERR( * ), D( * ), DF( * ),
     $                   DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
     $                   FERR( * ), WORK( * ), X( LDX, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            ITMAX
      parameter( itmax = 5 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      REAL               TWO
      parameter( two = 2.0e+0 )
      REAL               THREE
      parameter( three = 3.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOTRAN
      CHARACTER          TRANSN, TRANST
      INTEGER            COUNT, I, J, KASE, NZ
      REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           saxpy, scopy, sgttrs, slacn2, slagtm, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH
      EXTERNAL           lsame, slamch
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      notran = lsame( trans, 'N' )
      IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
     $    lsame( trans, 'C' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( nrhs.LT.0 ) THEN
         info = -3
      ELSE IF( ldb.LT.max( 1, n ) ) THEN
         info = -13
      ELSE IF( ldx.LT.max( 1, n ) ) THEN
         info = -15
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SGTRFS', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
         DO 10 j = 1, nrhs
            ferr( j ) = zero
            berr( j ) = zero
   10    CONTINUE
         RETURN
      END IF
*
      IF( notran ) THEN
         transn = 'N'
         transt = 'T'
      ELSE
         transn = 'T'
         transt = 'N'
      END IF
*
*     NZ = maximum number of nonzero elements in each row of A, plus 1
*
      nz = 4
      eps = slamch( 'Epsilon' )
      safmin = slamch( 'Safe minimum' )
      safe1 = nz*safmin
      safe2 = safe1 / eps
*
*     Do for each right hand side
*
      DO 110 j = 1, nrhs
*
         count = 1
         lstres = three
   20    CONTINUE
*
*        Loop until stopping criterion is satisfied.
*
*        Compute residual R = B - op(A) * X,
*        where op(A) = A, A**T, or A**H, depending on TRANS.
*
         CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
         CALL slagtm( trans, n, 1, -one, dl, d, du, x( 1, j ), ldx, one,
     $                work( n+1 ), n )
*
*        Compute abs(op(A))*abs(x) + abs(b) for use in the backward
*        error bound.
*
         IF( notran ) THEN
            IF( n.EQ.1 ) THEN
               work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) )
            ELSE
               work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) ) +
     $                     abs( du( 1 )*x( 2, j ) )
               DO 30 i = 2, n - 1
                  work( i ) = abs( b( i, j ) ) +
     $                        abs( dl( i-1 )*x( i-1, j ) ) +
     $                        abs( d( i )*x( i, j ) ) +
     $                        abs( du( i )*x( i+1, j ) )
   30          CONTINUE
               work( n ) = abs( b( n, j ) ) +
     $                     abs( dl( n-1 )*x( n-1, j ) ) +
     $                     abs( d( n )*x( n, j ) )
            END IF
         ELSE
            IF( n.EQ.1 ) THEN
               work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) )
            ELSE
               work( 1 ) = abs( b( 1, j ) ) + abs( d( 1 )*x( 1, j ) ) +
     $                     abs( dl( 1 )*x( 2, j ) )
               DO 40 i = 2, n - 1
                  work( i ) = abs( b( i, j ) ) +
     $                        abs( du( i-1 )*x( i-1, j ) ) +
     $                        abs( d( i )*x( i, j ) ) +
     $                        abs( dl( i )*x( i+1, j ) )
   40          CONTINUE
               work( n ) = abs( b( n, j ) ) +
     $                     abs( du( n-1 )*x( n-1, j ) ) +
     $                     abs( d( n )*x( n, j ) )
            END IF
         END IF
*
*        Compute componentwise relative backward error from formula
*
*        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
*
*        where abs(Z) is the componentwise absolute value of the matrix
*        or vector Z.  If the i-th component of the denominator is less
*        than SAFE2, then SAFE1 is added to the i-th components of the
*        numerator and denominator before dividing.
*
         s = zero
         DO 50 i = 1, n
            IF( work( i ).GT.safe2 ) THEN
               s = max( s, abs( work( n+i ) ) / work( i ) )
            ELSE
               s = max( s, ( abs( work( n+i ) )+safe1 ) /
     $             ( work( i )+safe1 ) )
            END IF
   50    CONTINUE
         berr( j ) = s
*
*        Test stopping criterion. Continue iterating if
*           1) The residual BERR(J) is larger than machine epsilon, and
*           2) BERR(J) decreased by at least a factor of 2 during the
*              last iteration, and
*           3) At most ITMAX iterations tried.
*
         IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
     $       count.LE.itmax ) THEN
*
*           Update solution and try again.
*
            CALL sgttrs( trans, n, 1, dlf, df, duf, du2, ipiv,
     $                   work( n+1 ), n, info )
            CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
            lstres = berr( j )
            count = count + 1
            GO TO 20
         END IF
*
*        Bound error from formula
*
*        norm(X - XTRUE) / norm(X) .le. FERR =
*        norm( abs(inv(op(A)))*
*           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
*
*        where
*          norm(Z) is the magnitude of the largest component of Z
*          inv(op(A)) is the inverse of op(A)
*          abs(Z) is the componentwise absolute value of the matrix or
*             vector Z
*          NZ is the maximum number of nonzeros in any row of A, plus 1
*          EPS is machine epsilon
*
*        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
*        is incremented by SAFE1 if the i-th component of
*        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
*
*        Use SLACN2 to estimate the infinity-norm of the matrix
*           inv(op(A)) * diag(W),
*        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
*
         DO 60 i = 1, n
            IF( work( i ).GT.safe2 ) THEN
               work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
            ELSE
               work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
            END IF
   60    CONTINUE
*
         kase = 0
   70    CONTINUE
         CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
     $                kase, isave )
         IF( kase.NE.0 ) THEN
            IF( kase.EQ.1 ) THEN
*
*              Multiply by diag(W)*inv(op(A)**T).
*
               CALL sgttrs( transt, n, 1, dlf, df, duf, du2, ipiv,
     $                      work( n+1 ), n, info )
               DO 80 i = 1, n
                  work( n+i ) = work( i )*work( n+i )
   80          CONTINUE
            ELSE
*
*              Multiply by inv(op(A))*diag(W).
*
               DO 90 i = 1, n
                  work( n+i ) = work( i )*work( n+i )
   90          CONTINUE
               CALL sgttrs( transn, n, 1, dlf, df, duf, du2, ipiv,
     $                      work( n+1 ), n, info )
            END IF
            GO TO 70
         END IF
*
*        Normalize error.
*
         lstres = zero
         DO 100 i = 1, n
            lstres = max( lstres, abs( x( i, j ) ) )
  100    CONTINUE
         IF( lstres.NE.zero )
     $      ferr( j ) = ferr( j ) / lstres
*
  110 CONTINUE
*
      RETURN
*
*     End of SGTRFS
*

sgttrf()

subroutine sgttrf	(	integer	n,
		real, dimension( * )	dl,
		real, dimension( * )	d,
		real, dimension( * )	du,
		real, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		integer	info )

SGTTRF

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

Purpose:

!>
!> SGTTRF computes an LU factorization of a real tridiagonal matrix A
!> using elimination with partial pivoting and row interchanges.
!>
!> The factorization has the form
!>    A = L * U
!> where L is a product of permutation and unit lower bidiagonal
!> matrices and U is upper triangular with nonzeros in only the main
!> diagonal and first two superdiagonals.
!> 

Parameters

[in]	N	!> N is INTEGER !> The order of the matrix A. !>
[in,out]	DL	!> DL is REAL array, dimension (N-1) !> On entry, DL must contain the (n-1) sub-diagonal elements of !> A. !> !> On exit, DL is overwritten by the (n-1) multipliers that !> define the matrix L from the LU factorization of A. !>
[in,out]	D	!> D is REAL array, dimension (N) !> On entry, D must contain the diagonal elements of A. !> !> On exit, D is overwritten by the n diagonal elements of the !> upper triangular matrix U from the LU factorization of A. !>
[in,out]	DU	!> DU is REAL array, dimension (N-1) !> On entry, DU must contain the (n-1) super-diagonal elements !> of A. !> !> On exit, DU is overwritten by the (n-1) elements of the first !> super-diagonal of U. !>
[out]	DU2	!> DU2 is REAL array, dimension (N-2) !> On exit, DU2 is overwritten by the (n-2) elements of the !> second super-diagonal of U. !>
[out]	IPIV	!> IPIV is INTEGER array, dimension (N) !> The pivot indices; for 1 <= i <= n, row i of the matrix was !> interchanged with row IPIV(i). IPIV(i) will always be either !> i or i+1; IPIV(i) = i indicates a row interchange was not !> required. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -k, the k-th argument had an illegal value !> > 0: if INFO = k, U(k,k) is exactly zero. The factorization !> has been completed, but the factor U 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.

Definition at line 123 of file sgttrf.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 ..
      INTEGER            INFO, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      REAL               FACT, TEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla
*     ..
*     .. Executable Statements ..
*
      info = 0
      IF( n.LT.0 ) THEN
         info = -1
         CALL xerbla( 'SGTTRF', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Initialize IPIV(i) = i and DU2(I) = 0
*
      DO 10 i = 1, n
         ipiv( i ) = i
   10 CONTINUE
      DO 20 i = 1, n - 2
         du2( i ) = zero
   20 CONTINUE
*
      DO 30 i = 1, n - 2
         IF( abs( d( i ) ).GE.abs( dl( i ) ) ) THEN
*
*           No row interchange required, eliminate DL(I)
*
            IF( d( i ).NE.zero ) THEN
               fact = dl( i ) / d( i )
               dl( i ) = fact
               d( i+1 ) = d( i+1 ) - fact*du( i )
            END IF
         ELSE
*
*           Interchange rows I and I+1, eliminate DL(I)
*
            fact = d( i ) / dl( i )
            d( i ) = dl( i )
            dl( i ) = fact
            temp = du( i )
            du( i ) = d( i+1 )
            d( i+1 ) = temp - fact*d( i+1 )
            du2( i ) = du( i+1 )
            du( i+1 ) = -fact*du( i+1 )
            ipiv( i ) = i + 1
         END IF
   30 CONTINUE
      IF( n.GT.1 ) THEN
         i = n - 1
         IF( abs( d( i ) ).GE.abs( dl( i ) ) ) THEN
            IF( d( i ).NE.zero ) THEN
               fact = dl( i ) / d( i )
               dl( i ) = fact
               d( i+1 ) = d( i+1 ) - fact*du( i )
            END IF
         ELSE
            fact = d( i ) / dl( i )
            d( i ) = dl( i )
            dl( i ) = fact
            temp = du( i )
            du( i ) = d( i+1 )
            d( i+1 ) = temp - fact*d( i+1 )
            ipiv( i ) = i + 1
         END IF
      END IF
*
*     Check for a zero on the diagonal of U.
*
      DO 40 i = 1, n
         IF( d( i ).EQ.zero ) THEN
            info = i
            GO TO 50
         END IF
   40 CONTINUE
   50 CONTINUE
*
      RETURN
*
*     End of SGTTRF
*

sgttrs()

subroutine sgttrs	(	character	trans,
		integer	n,
		integer	nrhs,
		real, dimension( * )	dl,
		real, dimension( * )	d,
		real, dimension( * )	du,
		real, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		real, dimension( ldb, * )	b,
		integer	ldb,
		integer	info )

SGTTRS

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

Purpose:

!>
!> SGTTRS solves one of the systems of equations
!>    A*X = B  or  A**T*X = B,
!> with a tridiagonal matrix A using the LU factorization computed
!> by SGTTRF.
!> 

Parameters

[in]	TRANS	!> TRANS is CHARACTER*1 !> Specifies the form of the system of equations. !> = 'N': A * X = B (No transpose) !> = 'T': A**T* X = B (Transpose) !> = 'C': A**T* X = B (Conjugate transpose = Transpose) !>
[in]	N	!> N is INTEGER !> The order of the matrix A. !>
[in]	NRHS	!> NRHS is INTEGER !> The number of right hand sides, i.e., the number of columns !> of the matrix B. NRHS >= 0. !>
[in]	DL	!> DL is REAL array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A. !>
[in]	D	!> D is REAL array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is REAL array, dimension (N-1) !> The (n-1) elements of the first super-diagonal of U. !>
[in]	DU2	!> DU2 is REAL array, dimension (N-2) !> The (n-2) elements of the second super-diagonal of U. !>
[in]	IPIV	!> IPIV is INTEGER array, dimension (N) !> The pivot indices; for 1 <= i <= n, row i of the matrix was !> interchanged with row IPIV(i). IPIV(i) will always be either !> i or i+1; IPIV(i) = i indicates a row interchange was not !> required. !>
[in,out]	B	!> B is REAL array, dimension (LDB,NRHS) !> On entry, the matrix of right hand side vectors B. !> On exit, B is overwritten by the solution vectors X. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,N). !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 136 of file sgttrs.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          TRANS
      INTEGER            INFO, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            NOTRAN
      INTEGER            ITRANS, J, JB, NB
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      EXTERNAL           ilaenv
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgtts2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min
*     ..
*     .. Executable Statements ..
*
      info = 0
      notran = ( trans.EQ.'N' .OR. trans.EQ.'n' )
      IF( .NOT.notran .AND. .NOT.( trans.EQ.'T' .OR. trans.EQ.
     $    't' ) .AND. .NOT.( trans.EQ.'C' .OR. trans.EQ.'c' ) ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( nrhs.LT.0 ) THEN
         info = -3
      ELSE IF( ldb.LT.max( n, 1 ) ) THEN
         info = -10
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SGTTRS', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 .OR. nrhs.EQ.0 )
     $   RETURN
*
*     Decode TRANS
*
      IF( notran ) THEN
         itrans = 0
      ELSE
         itrans = 1
      END IF
*
*     Determine the number of right-hand sides to solve at a time.
*
      IF( nrhs.EQ.1 ) THEN
         nb = 1
      ELSE
         nb = max( 1, ilaenv( 1, 'SGTTRS', trans, n, nrhs, -1, -1 ) )
      END IF
*
      IF( nb.GE.nrhs ) THEN
         CALL sgtts2( itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb )
      ELSE
         DO 10 j = 1, nrhs, nb
            jb = min( nrhs-j+1, nb )
            CALL sgtts2( itrans, n, jb, dl, d, du, du2, ipiv, b( 1, j ),
     $                   ldb )
   10    CONTINUE
      END IF
*
*     End of SGTTRS
*

sgtts2()

subroutine sgtts2	(	integer	itrans,
		integer	n,
		integer	nrhs,
		real, dimension( * )	dl,
		real, dimension( * )	d,
		real, dimension( * )	du,
		real, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		real, dimension( ldb, * )	b,
		integer	ldb )

SGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.

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

Purpose:

!>
!> SGTTS2 solves one of the systems of equations
!>    A*X = B  or  A**T*X = B,
!> with a tridiagonal matrix A using the LU factorization computed
!> by SGTTRF.
!> 

Parameters

[in]	ITRANS	!> ITRANS is INTEGER !> Specifies the form of the system of equations. !> = 0: A * X = B (No transpose) !> = 1: A**T* X = B (Transpose) !> = 2: A**T* X = B (Conjugate transpose = Transpose) !>
[in]	N	!> N is INTEGER !> The order of the matrix A. !>
[in]	NRHS	!> NRHS is INTEGER !> The number of right hand sides, i.e., the number of columns !> of the matrix B. NRHS >= 0. !>
[in]	DL	!> DL is REAL array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A. !>
[in]	D	!> D is REAL array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is REAL array, dimension (N-1) !> The (n-1) elements of the first super-diagonal of U. !>
[in]	DU2	!> DU2 is REAL array, dimension (N-2) !> The (n-2) elements of the second super-diagonal of U. !>
[in]	IPIV	!> IPIV is INTEGER array, dimension (N) !> The pivot indices; for 1 <= i <= n, row i of the matrix was !> interchanged with row IPIV(i). IPIV(i) will always be either !> i or i+1; IPIV(i) = i indicates a row interchange was not !> required. !>
[in,out]	B	!> B is REAL array, dimension (LDB,NRHS) !> On entry, the matrix of right hand side vectors B. !> On exit, B is overwritten by the solution vectors X. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,N). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 127 of file sgtts2.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 ..
      INTEGER            ITRANS, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      REAL               B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, IP, J
      REAL               TEMP
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( n.EQ.0 .OR. nrhs.EQ.0 )
     $   RETURN
*
      IF( itrans.EQ.0 ) THEN
*
*        Solve A*X = B using the LU factorization of A,
*        overwriting each right hand side vector with its solution.
*
         IF( nrhs.LE.1 ) THEN
            j = 1
   10       CONTINUE
*
*           Solve L*x = b.
*
            DO 20 i = 1, n - 1
               ip = ipiv( i )
               temp = b( i+1-ip+i, j ) - dl( i )*b( ip, j )
               b( i, j ) = b( ip, j )
               b( i+1, j ) = temp
   20       CONTINUE
*
*           Solve U*x = b.
*
            b( n, j ) = b( n, j ) / d( n )
            IF( n.GT.1 )
     $         b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /
     $                       d( n-1 )
            DO 30 i = n - 2, 1, -1
               b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-du2( i )*
     $                     b( i+2, j ) ) / d( i )
   30       CONTINUE
            IF( j.LT.nrhs ) THEN
               j = j + 1
               GO TO 10
            END IF
         ELSE
            DO 60 j = 1, nrhs
*
*              Solve L*x = b.
*
               DO 40 i = 1, n - 1
                  IF( ipiv( i ).EQ.i ) THEN
                     b( i+1, j ) = b( i+1, j ) - dl( i )*b( i, j )
                  ELSE
                     temp = b( i, j )
                     b( i, j ) = b( i+1, j )
                     b( i+1, j ) = temp - dl( i )*b( i, j )
                  END IF
   40          CONTINUE
*
*              Solve U*x = b.
*
               b( n, j ) = b( n, j ) / d( n )
               IF( n.GT.1 )
     $            b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /
     $                          d( n-1 )
               DO 50 i = n - 2, 1, -1
                  b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-du2( i )*
     $                        b( i+2, j ) ) / d( i )
   50          CONTINUE
   60       CONTINUE
         END IF
      ELSE
*
*        Solve A**T * X = B.
*
         IF( nrhs.LE.1 ) THEN
*
*           Solve U**T*x = b.
*
            j = 1
   70       CONTINUE
            b( 1, j ) = b( 1, j ) / d( 1 )
            IF( n.GT.1 )
     $         b( 2, j ) = ( b( 2, j )-du( 1 )*b( 1, j ) ) / d( 2 )
            DO 80 i = 3, n
               b( i, j ) = ( b( i, j )-du( i-1 )*b( i-1, j )-du2( i-2 )*
     $                     b( i-2, j ) ) / d( i )
   80       CONTINUE
*
*           Solve L**T*x = b.
*
            DO 90 i = n - 1, 1, -1
               ip = ipiv( i )
               temp = b( i, j ) - dl( i )*b( i+1, j )
               b( i, j ) = b( ip, j )
               b( ip, j ) = temp
   90       CONTINUE
            IF( j.LT.nrhs ) THEN
               j = j + 1
               GO TO 70
            END IF
*
         ELSE
            DO 120 j = 1, nrhs
*
*              Solve U**T*x = b.
*
               b( 1, j ) = b( 1, j ) / d( 1 )
               IF( n.GT.1 )
     $            b( 2, j ) = ( b( 2, j )-du( 1 )*b( 1, j ) ) / d( 2 )
               DO 100 i = 3, n
                  b( i, j ) = ( b( i, j )-du( i-1 )*b( i-1, j )-
     $                        du2( i-2 )*b( i-2, j ) ) / d( i )
  100          CONTINUE
               DO 110 i = n - 1, 1, -1
                  IF( ipiv( i ).EQ.i ) THEN
                     b( i, j ) = b( i, j ) - dl( i )*b( i+1, j )
                  ELSE
                     temp = b( i+1, j )
                     b( i+1, j ) = b( i, j ) - dl( i )*temp
                     b( i, j ) = temp
                  END IF
  110          CONTINUE
  120       CONTINUE
         END IF
      END IF
*
*     End of SGTTS2
*