1       SUBROUTINE DLAGTM( TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA,
  2      $                   B, LDB )
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.3.1) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *  -- April 2011                                                      --
  8 *
  9 *     .. Scalar Arguments ..
 10       CHARACTER          TRANS
 11       INTEGER            LDB, LDX, N, NRHS
 12       DOUBLE PRECISION   ALPHA, BETA
 13 *     ..
 14 *     .. Array Arguments ..
 15       DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * ),
 16      $                   X( LDX, * )
 17 *     ..
 18 *
 19 *  Purpose
 20 *  =======
 21 *
 22 *  DLAGTM performs a matrix-vector product of the form
 23 *
 24 *     B := alpha * A * X + beta * B
 25 *
 26 *  where A is a tridiagonal matrix of order N, B and X are N by NRHS
 27 *  matrices, and alpha and beta are real scalars, each of which may be
 28 *  0., 1., or -1.
 29 *
 30 *  Arguments
 31 *  =========
 32 *
 33 *  TRANS   (input) CHARACTER*1
 34 *          Specifies the operation applied to A.
 35 *          = 'N':  No transpose, B := alpha * A * X + beta * B
 36 *          = 'T':  Transpose,    B := alpha * A'* X + beta * B
 37 *          = 'C':  Conjugate transpose = Transpose
 38 *
 39 *  N       (input) INTEGER
 40 *          The order of the matrix A.  N >= 0.
 41 *
 42 *  NRHS    (input) INTEGER
 43 *          The number of right hand sides, i.e., the number of columns
 44 *          of the matrices X and B.
 45 *
 46 *  ALPHA   (input) DOUBLE PRECISION
 47 *          The scalar alpha.  ALPHA must be 0., 1., or -1.; otherwise,
 48 *          it is assumed to be 0.
 49 *
 50 *  DL      (input) DOUBLE PRECISION array, dimension (N-1)
 51 *          The (n-1) sub-diagonal elements of T.
 52 *
 53 *  D       (input) DOUBLE PRECISION array, dimension (N)
 54 *          The diagonal elements of T.
 55 *
 56 *  DU      (input) DOUBLE PRECISION array, dimension (N-1)
 57 *          The (n-1) super-diagonal elements of T.
 58 *
 59 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
 60 *          The N by NRHS matrix X.
 61 *  LDX     (input) INTEGER
 62 *          The leading dimension of the array X.  LDX >= max(N,1).
 63 *
 64 *  BETA    (input) DOUBLE PRECISION
 65 *          The scalar beta.  BETA must be 0., 1., or -1.; otherwise,
 66 *          it is assumed to be 1.
 67 *
 68 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 69 *          On entry, the N by NRHS matrix B.
 70 *          On exit, B is overwritten by the matrix expression
 71 *          B := alpha * A * X + beta * B.
 72 *
 73 *  LDB     (input) INTEGER
 74 *          The leading dimension of the array B.  LDB >= max(N,1).
 75 *
 76 *  =====================================================================
 77 *
 78 *     .. Parameters ..
 79       DOUBLE PRECISION   ONE, ZERO
 80       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 81 *     ..
 82 *     .. Local Scalars ..
 83       INTEGER            I, J
 84 *     ..
 85 *     .. External Functions ..
 86       LOGICAL            LSAME
 87       EXTERNAL           LSAME
 88 *     ..
 89 *     .. Executable Statements ..
 90 *
 91       IF( N.EQ.0 )
 92      $   RETURN
 93 *
 94 *     Multiply B by BETA if BETA.NE.1.
 95 *
 96       IF( BETA.EQ.ZERO ) THEN
 97          DO 20 J = 1, NRHS
 98             DO 10 I = 1, N
 99                B( I, J ) = ZERO
100    10       CONTINUE
101    20    CONTINUE
102       ELSE IF( BETA.EQ.-ONE ) THEN
103          DO 40 J = 1, NRHS
104             DO 30 I = 1, N
105                B( I, J ) = -B( I, J )
106    30       CONTINUE
107    40    CONTINUE
108       END IF
109 *
110       IF( ALPHA.EQ.ONE ) THEN
111          IF( LSAME( TRANS, 'N' ) ) THEN
112 *
113 *           Compute B := B + A*X
114 *
115             DO 60 J = 1, NRHS
116                IF( N.EQ.1 ) THEN
117                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
118                ELSE
119                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
120      $                        DU( 1 )*X( 2, J )
121                   B( N, J ) = B( N, J ) + DL( N-1 )*X( N-1, J ) +
122      $                        D( N )*X( N, J )
123                   DO 50 I = 2, N - 1
124                      B( I, J ) = B( I, J ) + DL( I-1 )*X( I-1, J ) +
125      $                           D( I )*X( I, J ) + DU( I )*X( I+1, J )
126    50             CONTINUE
127                END IF
128    60       CONTINUE
129          ELSE
130 *
131 *           Compute B := B + A**T*X
132 *
133             DO 80 J = 1, NRHS
134                IF( N.EQ.1 ) THEN
135                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J )
136                ELSE
137                   B( 1, J ) = B( 1, J ) + D( 1 )*X( 1, J ) +
138      $                        DL( 1 )*X( 2, J )
139                   B( N, J ) = B( N, J ) + DU( N-1 )*X( N-1, J ) +
140      $                        D( N )*X( N, J )
141                   DO 70 I = 2, N - 1
142                      B( I, J ) = B( I, J ) + DU( I-1 )*X( I-1, J ) +
143      $                           D( I )*X( I, J ) + DL( I )*X( I+1, J )
144    70             CONTINUE
145                END IF
146    80       CONTINUE
147          END IF
148       ELSE IF( ALPHA.EQ.-ONE ) THEN
149          IF( LSAME( TRANS, 'N' ) ) THEN
150 *
151 *           Compute B := B - A*X
152 *
153             DO 100 J = 1, NRHS
154                IF( N.EQ.1 ) THEN
155                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
156                ELSE
157                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
158      $                        DU( 1 )*X( 2, J )
159                   B( N, J ) = B( N, J ) - DL( N-1 )*X( N-1, J ) -
160      $                        D( N )*X( N, J )
161                   DO 90 I = 2, N - 1
162                      B( I, J ) = B( I, J ) - DL( I-1 )*X( I-1, J ) -
163      $                           D( I )*X( I, J ) - DU( I )*X( I+1, J )
164    90             CONTINUE
165                END IF
166   100       CONTINUE
167          ELSE
168 *
169 *           Compute B := B - A**T*X
170 *
171             DO 120 J = 1, NRHS
172                IF( N.EQ.1 ) THEN
173                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J )
174                ELSE
175                   B( 1, J ) = B( 1, J ) - D( 1 )*X( 1, J ) -
176      $                        DL( 1 )*X( 2, J )
177                   B( N, J ) = B( N, J ) - DU( N-1 )*X( N-1, J ) -
178      $                        D( N )*X( N, J )
179                   DO 110 I = 2, N - 1
180                      B( I, J ) = B( I, J ) - DU( I-1 )*X( I-1, J ) -
181      $                           D( I )*X( I, J ) - DL( I )*X( I+1, J )
182   110             CONTINUE
183                END IF
184   120       CONTINUE
185          END IF
186       END IF
187       RETURN
188 *
189 *     End of DLAGTM
190 *
191       END