1       SUBROUTINE ZPTTS2( IUPLO, N, NRHS, D, E, B, LDB )
  2 *
  3 *  -- LAPACK routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            IUPLO, LDB, N, NRHS
 10 *     ..
 11 *     .. Array Arguments ..
 12       DOUBLE PRECISION   D( * )
 13       COMPLEX*16         B( LDB, * ), E( * )
 14 *     ..
 15 *
 16 *  Purpose
 17 *  =======
 18 *
 19 *  ZPTTS2 solves a tridiagonal system of the form
 20 *     A * X = B
 21 *  using the factorization A = U**H *D*U or A = L*D*L**H computed by ZPTTRF.
 22 *  D is a diagonal matrix specified in the vector D, U (or L) is a unit
 23 *  bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
 24 *  the vector E, and X and B are N by NRHS matrices.
 25 *
 26 *  Arguments
 27 *  =========
 28 *
 29 *  IUPLO   (input) INTEGER
 30 *          Specifies the form of the factorization and whether the
 31 *          vector E is the superdiagonal of the upper bidiagonal factor
 32 *          U or the subdiagonal of the lower bidiagonal factor L.
 33 *          = 1:  A = U**H *D*U, E is the superdiagonal of U
 34 *          = 0:  A = L*D*L**H, E is the subdiagonal of L
 35 *
 36 *  N       (input) INTEGER
 37 *          The order of the tridiagonal matrix A.  N >= 0.
 38 *
 39 *  NRHS    (input) INTEGER
 40 *          The number of right hand sides, i.e., the number of columns
 41 *          of the matrix B.  NRHS >= 0.
 42 *
 43 *  D       (input) DOUBLE PRECISION array, dimension (N)
 44 *          The n diagonal elements of the diagonal matrix D from the
 45 *          factorization A = U**H *D*U or A = L*D*L**H.
 46 *
 47 *  E       (input) COMPLEX*16 array, dimension (N-1)
 48 *          If IUPLO = 1, the (n-1) superdiagonal elements of the unit
 49 *          bidiagonal factor U from the factorization A = U**H*D*U.
 50 *          If IUPLO = 0, the (n-1) subdiagonal elements of the unit
 51 *          bidiagonal factor L from the factorization A = L*D*L**H.
 52 *
 53 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 54 *          On entry, the right hand side vectors B for the system of
 55 *          linear equations.
 56 *          On exit, the solution vectors, X.
 57 *
 58 *  LDB     (input) INTEGER
 59 *          The leading dimension of the array B.  LDB >= max(1,N).
 60 *
 61 *  =====================================================================
 62 *
 63 *     .. Local Scalars ..
 64       INTEGER            I, J
 65 *     ..
 66 *     .. External Subroutines ..
 67       EXTERNAL           ZDSCAL
 68 *     ..
 69 *     .. Intrinsic Functions ..
 70       INTRINSIC          DCONJG
 71 *     ..
 72 *     .. Executable Statements ..
 73 *
 74 *     Quick return if possible
 75 *
 76       IF( N.LE.1 ) THEN
 77          IF( N.EQ.1 )
 78      $      CALL ZDSCAL( NRHS, 1.D0 / D( 1 ), B, LDB )
 79          RETURN
 80       END IF
 81 *
 82       IF( IUPLO.EQ.1 ) THEN
 83 *
 84 *        Solve A * X = B using the factorization A = U**H *D*U,
 85 *        overwriting each right hand side vector with its solution.
 86 *
 87          IF( NRHS.LE.2 ) THEN
 88             J = 1
 89    10       CONTINUE
 90 *
 91 *           Solve U**H * x = b.
 92 *
 93             DO 20 I = 2, N
 94                B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
 95    20       CONTINUE
 96 *
 97 *           Solve D * U * x = b.
 98 *
 99             DO 30 I = 1, N
100                B( I, J ) = B( I, J ) / D( I )
101    30       CONTINUE
102             DO 40 I = N - 11-1
103                B( I, J ) = B( I, J ) - B( I+1, J )*E( I )
104    40       CONTINUE
105             IF( J.LT.NRHS ) THEN
106                J = J + 1
107                GO TO 10
108             END IF
109          ELSE
110             DO 70 J = 1, NRHS
111 *
112 *              Solve U**H * x = b.
113 *
114                DO 50 I = 2, N
115                   B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
116    50          CONTINUE
117 *
118 *              Solve D * U * x = b.
119 *
120                B( N, J ) = B( N, J ) / D( N )
121                DO 60 I = N - 11-1
122                   B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I )
123    60          CONTINUE
124    70       CONTINUE
125          END IF
126       ELSE
127 *
128 *        Solve A * X = B using the factorization A = L*D*L**H,
129 *        overwriting each right hand side vector with its solution.
130 *
131          IF( NRHS.LE.2 ) THEN
132             J = 1
133    80       CONTINUE
134 *
135 *           Solve L * x = b.
136 *
137             DO 90 I = 2, N
138                B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
139    90       CONTINUE
140 *
141 *           Solve D * L**H * x = b.
142 *
143             DO 100 I = 1, N
144                B( I, J ) = B( I, J ) / D( I )
145   100       CONTINUE
146             DO 110 I = N - 11-1
147                B( I, J ) = B( I, J ) - B( I+1, J )*DCONJG( E( I ) )
148   110       CONTINUE
149             IF( J.LT.NRHS ) THEN
150                J = J + 1
151                GO TO 80
152             END IF
153          ELSE
154             DO 140 J = 1, NRHS
155 *
156 *              Solve L * x = b.
157 *
158                DO 120 I = 2, N
159                   B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
160   120          CONTINUE
161 *
162 *              Solve D * L**H * x = b.
163 *
164                B( N, J ) = B( N, J ) / D( N )
165                DO 130 I = N - 11-1
166                   B( I, J ) = B( I, J ) / D( I ) -
167      $                        B( I+1, J )*DCONJG( E( I ) )
168   130          CONTINUE
169   140       CONTINUE
170          END IF
171       END IF
172 *
173       RETURN
174 *
175 *     End of ZPTTS2
176 *
177       END