complex

Functions
subroutine	cgtcon (norm, n, dl, d, du, du2, ipiv, anorm, rcond, work, info)
	CGTCON
subroutine	cgtrfs (trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
	CGTRFS
subroutine	cgttrf (n, dl, d, du, du2, ipiv, info)
	CGTTRF
subroutine	cgttrs (trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
	CGTTRS
subroutine	cgtts2 (itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb)
	CGTTS2 solves a system of linear equations with a tridiagonal matrix using the LU factorization computed by sgttrf.

Detailed Description

This is the group of complex computational functions for GT matrices

Function Documentation

cgtcon()

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

CGTCON

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

Purpose:

!>
!> CGTCON estimates the reciprocal of the condition number of a complex
!> tridiagonal matrix A using the LU factorization as computed by
!> CGTTRF.
!>
!> 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 COMPLEX array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A as computed by CGTTRF. !>
[in]	D	!> D is COMPLEX array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is COMPLEX array, dimension (N-1) !> The (n-1) elements of the first superdiagonal of U. !>
[in]	DU2	!> DU2 is COMPLEX 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 COMPLEX array, dimension (2*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 139 of file cgtcon.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( * )
      COMPLEX            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           cgttrs, clacn2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cmplx
*     ..
*     .. 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( 'CGTCON', -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.cmplx( zero ) )
     $      RETURN
   10 CONTINUE
*
      ainvnm = zero
      IF( onenrm ) THEN
         kase1 = 1
      ELSE
         kase1 = 2
      END IF
      kase = 0
   20 CONTINUE
      CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
      IF( kase.NE.0 ) THEN
         IF( kase.EQ.kase1 ) THEN
*
*           Multiply by inv(U)*inv(L).
*
            CALL cgttrs( 'No transpose', n, 1, dl, d, du, du2, ipiv,
     $                   work, n, info )
         ELSE
*
*           Multiply by inv(L**H)*inv(U**H).
*
            CALL cgttrs( 'Conjugate 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 CGTCON
*

cgtrfs()

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

CGTRFS

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

Purpose:

!>
!> CGTRFS 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) !>
[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 COMPLEX array, dimension (N-1) !> The (n-1) subdiagonal elements of A. !>
[in]	D	!> D is COMPLEX array, dimension (N) !> The diagonal elements of A. !>
[in]	DU	!> DU is COMPLEX array, dimension (N-1) !> The (n-1) superdiagonal elements of A. !>
[in]	DLF	!> DLF is COMPLEX array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A as computed by CGTTRF. !>
[in]	DF	!> DF is COMPLEX array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DUF	!> DUF is COMPLEX array, dimension (N-1) !> The (n-1) elements of the first superdiagonal of U. !>
[in]	DU2	!> DU2 is COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDX,NRHS) !> On entry, the solution matrix X, as computed by CGTTRS. !> 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 COMPLEX array, dimension (2*N) !>
[out]	RWORK	!> RWORK is REAL 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 207 of file cgtrfs.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( * )
      REAL               BERR( * ), FERR( * ), RWORK( * )
      COMPLEX            B( LDB, * ), D( * ), DF( * ), DL( * ),
     $                   DLF( * ), DU( * ), DU2( * ), DUF( * ),
     $                   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
      COMPLEX            ZDUM
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           caxpy, ccopy, cgttrs, clacn2, clagtm, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, cmplx, max, real
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH
      EXTERNAL           lsame, slamch
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
*     ..
*     .. 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( 'CGTRFS', -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 = 'C'
      ELSE
         transn = 'C'
         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 ccopy( n, b( 1, j ), 1, work, 1 )
         CALL clagtm( trans, n, 1, -one, dl, d, du, x( 1, j ), ldx, one,
     $                work, n )
*
*        Compute abs(op(A))*abs(x) + abs(b) for use in the backward
*        error bound.
*
         IF( notran ) THEN
            IF( n.EQ.1 ) THEN
               rwork( 1 ) = cabs1( b( 1, j ) ) +
     $                      cabs1( d( 1 ) )*cabs1( x( 1, j ) )
            ELSE
               rwork( 1 ) = cabs1( b( 1, j ) ) +
     $                      cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
     $                      cabs1( du( 1 ) )*cabs1( x( 2, j ) )
               DO 30 i = 2, n - 1
                  rwork( i ) = cabs1( b( i, j ) ) +
     $                         cabs1( dl( i-1 ) )*cabs1( x( i-1, j ) ) +
     $                         cabs1( d( i ) )*cabs1( x( i, j ) ) +
     $                         cabs1( du( i ) )*cabs1( x( i+1, j ) )
   30          CONTINUE
               rwork( n ) = cabs1( b( n, j ) ) +
     $                      cabs1( dl( n-1 ) )*cabs1( x( n-1, j ) ) +
     $                      cabs1( d( n ) )*cabs1( x( n, j ) )
            END IF
         ELSE
            IF( n.EQ.1 ) THEN
               rwork( 1 ) = cabs1( b( 1, j ) ) +
     $                      cabs1( d( 1 ) )*cabs1( x( 1, j ) )
            ELSE
               rwork( 1 ) = cabs1( b( 1, j ) ) +
     $                      cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
     $                      cabs1( dl( 1 ) )*cabs1( x( 2, j ) )
               DO 40 i = 2, n - 1
                  rwork( i ) = cabs1( b( i, j ) ) +
     $                         cabs1( du( i-1 ) )*cabs1( x( i-1, j ) ) +
     $                         cabs1( d( i ) )*cabs1( x( i, j ) ) +
     $                         cabs1( dl( i ) )*cabs1( x( i+1, j ) )
   40          CONTINUE
               rwork( n ) = cabs1( b( n, j ) ) +
     $                      cabs1( du( n-1 ) )*cabs1( x( n-1, j ) ) +
     $                      cabs1( d( n ) )*cabs1( 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( rwork( i ).GT.safe2 ) THEN
               s = max( s, cabs1( work( i ) ) / rwork( i ) )
            ELSE
               s = max( s, ( cabs1( work( i ) )+safe1 ) /
     $             ( rwork( 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 cgttrs( trans, n, 1, dlf, df, duf, du2, ipiv, work, n,
     $                   info )
            CALL caxpy( n, cmplx( one ), work, 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 CLACN2 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( rwork( i ).GT.safe2 ) THEN
               rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
            ELSE
               rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
     $                      safe1
            END IF
   60    CONTINUE
*
         kase = 0
   70    CONTINUE
         CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
         IF( kase.NE.0 ) THEN
            IF( kase.EQ.1 ) THEN
*
*              Multiply by diag(W)*inv(op(A)**H).
*
               CALL cgttrs( transt, n, 1, dlf, df, duf, du2, ipiv, work,
     $                      n, info )
               DO 80 i = 1, n
                  work( i ) = rwork( i )*work( i )
   80          CONTINUE
            ELSE
*
*              Multiply by inv(op(A))*diag(W).
*
               DO 90 i = 1, n
                  work( i ) = rwork( i )*work( i )
   90          CONTINUE
               CALL cgttrs( transn, n, 1, dlf, df, duf, du2, ipiv, work,
     $                      n, info )
            END IF
            GO TO 70
         END IF
*
*        Normalize error.
*
         lstres = zero
         DO 100 i = 1, n
            lstres = max( lstres, cabs1( x( i, j ) ) )
  100    CONTINUE
         IF( lstres.NE.zero )
     $      ferr( j ) = ferr( j ) / lstres
*
  110 CONTINUE
*
      RETURN
*
*     End of CGTRFS
*

cgttrf()

subroutine cgttrf	(	integer	n,
		complex, dimension( * )	dl,
		complex, dimension( * )	d,
		complex, dimension( * )	du,
		complex, dimension( * )	du2,
		integer, dimension( * )	ipiv,
		integer	info )

CGTTRF

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

Purpose:

!>
!> CGTTRF computes an LU factorization of a complex 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 cgttrf.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( * )
      COMPLEX            D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      COMPLEX            FACT, TEMP, ZDUM
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, real
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
*     ..
*     .. Executable Statements ..
*
      info = 0
      IF( n.LT.0 ) THEN
         info = -1
         CALL xerbla( 'CGTTRF', -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( cabs1( d( i ) ).GE.cabs1( dl( i ) ) ) THEN
*
*           No row interchange required, eliminate DL(I)
*
            IF( cabs1( 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( cabs1( d( i ) ).GE.cabs1( dl( i ) ) ) THEN
            IF( cabs1( 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( cabs1( d( i ) ).EQ.zero ) THEN
            info = i
            GO TO 50
         END IF
   40 CONTINUE
   50 CONTINUE
*
      RETURN
*
*     End of CGTTRF
*

cgttrs()

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

CGTTRS

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

Purpose:

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

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) !>
[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 COMPLEX array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A. !>
[in]	D	!> D is COMPLEX array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is COMPLEX array, dimension (N-1) !> The (n-1) elements of the first super-diagonal of U. !>
[in]	DU2	!> DU2 is COMPLEX 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 COMPLEX 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 = -k, the k-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 cgttrs.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( * )
      COMPLEX            B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            NOTRAN
      INTEGER            ITRANS, J, JB, NB
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      EXTERNAL           ilaenv
*     ..
*     .. External Subroutines ..
      EXTERNAL           cgtts2, 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( 'CGTTRS', -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 IF( trans.EQ.'T' .OR. trans.EQ.'t' ) THEN
         itrans = 1
      ELSE
         itrans = 2
      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, 'CGTTRS', trans, n, nrhs, -1, -1 ) )
      END IF
*
      IF( nb.GE.nrhs ) THEN
         CALL cgtts2( itrans, n, nrhs, dl, d, du, du2, ipiv, b, ldb )
      ELSE
         DO 10 j = 1, nrhs, nb
            jb = min( nrhs-j+1, nb )
            CALL cgtts2( itrans, n, jb, dl, d, du, du2, ipiv, b( 1, j ),
     $                   ldb )
   10    CONTINUE
      END IF
*
*     End of CGTTRS
*

cgtts2()

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

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

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

Purpose:

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

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**H * X = B (Conjugate 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 COMPLEX array, dimension (N-1) !> The (n-1) multipliers that define the matrix L from the !> LU factorization of A. !>
[in]	D	!> D is COMPLEX array, dimension (N) !> The n diagonal elements of the upper triangular matrix U from !> the LU factorization of A. !>
[in]	DU	!> DU is COMPLEX array, dimension (N-1) !> The (n-1) elements of the first super-diagonal of U. !>
[in]	DU2	!> DU2 is COMPLEX 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 COMPLEX 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 cgtts2.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( * )
      COMPLEX            B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, J
      COMPLEX            TEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          conjg
*     ..
*     .. 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
               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
   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 IF( itrans.EQ.1 ) THEN
*
*        Solve A**T * X = B.
*
         IF( nrhs.LE.1 ) THEN
            j = 1
   70       CONTINUE
*
*           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 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
               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
   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
*
*           Solve L**T * x = b.
*
               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
      ELSE
*
*        Solve A**H * X = B.
*
         IF( nrhs.LE.1 ) THEN
            j = 1
  130       CONTINUE
*
*           Solve U**H * x = b.
*
            b( 1, j ) = b( 1, j ) / conjg( d( 1 ) )
            IF( n.GT.1 )
     $         b( 2, j ) = ( b( 2, j )-conjg( du( 1 ) )*b( 1, j ) ) /
     $                     conjg( d( 2 ) )
            DO 140 i = 3, n
               b( i, j ) = ( b( i, j )-conjg( du( i-1 ) )*b( i-1, j )-
     $                     conjg( du2( i-2 ) )*b( i-2, j ) ) /
     $                     conjg( d( i ) )
  140       CONTINUE
*
*           Solve L**H * x = b.
*
            DO 150 i = n - 1, 1, -1
               IF( ipiv( i ).EQ.i ) THEN
                  b( i, j ) = b( i, j ) - conjg( dl( i ) )*b( i+1, j )
               ELSE
                  temp = b( i+1, j )
                  b( i+1, j ) = b( i, j ) - conjg( dl( i ) )*temp
                  b( i, j ) = temp
               END IF
  150       CONTINUE
            IF( j.LT.nrhs ) THEN
               j = j + 1
               GO TO 130
            END IF
         ELSE
            DO 180 j = 1, nrhs
*
*           Solve U**H * x = b.
*
               b( 1, j ) = b( 1, j ) / conjg( d( 1 ) )
               IF( n.GT.1 )
     $            b( 2, j ) = ( b( 2, j )-conjg( du( 1 ) )*b( 1, j ) ) /
     $                        conjg( d( 2 ) )
               DO 160 i = 3, n
                  b( i, j ) = ( b( i, j )-conjg( du( i-1 ) )*
     $                        b( i-1, j )-conjg( du2( i-2 ) )*
     $                        b( i-2, j ) ) / conjg( d( i ) )
  160          CONTINUE
*
*           Solve L**H * x = b.
*
               DO 170 i = n - 1, 1, -1
                  IF( ipiv( i ).EQ.i ) THEN
                     b( i, j ) = b( i, j ) - conjg( dl( i ) )*
     $                           b( i+1, j )
                  ELSE
                     temp = b( i+1, j )
                     b( i+1, j ) = b( i, j ) - conjg( dl( i ) )*temp
                     b( i, j ) = temp
                  END IF
  170          CONTINUE
  180       CONTINUE
         END IF
      END IF
*
*     End of CGTTS2
*