SRC/zgtsv.f

   1 *> \brief <b> ZGTSV computes the solution to system of linear equations A * X = B for GT matrices </b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGTSV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgtsv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgtsv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgtsv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, LDB, N, NRHS
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       COMPLEX*16         B( LDB, * ), D( * ), DL( * ), DU( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> ZGTSV  solves the equation
  37 *>
  38 *>    A*X = B,
  39 *>
  40 *> where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
  41 *> partial pivoting.
  42 *>
  43 *> Note that the equation  A**T *X = B  may be solved by interchanging the
  44 *> order of the arguments DU and DL.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] N
  51 *> \verbatim
  52 *>          N is INTEGER
  53 *>          The order of the matrix A.  N >= 0.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] NRHS
  57 *> \verbatim
  58 *>          NRHS is INTEGER
  59 *>          The number of right hand sides, i.e., the number of columns
  60 *>          of the matrix B.  NRHS >= 0.
  61 *> \endverbatim
  62 *>
  63 *> \param[in,out] DL
  64 *> \verbatim
  65 *>          DL is COMPLEX*16 array, dimension (N-1)
  66 *>          On entry, DL must contain the (n-1) subdiagonal elements of
  67 *>          A.
  68 *>          On exit, DL is overwritten by the (n-2) elements of the
  69 *>          second superdiagonal of the upper triangular matrix U from
  70 *>          the LU factorization of A, in DL(1), ..., DL(n-2).
  71 *> \endverbatim
  72 *>
  73 *> \param[in,out] D
  74 *> \verbatim
  75 *>          D is COMPLEX*16 array, dimension (N)
  76 *>          On entry, D must contain the diagonal elements of A.
  77 *>          On exit, D is overwritten by the n diagonal elements of U.
  78 *> \endverbatim
  79 *>
  80 *> \param[in,out] DU
  81 *> \verbatim
  82 *>          DU is COMPLEX*16 array, dimension (N-1)
  83 *>          On entry, DU must contain the (n-1) superdiagonal elements
  84 *>          of A.
  85 *>          On exit, DU is overwritten by the (n-1) elements of the first
  86 *>          superdiagonal of U.
  87 *> \endverbatim
  88 *>
  89 *> \param[in,out] B
  90 *> \verbatim
  91 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
  92 *>          On entry, the N-by-NRHS right hand side matrix B.
  93 *>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] LDB
  97 *> \verbatim
  98 *>          LDB is INTEGER
  99 *>          The leading dimension of the array B.  LDB >= max(1,N).
 100 *> \endverbatim
 101 *>
 102 *> \param[out] INFO
 103 *> \verbatim
 104 *>          INFO is INTEGER
 105 *>          = 0:  successful exit
 106 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 107 *>          > 0:  if INFO = i, U(i,i) is exactly zero, and the solution
 108 *>                has not been computed.  The factorization has not been
 109 *>                completed unless i = N.
 110 *> \endverbatim
 111 *
 112 *  Authors:
 113 *  ========
 114 *
 115 *> \author Univ. of Tennessee
 116 *> \author Univ. of California Berkeley
 117 *> \author Univ. of Colorado Denver
 118 *> \author NAG Ltd.
 119 *
 120 *> \date September 2012
 121 *
 122 *> \ingroup complex16GTsolve
 123 *
 124 *  =====================================================================
 125       SUBROUTINE ZGTSV( N, NRHS, DL, D, DU, B, LDB, INFO )
 126 *
 127 *  -- LAPACK driver routine (version 3.4.2) --
 128 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 129 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 130 *     September 2012
 131 *
 132 *     .. Scalar Arguments ..
 133       INTEGER            INFO, LDB, N, NRHS
 134 *     ..
 135 *     .. Array Arguments ..
 136       COMPLEX*16         B( LDB, * ), D( * ), DL( * ), DU( * )
 137 *     ..
 138 *
 139 *  =====================================================================
 140 *
 141 *     .. Parameters ..
 142       COMPLEX*16         ZERO
 143       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 144 *     ..
 145 *     .. Local Scalars ..
 146       INTEGER            J, K
 147       COMPLEX*16         MULT, TEMP, ZDUM
 148 *     ..
 149 *     .. Intrinsic Functions ..
 150       INTRINSIC          ABS, DBLE, DIMAG, MAX
 151 *     ..
 152 *     .. External Subroutines ..
 153       EXTERNAL           XERBLA
 154 *     ..
 155 *     .. Statement Functions ..
 156       DOUBLE PRECISION   CABS1
 157 *     ..
 158 *     .. Statement Function definitions ..
 159       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 160 *     ..
 161 *     .. Executable Statements ..
 162 *
 163       INFO = 0
 164       IF( N.LT.0 ) THEN
 165          INFO = -1
 166       ELSE IF( NRHS.LT.0 ) THEN
 167          INFO = -2
 168       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 169          INFO = -7
 170       END IF
 171       IF( INFO.NE.0 ) THEN
 172          CALL XERBLA( 'ZGTSV ', -INFO )
 173          RETURN
 174       END IF
 175 *
 176       IF( N.EQ.0 )
 177      $   RETURN
 178 *
 179       DO 30 K = 1, N - 1
 180          IF( DL( K ).EQ.ZERO ) THEN
 181 *
 182 *           Subdiagonal is zero, no elimination is required.
 183 *
 184             IF( D( K ).EQ.ZERO ) THEN
 185 *
 186 *              Diagonal is zero: set INFO = K and return; a unique
 187 *              solution can not be found.
 188 *
 189                INFO = K
 190                RETURN
 191             END IF
 192          ELSE IF( CABS1( D( K ) ).GE.CABS1( DL( K ) ) ) THEN
 193 *
 194 *           No row interchange required
 195 *
 196             MULT = DL( K ) / D( K )
 197             D( K+1 ) = D( K+1 ) - MULT*DU( K )
 198             DO 10 J = 1, NRHS
 199                B( K+1, J ) = B( K+1, J ) - MULT*B( K, J )
 200    10       CONTINUE
 201             IF( K.LT.( N-1 ) )
 202      $         DL( K ) = ZERO
 203          ELSE
 204 *
 205 *           Interchange rows K and K+1
 206 *
 207             MULT = D( K ) / DL( K )
 208             D( K ) = DL( K )
 209             TEMP = D( K+1 )
 210             D( K+1 ) = DU( K ) - MULT*TEMP
 211             IF( K.LT.( N-1 ) ) THEN
 212                DL( K ) = DU( K+1 )
 213                DU( K+1 ) = -MULT*DL( K )
 214             END IF
 215             DU( K ) = TEMP
 216             DO 20 J = 1, NRHS
 217                TEMP = B( K, J )
 218                B( K, J ) = B( K+1, J )
 219                B( K+1, J ) = TEMP - MULT*B( K+1, J )
 220    20       CONTINUE
 221          END IF
 222    30 CONTINUE
 223       IF( D( N ).EQ.ZERO ) THEN
 224          INFO = N
 225          RETURN
 226       END IF
 227 *
 228 *     Back solve with the matrix U from the factorization.
 229 *
 230       DO 50 J = 1, NRHS
 231          B( N, J ) = B( N, J ) / D( N )
 232          IF( N.GT.1 )
 233      $      B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D( N-1 )
 234          DO 40 K = N - 2, 1, -1
 235             B( K, J ) = ( B( K, J )-DU( K )*B( K+1, J )-DL( K )*
 236      $                  B( K+2, J ) ) / D( K )
 237    40    CONTINUE
 238    50 CONTINUE
 239 *
 240       RETURN
 241 *
 242 *     End of ZGTSV
 243 *
 244       END