001:       SUBROUTINE ZTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
002: *     .. Scalar Arguments ..
003:       INTEGER INCX,K,LDA,N
004:       CHARACTER DIAG,TRANS,UPLO
005: *     ..
006: *     .. Array Arguments ..
007:       DOUBLE COMPLEX A(LDA,*),X(*)
008: *     ..
009: *
010: *  Purpose
011: *  =======
012: *
013: *  ZTBMV  performs one of the matrix-vector operations
014: *
015: *     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x,
016: *
017: *  where x is an n element vector and  A is an n by n unit, or non-unit,
018: *  upper or lower triangular band matrix, with ( k + 1 ) diagonals.
019: *
020: *  Arguments
021: *  ==========
022: *
023: *  UPLO   - CHARACTER*1.
024: *           On entry, UPLO specifies whether the matrix is an upper or
025: *           lower triangular matrix as follows:
026: *
027: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
028: *
029: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
030: *
031: *           Unchanged on exit.
032: *
033: *  TRANS  - CHARACTER*1.
034: *           On entry, TRANS specifies the operation to be performed as
035: *           follows:
036: *
037: *              TRANS = 'N' or 'n'   x := A*x.
038: *
039: *              TRANS = 'T' or 't'   x := A'*x.
040: *
041: *              TRANS = 'C' or 'c'   x := conjg( A' )*x.
042: *
043: *           Unchanged on exit.
044: *
045: *  DIAG   - CHARACTER*1.
046: *           On entry, DIAG specifies whether or not A is unit
047: *           triangular as follows:
048: *
049: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
050: *
051: *              DIAG = 'N' or 'n'   A is not assumed to be unit
052: *                                  triangular.
053: *
054: *           Unchanged on exit.
055: *
056: *  N      - INTEGER.
057: *           On entry, N specifies the order of the matrix A.
058: *           N must be at least zero.
059: *           Unchanged on exit.
060: *
061: *  K      - INTEGER.
062: *           On entry with UPLO = 'U' or 'u', K specifies the number of
063: *           super-diagonals of the matrix A.
064: *           On entry with UPLO = 'L' or 'l', K specifies the number of
065: *           sub-diagonals of the matrix A.
066: *           K must satisfy  0 .le. K.
067: *           Unchanged on exit.
068: *
069: *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
070: *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
071: *           by n part of the array A must contain the upper triangular
072: *           band part of the matrix of coefficients, supplied column by
073: *           column, with the leading diagonal of the matrix in row
074: *           ( k + 1 ) of the array, the first super-diagonal starting at
075: *           position 2 in row k, and so on. The top left k by k triangle
076: *           of the array A is not referenced.
077: *           The following program segment will transfer an upper
078: *           triangular band matrix from conventional full matrix storage
079: *           to band storage:
080: *
081: *                 DO 20, J = 1, N
082: *                    M = K + 1 - J
083: *                    DO 10, I = MAX( 1, J - K ), J
084: *                       A( M + I, J ) = matrix( I, J )
085: *              10    CONTINUE
086: *              20 CONTINUE
087: *
088: *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
089: *           by n part of the array A must contain the lower triangular
090: *           band part of the matrix of coefficients, supplied column by
091: *           column, with the leading diagonal of the matrix in row 1 of
092: *           the array, the first sub-diagonal starting at position 1 in
093: *           row 2, and so on. The bottom right k by k triangle of the
094: *           array A is not referenced.
095: *           The following program segment will transfer a lower
096: *           triangular band matrix from conventional full matrix storage
097: *           to band storage:
098: *
099: *                 DO 20, J = 1, N
100: *                    M = 1 - J
101: *                    DO 10, I = J, MIN( N, J + K )
102: *                       A( M + I, J ) = matrix( I, J )
103: *              10    CONTINUE
104: *              20 CONTINUE
105: *
106: *           Note that when DIAG = 'U' or 'u' the elements of the array A
107: *           corresponding to the diagonal elements of the matrix are not
108: *           referenced, but are assumed to be unity.
109: *           Unchanged on exit.
110: *
111: *  LDA    - INTEGER.
112: *           On entry, LDA specifies the first dimension of A as declared
113: *           in the calling (sub) program. LDA must be at least
114: *           ( k + 1 ).
115: *           Unchanged on exit.
116: *
117: *  X      - COMPLEX*16       array of dimension at least
118: *           ( 1 + ( n - 1 )*abs( INCX ) ).
119: *           Before entry, the incremented array X must contain the n
120: *           element vector x. On exit, X is overwritten with the
121: *           tranformed vector x.
122: *
123: *  INCX   - INTEGER.
124: *           On entry, INCX specifies the increment for the elements of
125: *           X. INCX must not be zero.
126: *           Unchanged on exit.
127: *
128: *  Further Details
129: *  ===============
130: *
131: *  Level 2 Blas routine.
132: *
133: *  -- Written on 22-October-1986.
134: *     Jack Dongarra, Argonne National Lab.
135: *     Jeremy Du Croz, Nag Central Office.
136: *     Sven Hammarling, Nag Central Office.
137: *     Richard Hanson, Sandia National Labs.
138: *
139: *  =====================================================================
140: *
141: *     .. Parameters ..
142:       DOUBLE COMPLEX ZERO
143:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
144: *     ..
145: *     .. Local Scalars ..
146:       DOUBLE COMPLEX TEMP
147:       INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
148:       LOGICAL NOCONJ,NOUNIT
149: *     ..
150: *     .. External Functions ..
151:       LOGICAL LSAME
152:       EXTERNAL LSAME
153: *     ..
154: *     .. External Subroutines ..
155:       EXTERNAL XERBLA
156: *     ..
157: *     .. Intrinsic Functions ..
158:       INTRINSIC DCONJG,MAX,MIN
159: *     ..
160: *
161: *     Test the input parameters.
162: *
163:       INFO = 0
164:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
165:           INFO = 1
166:       ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
167:      +         .NOT.LSAME(TRANS,'C')) THEN
168:           INFO = 2
169:       ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
170:           INFO = 3
171:       ELSE IF (N.LT.0) THEN
172:           INFO = 4
173:       ELSE IF (K.LT.0) THEN
174:           INFO = 5
175:       ELSE IF (LDA.LT. (K+1)) THEN
176:           INFO = 7
177:       ELSE IF (INCX.EQ.0) THEN
178:           INFO = 9
179:       END IF
180:       IF (INFO.NE.0) THEN
181:           CALL XERBLA('ZTBMV ',INFO)
182:           RETURN
183:       END IF
184: *
185: *     Quick return if possible.
186: *
187:       IF (N.EQ.0) RETURN
188: *
189:       NOCONJ = LSAME(TRANS,'T')
190:       NOUNIT = LSAME(DIAG,'N')
191: *
192: *     Set up the start point in X if the increment is not unity. This
193: *     will be  ( N - 1 )*INCX   too small for descending loops.
194: *
195:       IF (INCX.LE.0) THEN
196:           KX = 1 - (N-1)*INCX
197:       ELSE IF (INCX.NE.1) THEN
198:           KX = 1
199:       END IF
200: *
201: *     Start the operations. In this version the elements of A are
202: *     accessed sequentially with one pass through A.
203: *
204:       IF (LSAME(TRANS,'N')) THEN
205: *
206: *         Form  x := A*x.
207: *
208:           IF (LSAME(UPLO,'U')) THEN
209:               KPLUS1 = K + 1
210:               IF (INCX.EQ.1) THEN
211:                   DO 20 J = 1,N
212:                       IF (X(J).NE.ZERO) THEN
213:                           TEMP = X(J)
214:                           L = KPLUS1 - J
215:                           DO 10 I = MAX(1,J-K),J - 1
216:                               X(I) = X(I) + TEMP*A(L+I,J)
217:    10                     CONTINUE
218:                           IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
219:                       END IF
220:    20             CONTINUE
221:               ELSE
222:                   JX = KX
223:                   DO 40 J = 1,N
224:                       IF (X(JX).NE.ZERO) THEN
225:                           TEMP = X(JX)
226:                           IX = KX
227:                           L = KPLUS1 - J
228:                           DO 30 I = MAX(1,J-K),J - 1
229:                               X(IX) = X(IX) + TEMP*A(L+I,J)
230:                               IX = IX + INCX
231:    30                     CONTINUE
232:                           IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
233:                       END IF
234:                       JX = JX + INCX
235:                       IF (J.GT.K) KX = KX + INCX
236:    40             CONTINUE
237:               END IF
238:           ELSE
239:               IF (INCX.EQ.1) THEN
240:                   DO 60 J = N,1,-1
241:                       IF (X(J).NE.ZERO) THEN
242:                           TEMP = X(J)
243:                           L = 1 - J
244:                           DO 50 I = MIN(N,J+K),J + 1,-1
245:                               X(I) = X(I) + TEMP*A(L+I,J)
246:    50                     CONTINUE
247:                           IF (NOUNIT) X(J) = X(J)*A(1,J)
248:                       END IF
249:    60             CONTINUE
250:               ELSE
251:                   KX = KX + (N-1)*INCX
252:                   JX = KX
253:                   DO 80 J = N,1,-1
254:                       IF (X(JX).NE.ZERO) THEN
255:                           TEMP = X(JX)
256:                           IX = KX
257:                           L = 1 - J
258:                           DO 70 I = MIN(N,J+K),J + 1,-1
259:                               X(IX) = X(IX) + TEMP*A(L+I,J)
260:                               IX = IX - INCX
261:    70                     CONTINUE
262:                           IF (NOUNIT) X(JX) = X(JX)*A(1,J)
263:                       END IF
264:                       JX = JX - INCX
265:                       IF ((N-J).GE.K) KX = KX - INCX
266:    80             CONTINUE
267:               END IF
268:           END IF
269:       ELSE
270: *
271: *        Form  x := A'*x  or  x := conjg( A' )*x.
272: *
273:           IF (LSAME(UPLO,'U')) THEN
274:               KPLUS1 = K + 1
275:               IF (INCX.EQ.1) THEN
276:                   DO 110 J = N,1,-1
277:                       TEMP = X(J)
278:                       L = KPLUS1 - J
279:                       IF (NOCONJ) THEN
280:                           IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
281:                           DO 90 I = J - 1,MAX(1,J-K),-1
282:                               TEMP = TEMP + A(L+I,J)*X(I)
283:    90                     CONTINUE
284:                       ELSE
285:                           IF (NOUNIT) TEMP = TEMP*DCONJG(A(KPLUS1,J))
286:                           DO 100 I = J - 1,MAX(1,J-K),-1
287:                               TEMP = TEMP + DCONJG(A(L+I,J))*X(I)
288:   100                     CONTINUE
289:                       END IF
290:                       X(J) = TEMP
291:   110             CONTINUE
292:               ELSE
293:                   KX = KX + (N-1)*INCX
294:                   JX = KX
295:                   DO 140 J = N,1,-1
296:                       TEMP = X(JX)
297:                       KX = KX - INCX
298:                       IX = KX
299:                       L = KPLUS1 - J
300:                       IF (NOCONJ) THEN
301:                           IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
302:                           DO 120 I = J - 1,MAX(1,J-K),-1
303:                               TEMP = TEMP + A(L+I,J)*X(IX)
304:                               IX = IX - INCX
305:   120                     CONTINUE
306:                       ELSE
307:                           IF (NOUNIT) TEMP = TEMP*DCONJG(A(KPLUS1,J))
308:                           DO 130 I = J - 1,MAX(1,J-K),-1
309:                               TEMP = TEMP + DCONJG(A(L+I,J))*X(IX)
310:                               IX = IX - INCX
311:   130                     CONTINUE
312:                       END IF
313:                       X(JX) = TEMP
314:                       JX = JX - INCX
315:   140             CONTINUE
316:               END IF
317:           ELSE
318:               IF (INCX.EQ.1) THEN
319:                   DO 170 J = 1,N
320:                       TEMP = X(J)
321:                       L = 1 - J
322:                       IF (NOCONJ) THEN
323:                           IF (NOUNIT) TEMP = TEMP*A(1,J)
324:                           DO 150 I = J + 1,MIN(N,J+K)
325:                               TEMP = TEMP + A(L+I,J)*X(I)
326:   150                     CONTINUE
327:                       ELSE
328:                           IF (NOUNIT) TEMP = TEMP*DCONJG(A(1,J))
329:                           DO 160 I = J + 1,MIN(N,J+K)
330:                               TEMP = TEMP + DCONJG(A(L+I,J))*X(I)
331:   160                     CONTINUE
332:                       END IF
333:                       X(J) = TEMP
334:   170             CONTINUE
335:               ELSE
336:                   JX = KX
337:                   DO 200 J = 1,N
338:                       TEMP = X(JX)
339:                       KX = KX + INCX
340:                       IX = KX
341:                       L = 1 - J
342:                       IF (NOCONJ) THEN
343:                           IF (NOUNIT) TEMP = TEMP*A(1,J)
344:                           DO 180 I = J + 1,MIN(N,J+K)
345:                               TEMP = TEMP + A(L+I,J)*X(IX)
346:                               IX = IX + INCX
347:   180                     CONTINUE
348:                       ELSE
349:                           IF (NOUNIT) TEMP = TEMP*DCONJG(A(1,J))
350:                           DO 190 I = J + 1,MIN(N,J+K)
351:                               TEMP = TEMP + DCONJG(A(L+I,J))*X(IX)
352:                               IX = IX + INCX
353:   190                     CONTINUE
354:                       END IF
355:                       X(JX) = TEMP
356:                       JX = JX + INCX
357:   200             CONTINUE
358:               END IF
359:           END IF
360:       END IF
361: *
362:       RETURN
363: *
364: *     End of ZTBMV .
365: *
366:       END
367: