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