001:       SUBROUTINE CPTTS2( IUPLO, N, NRHS, D, E, B, LDB )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            IUPLO, LDB, N, NRHS
009: *     ..
010: *     .. Array Arguments ..
011:       REAL               D( * )
012:       COMPLEX            B( LDB, * ), E( * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  CPTTS2 solves a tridiagonal system of the form
019: *     A * X = B
020: *  using the factorization A = U'*D*U or A = L*D*L' computed by CPTTRF.
021: *  D is a diagonal matrix specified in the vector D, U (or L) is a unit
022: *  bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
023: *  the vector E, and X and B are N by NRHS matrices.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  IUPLO   (input) INTEGER
029: *          Specifies the form of the factorization and whether the
030: *          vector E is the superdiagonal of the upper bidiagonal factor
031: *          U or the subdiagonal of the lower bidiagonal factor L.
032: *          = 1:  A = U'*D*U, E is the superdiagonal of U
033: *          = 0:  A = L*D*L', E is the subdiagonal of L
034: *
035: *  N       (input) INTEGER
036: *          The order of the tridiagonal matrix A.  N >= 0.
037: *
038: *  NRHS    (input) INTEGER
039: *          The number of right hand sides, i.e., the number of columns
040: *          of the matrix B.  NRHS >= 0.
041: *
042: *  D       (input) REAL array, dimension (N)
043: *          The n diagonal elements of the diagonal matrix D from the
044: *          factorization A = U'*D*U or A = L*D*L'.
045: *
046: *  E       (input) COMPLEX array, dimension (N-1)
047: *          If IUPLO = 1, the (n-1) superdiagonal elements of the unit
048: *          bidiagonal factor U from the factorization A = U'*D*U.
049: *          If IUPLO = 0, the (n-1) subdiagonal elements of the unit
050: *          bidiagonal factor L from the factorization A = L*D*L'.
051: *
052: *  B       (input/output) REAL array, dimension (LDB,NRHS)
053: *          On entry, the right hand side vectors B for the system of
054: *          linear equations.
055: *          On exit, the solution vectors, X.
056: *
057: *  LDB     (input) INTEGER
058: *          The leading dimension of the array B.  LDB >= max(1,N).
059: *
060: *  =====================================================================
061: *
062: *     .. Local Scalars ..
063:       INTEGER            I, J
064: *     ..
065: *     .. External Subroutines ..
066:       EXTERNAL           CSSCAL
067: *     ..
068: *     .. Intrinsic Functions ..
069:       INTRINSIC          CONJG
070: *     ..
071: *     .. Executable Statements ..
072: *
073: *     Quick return if possible
074: *
075:       IF( N.LE.1 ) THEN
076:          IF( N.EQ.1 )
077:      $      CALL CSSCAL( NRHS, 1. / D( 1 ), B, LDB )
078:          RETURN
079:       END IF
080: *
081:       IF( IUPLO.EQ.1 ) THEN
082: *
083: *        Solve A * X = B using the factorization A = U'*D*U,
084: *        overwriting each right hand side vector with its solution.
085: *
086:          IF( NRHS.LE.2 ) THEN
087:             J = 1
088:     5       CONTINUE
089: *
090: *           Solve U' * x = b.
091: *
092:             DO 10 I = 2, N
093:                B( I, J ) = B( I, J ) - B( I-1, J )*CONJG( E( I-1 ) )
094:    10       CONTINUE
095: *
096: *           Solve D * U * x = b.
097: *
098:             DO 20 I = 1, N
099:                B( I, J ) = B( I, J ) / D( I )
100:    20       CONTINUE
101:             DO 30 I = N - 1, 1, -1
102:                B( I, J ) = B( I, J ) - B( I+1, J )*E( I )
103:    30       CONTINUE
104:             IF( J.LT.NRHS ) THEN
105:                J = J + 1
106:                GO TO 5
107:             END IF
108:          ELSE
109:             DO 60 J = 1, NRHS
110: *
111: *              Solve U' * x = b.
112: *
113:                DO 40 I = 2, N
114:                   B( I, J ) = B( I, J ) - B( I-1, J )*CONJG( E( I-1 ) )
115:    40          CONTINUE
116: *
117: *              Solve D * U * x = b.
118: *
119:                B( N, J ) = B( N, J ) / D( N )
120:                DO 50 I = N - 1, 1, -1
121:                   B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I )
122:    50          CONTINUE
123:    60       CONTINUE
124:          END IF
125:       ELSE
126: *
127: *        Solve A * X = B using the factorization A = L*D*L',
128: *        overwriting each right hand side vector with its solution.
129: *
130:          IF( NRHS.LE.2 ) THEN
131:             J = 1
132:    65       CONTINUE
133: *
134: *           Solve L * x = b.
135: *
136:             DO 70 I = 2, N
137:                B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
138:    70       CONTINUE
139: *
140: *           Solve D * L' * x = b.
141: *
142:             DO 80 I = 1, N
143:                B( I, J ) = B( I, J ) / D( I )
144:    80       CONTINUE
145:             DO 90 I = N - 1, 1, -1
146:                B( I, J ) = B( I, J ) - B( I+1, J )*CONJG( E( I ) )
147:    90       CONTINUE
148:             IF( J.LT.NRHS ) THEN
149:                J = J + 1
150:                GO TO 65
151:             END IF
152:          ELSE
153:             DO 120 J = 1, NRHS
154: *
155: *              Solve L * x = b.
156: *
157:                DO 100 I = 2, N
158:                   B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
159:   100          CONTINUE
160: *
161: *              Solve D * L' * x = b.
162: *
163:                B( N, J ) = B( N, J ) / D( N )
164:                DO 110 I = N - 1, 1, -1
165:                   B( I, J ) = B( I, J ) / D( I ) -
166:      $                        B( I+1, J )*CONJG( E( I ) )
167:   110          CONTINUE
168:   120       CONTINUE
169:          END IF
170:       END IF
171: *
172:       RETURN
173: *
174: *     End of CPTTS2
175: *
176:       END
177: