LAPACK 3.3.1
Linear Algebra PACKage

cgtsv.f

Go to the documentation of this file.
00001       SUBROUTINE CGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, LDB, N, NRHS
00010 *     ..
00011 *     .. Array Arguments ..
00012       COMPLEX            B( LDB, * ), D( * ), DL( * ), DU( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  CGTSV  solves the equation
00019 *
00020 *     A*X = B,
00021 *
00022 *  where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
00023 *  partial pivoting.
00024 *
00025 *  Note that the equation  A**H *X = B  may be solved by interchanging the
00026 *  order of the arguments DU and DL.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.  N >= 0.
00033 *
00034 *  NRHS    (input) INTEGER
00035 *          The number of right hand sides, i.e., the number of columns
00036 *          of the matrix B.  NRHS >= 0.
00037 *
00038 *  DL      (input/output) COMPLEX array, dimension (N-1)
00039 *          On entry, DL must contain the (n-1) subdiagonal elements of
00040 *          A.
00041 *          On exit, DL is overwritten by the (n-2) elements of the
00042 *          second superdiagonal of the upper triangular matrix U from
00043 *          the LU factorization of A, in DL(1), ..., DL(n-2).
00044 *
00045 *  D       (input/output) COMPLEX array, dimension (N)
00046 *          On entry, D must contain the diagonal elements of A.
00047 *          On exit, D is overwritten by the n diagonal elements of U.
00048 *
00049 *  DU      (input/output) COMPLEX array, dimension (N-1)
00050 *          On entry, DU must contain the (n-1) superdiagonal elements
00051 *          of A.
00052 *          On exit, DU is overwritten by the (n-1) elements of the first
00053 *          superdiagonal of U.
00054 *
00055 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00056 *          On entry, the N-by-NRHS right hand side matrix B.
00057 *          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
00058 *
00059 *  LDB     (input) INTEGER
00060 *          The leading dimension of the array B.  LDB >= max(1,N).
00061 *
00062 *  INFO    (output) INTEGER
00063 *          = 0:  successful exit
00064 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00065 *          > 0:  if INFO = i, U(i,i) is exactly zero, and the solution
00066 *                has not been computed.  The factorization has not been
00067 *                completed unless i = N.
00068 *
00069 *  =====================================================================
00070 *
00071 *     .. Parameters ..
00072       COMPLEX            ZERO
00073       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
00074 *     ..
00075 *     .. Local Scalars ..
00076       INTEGER            J, K
00077       COMPLEX            MULT, TEMP, ZDUM
00078 *     ..
00079 *     .. Intrinsic Functions ..
00080       INTRINSIC          ABS, AIMAG, MAX, REAL
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           XERBLA
00084 *     ..
00085 *     .. Statement Functions ..
00086       REAL               CABS1
00087 *     ..
00088 *     .. Statement Function definitions ..
00089       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00090 *     ..
00091 *     .. Executable Statements ..
00092 *
00093       INFO = 0
00094       IF( N.LT.0 ) THEN
00095          INFO = -1
00096       ELSE IF( NRHS.LT.0 ) THEN
00097          INFO = -2
00098       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00099          INFO = -7
00100       END IF
00101       IF( INFO.NE.0 ) THEN
00102          CALL XERBLA( 'CGTSV ', -INFO )
00103          RETURN
00104       END IF
00105 *
00106       IF( N.EQ.0 )
00107      $   RETURN
00108 *
00109       DO 30 K = 1, N - 1
00110          IF( DL( K ).EQ.ZERO ) THEN
00111 *
00112 *           Subdiagonal is zero, no elimination is required.
00113 *
00114             IF( D( K ).EQ.ZERO ) THEN
00115 *
00116 *              Diagonal is zero: set INFO = K and return; a unique
00117 *              solution can not be found.
00118 *
00119                INFO = K
00120                RETURN
00121             END IF
00122          ELSE IF( CABS1( D( K ) ).GE.CABS1( DL( K ) ) ) THEN
00123 *
00124 *           No row interchange required
00125 *
00126             MULT = DL( K ) / D( K )
00127             D( K+1 ) = D( K+1 ) - MULT*DU( K )
00128             DO 10 J = 1, NRHS
00129                B( K+1, J ) = B( K+1, J ) - MULT*B( K, J )
00130    10       CONTINUE
00131             IF( K.LT.( N-1 ) )
00132      $         DL( K ) = ZERO
00133          ELSE
00134 *
00135 *           Interchange rows K and K+1
00136 *
00137             MULT = D( K ) / DL( K )
00138             D( K ) = DL( K )
00139             TEMP = D( K+1 )
00140             D( K+1 ) = DU( K ) - MULT*TEMP
00141             IF( K.LT.( N-1 ) ) THEN
00142                DL( K ) = DU( K+1 )
00143                DU( K+1 ) = -MULT*DL( K )
00144             END IF
00145             DU( K ) = TEMP
00146             DO 20 J = 1, NRHS
00147                TEMP = B( K, J )
00148                B( K, J ) = B( K+1, J )
00149                B( K+1, J ) = TEMP - MULT*B( K+1, J )
00150    20       CONTINUE
00151          END IF
00152    30 CONTINUE
00153       IF( D( N ).EQ.ZERO ) THEN
00154          INFO = N
00155          RETURN
00156       END IF
00157 *
00158 *     Back solve with the matrix U from the factorization.
00159 *
00160       DO 50 J = 1, NRHS
00161          B( N, J ) = B( N, J ) / D( N )
00162          IF( N.GT.1 )
00163      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
00164          DO 40 K = N - 2, 1, -1
00165             B( K, J ) = ( B( K, J )-DU( K )*B( K+1, J )-DL( K )*
00166      $                  B( K+2, J ) ) / D( K )
00167    40    CONTINUE
00168    50 CONTINUE
00169 *
00170       RETURN
00171 *
00172 *     End of CGTSV
00173 *
00174       END
 All Files Functions