LAPACK: SRC/dgtsv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, LDB, N, NRHS
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DGTSV  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'*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) DOUBLE PRECISION array, dimension (N-1)
00039 *          On entry, DL must contain the (n-1) sub-diagonal elements of
00040 *          A.
00041 *
00042 *          On exit, DL is overwritten by the (n-2) elements of the
00043 *          second super-diagonal of the upper triangular matrix U from
00044 *          the LU factorization of A, in DL(1), ..., DL(n-2).
00045 *
00046 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00047 *          On entry, D must contain the diagonal elements of A.
00048 *
00049 *          On exit, D is overwritten by the n diagonal elements of U.
00050 *
00051 *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)
00052 *          On entry, DU must contain the (n-1) super-diagonal elements
00053 *          of A.
00054 *
00055 *          On exit, DU is overwritten by the (n-1) elements of the first
00056 *          super-diagonal of U.
00057 *
00058 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00059 *          On entry, the N by NRHS matrix of right hand side matrix B.
00060 *          On exit, if INFO = 0, the N by NRHS solution matrix X.
00061 *
00062 *  LDB     (input) INTEGER
00063 *          The leading dimension of the array B.  LDB >= max(1,N).
00064 *
00065 *  INFO    (output) INTEGER
00066 *          = 0: successful exit
00067 *          < 0: if INFO = -i, the i-th argument had an illegal value
00068 *          > 0: if INFO = i, U(i,i) is exactly zero, and the solution
00069 *               has not been computed.  The factorization has not been
00070 *               completed unless i = N.
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       DOUBLE PRECISION   ZERO
00076       PARAMETER          ( ZERO = 0.0D+0 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       INTEGER            I, J
00080       DOUBLE PRECISION   FACT, TEMP
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          ABS, MAX
00084 *     ..
00085 *     .. External Subroutines ..
00086       EXTERNAL           XERBLA
00087 *     ..
00088 *     .. Executable Statements ..
00089 *
00090       INFO = 0
00091       IF( N.LT.0 ) THEN
00092          INFO = -1
00093       ELSE IF( NRHS.LT.0 ) THEN
00094          INFO = -2
00095       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00096          INFO = -7
00097       END IF
00098       IF( INFO.NE.0 ) THEN
00099          CALL XERBLA( 'DGTSV ', -INFO )
00100          RETURN
00101       END IF
00102 *
00103       IF( N.EQ.0 )
00104      $   RETURN
00105 *
00106       IF( NRHS.EQ.1 ) THEN
00107          DO 10 I = 1, N - 2
00108             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00109 *
00110 *              No row interchange required
00111 *
00112                IF( D( I ).NE.ZERO ) THEN
00113                   FACT = DL( I ) / D( I )
00114                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
00115                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
00116                ELSE
00117                   INFO = I
00118                   RETURN
00119                END IF
00120                DL( I ) = ZERO
00121             ELSE
00122 *
00123 *              Interchange rows I and I+1
00124 *
00125                FACT = D( I ) / DL( I )
00126                D( I ) = DL( I )
00127                TEMP = D( I+1 )
00128                D( I+1 ) = DU( I ) - FACT*TEMP
00129                DL( I ) = DU( I+1 )
00130                DU( I+1 ) = -FACT*DL( I )
00131                DU( I ) = TEMP
00132                TEMP = B( I, 1 )
00133                B( I, 1 ) = B( I+1, 1 )
00134                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
00135             END IF
00136    10    CONTINUE
00137          IF( N.GT.1 ) THEN
00138             I = N - 1
00139             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00140                IF( D( I ).NE.ZERO ) THEN
00141                   FACT = DL( I ) / D( I )
00142                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
00143                   B( I+1, 1 ) = B( I+1, 1 ) - FACT*B( I, 1 )
00144                ELSE
00145                   INFO = I
00146                   RETURN
00147                END IF
00148             ELSE
00149                FACT = D( I ) / DL( I )
00150                D( I ) = DL( I )
00151                TEMP = D( I+1 )
00152                D( I+1 ) = DU( I ) - FACT*TEMP
00153                DU( I ) = TEMP
00154                TEMP = B( I, 1 )
00155                B( I, 1 ) = B( I+1, 1 )
00156                B( I+1, 1 ) = TEMP - FACT*B( I+1, 1 )
00157             END IF
00158          END IF
00159          IF( D( N ).EQ.ZERO ) THEN
00160             INFO = N
00161             RETURN
00162          END IF
00163       ELSE
00164          DO 40 I = 1, N - 2
00165             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00166 *
00167 *              No row interchange required
00168 *
00169                IF( D( I ).NE.ZERO ) THEN
00170                   FACT = DL( I ) / D( I )
00171                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
00172                   DO 20 J = 1, NRHS
00173                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
00174    20             CONTINUE
00175                ELSE
00176                   INFO = I
00177                   RETURN
00178                END IF
00179                DL( I ) = ZERO
00180             ELSE
00181 *
00182 *              Interchange rows I and I+1
00183 *
00184                FACT = D( I ) / DL( I )
00185                D( I ) = DL( I )
00186                TEMP = D( I+1 )
00187                D( I+1 ) = DU( I ) - FACT*TEMP
00188                DL( I ) = DU( I+1 )
00189                DU( I+1 ) = -FACT*DL( I )
00190                DU( I ) = TEMP
00191                DO 30 J = 1, NRHS
00192                   TEMP = B( I, J )
00193                   B( I, J ) = B( I+1, J )
00194                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
00195    30          CONTINUE
00196             END IF
00197    40    CONTINUE
00198          IF( N.GT.1 ) THEN
00199             I = N - 1
00200             IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN
00201                IF( D( I ).NE.ZERO ) THEN
00202                   FACT = DL( I ) / D( I )
00203                   D( I+1 ) = D( I+1 ) - FACT*DU( I )
00204                   DO 50 J = 1, NRHS
00205                      B( I+1, J ) = B( I+1, J ) - FACT*B( I, J )
00206    50             CONTINUE
00207                ELSE
00208                   INFO = I
00209                   RETURN
00210                END IF
00211             ELSE
00212                FACT = D( I ) / DL( I )
00213                D( I ) = DL( I )
00214                TEMP = D( I+1 )
00215                D( I+1 ) = DU( I ) - FACT*TEMP
00216                DU( I ) = TEMP
00217                DO 60 J = 1, NRHS
00218                   TEMP = B( I, J )
00219                   B( I, J ) = B( I+1, J )
00220                   B( I+1, J ) = TEMP - FACT*B( I+1, J )
00221    60          CONTINUE
00222             END IF
00223          END IF
00224          IF( D( N ).EQ.ZERO ) THEN
00225             INFO = N
00226             RETURN
00227          END IF
00228       END IF
00229 *
00230 *     Back solve with the matrix U from the factorization.
00231 *
00232       IF( NRHS.LE.2 ) THEN
00233          J = 1
00234    70    CONTINUE
00235          B( N, J ) = B( N, J ) / D( N )
00236          IF( N.GT.1 )
00237      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
00238          DO 80 I = N - 2, 1, -1
00239             B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
00240      $                  B( I+2, J ) ) / D( I )
00241    80    CONTINUE
00242          IF( J.LT.NRHS ) THEN
00243             J = J + 1
00244             GO TO 70
00245          END IF
00246       ELSE
00247          DO 100 J = 1, NRHS
00248             B( N, J ) = B( N, J ) / D( N )
00249             IF( N.GT.1 )
00250      $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
00251      $                       D( N-1 )
00252             DO 90 I = N - 2, 1, -1
00253                B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DL( I )*
00254      $                     B( I+2, J ) ) / D( I )
00255    90       CONTINUE
00256   100    CONTINUE
00257       END IF
00258 *
00259       RETURN
00260 *
00261 *     End of DGTSV
00262 *
00263       END