1 SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2 * .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA,BETA
4 INTEGER INCX,INCY,KL,KU,LDA,M,N
5 CHARACTER TRANS
6 * ..
7 * .. Array Arguments ..
8 DOUBLE PRECISION A(LDA,*),X(*),Y(*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DGBMV performs one of the matrix-vector operations
15 *
16 * y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y,
17 *
18 * where alpha and beta are scalars, x and y are vectors and A is an
19 * m by n band matrix, with kl sub-diagonals and ku super-diagonals.
20 *
21 * Arguments
22 * ==========
23 *
24 * TRANS - CHARACTER*1.
25 * On entry, TRANS specifies the operation to be performed as
26 * follows:
27 *
28 * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
29 *
30 * TRANS = 'T' or 't' y := alpha*A**T*x + beta*y.
31 *
32 * TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y.
33 *
34 * Unchanged on exit.
35 *
36 * M - INTEGER.
37 * On entry, M specifies the number of rows of the matrix A.
38 * M must be at least zero.
39 * Unchanged on exit.
40 *
41 * N - INTEGER.
42 * On entry, N specifies the number of columns of the matrix A.
43 * N must be at least zero.
44 * Unchanged on exit.
45 *
46 * KL - INTEGER.
47 * On entry, KL specifies the number of sub-diagonals of the
48 * matrix A. KL must satisfy 0 .le. KL.
49 * Unchanged on exit.
50 *
51 * KU - INTEGER.
52 * On entry, KU specifies the number of super-diagonals of the
53 * matrix A. KU must satisfy 0 .le. KU.
54 * Unchanged on exit.
55 *
56 * ALPHA - DOUBLE PRECISION.
57 * On entry, ALPHA specifies the scalar alpha.
58 * Unchanged on exit.
59 *
60 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
61 * Before entry, the leading ( kl + ku + 1 ) by n part of the
62 * array A must contain the matrix of coefficients, supplied
63 * column by column, with the leading diagonal of the matrix in
64 * row ( ku + 1 ) of the array, the first super-diagonal
65 * starting at position 2 in row ku, the first sub-diagonal
66 * starting at position 1 in row ( ku + 2 ), and so on.
67 * Elements in the array A that do not correspond to elements
68 * in the band matrix (such as the top left ku by ku triangle)
69 * are not referenced.
70 * The following program segment will transfer a band matrix
71 * from conventional full matrix storage to band storage:
72 *
73 * DO 20, J = 1, N
74 * K = KU + 1 - J
75 * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
76 * A( K + I, J ) = matrix( I, J )
77 * 10 CONTINUE
78 * 20 CONTINUE
79 *
80 * Unchanged on exit.
81 *
82 * LDA - INTEGER.
83 * On entry, LDA specifies the first dimension of A as declared
84 * in the calling (sub) program. LDA must be at least
85 * ( kl + ku + 1 ).
86 * Unchanged on exit.
87 *
88 * X - DOUBLE PRECISION array of DIMENSION at least
89 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
90 * and at least
91 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
92 * Before entry, the incremented array X must contain the
93 * vector x.
94 * Unchanged on exit.
95 *
96 * INCX - INTEGER.
97 * On entry, INCX specifies the increment for the elements of
98 * X. INCX must not be zero.
99 * Unchanged on exit.
100 *
101 * BETA - DOUBLE PRECISION.
102 * On entry, BETA specifies the scalar beta. When BETA is
103 * supplied as zero then Y need not be set on input.
104 * Unchanged on exit.
105 *
106 * Y - DOUBLE PRECISION array of DIMENSION at least
107 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
108 * and at least
109 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
110 * Before entry, the incremented array Y must contain the
111 * vector y. On exit, Y is overwritten by the updated vector y.
112 *
113 * INCY - INTEGER.
114 * On entry, INCY specifies the increment for the elements of
115 * Y. INCY must not be zero.
116 * Unchanged on exit.
117 *
118 * Further Details
119 * ===============
120 *
121 * Level 2 Blas routine.
122 * The vector and matrix arguments are not referenced when N = 0, or M = 0
123 *
124 * -- Written on 22-October-1986.
125 * Jack Dongarra, Argonne National Lab.
126 * Jeremy Du Croz, Nag Central Office.
127 * Sven Hammarling, Nag Central Office.
128 * Richard Hanson, Sandia National Labs.
129 *
130 * =====================================================================
131 *
132 * .. Parameters ..
133 DOUBLE PRECISION ONE,ZERO
134 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
135 * ..
136 * .. Local Scalars ..
137 DOUBLE PRECISION TEMP
138 INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
139 * ..
140 * .. External Functions ..
141 LOGICAL LSAME
142 EXTERNAL LSAME
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC MAX,MIN
149 * ..
150 *
151 * Test the input parameters.
152 *
153 INFO = 0
154 IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
155 + .NOT.LSAME(TRANS,'C')) THEN
156 INFO = 1
157 ELSE IF (M.LT.0) THEN
158 INFO = 2
159 ELSE IF (N.LT.0) THEN
160 INFO = 3
161 ELSE IF (KL.LT.0) THEN
162 INFO = 4
163 ELSE IF (KU.LT.0) THEN
164 INFO = 5
165 ELSE IF (LDA.LT. (KL+KU+1)) THEN
166 INFO = 8
167 ELSE IF (INCX.EQ.0) THEN
168 INFO = 10
169 ELSE IF (INCY.EQ.0) THEN
170 INFO = 13
171 END IF
172 IF (INFO.NE.0) THEN
173 CALL XERBLA('DGBMV ',INFO)
174 RETURN
175 END IF
176 *
177 * Quick return if possible.
178 *
179 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
180 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
181 *
182 * Set LENX and LENY, the lengths of the vectors x and y, and set
183 * up the start points in X and Y.
184 *
185 IF (LSAME(TRANS,'N')) THEN
186 LENX = N
187 LENY = M
188 ELSE
189 LENX = M
190 LENY = N
191 END IF
192 IF (INCX.GT.0) THEN
193 KX = 1
194 ELSE
195 KX = 1 - (LENX-1)*INCX
196 END IF
197 IF (INCY.GT.0) THEN
198 KY = 1
199 ELSE
200 KY = 1 - (LENY-1)*INCY
201 END IF
202 *
203 * Start the operations. In this version the elements of A are
204 * accessed sequentially with one pass through the band part of A.
205 *
206 * First form y := beta*y.
207 *
208 IF (BETA.NE.ONE) THEN
209 IF (INCY.EQ.1) THEN
210 IF (BETA.EQ.ZERO) THEN
211 DO 10 I = 1,LENY
212 Y(I) = ZERO
213 10 CONTINUE
214 ELSE
215 DO 20 I = 1,LENY
216 Y(I) = BETA*Y(I)
217 20 CONTINUE
218 END IF
219 ELSE
220 IY = KY
221 IF (BETA.EQ.ZERO) THEN
222 DO 30 I = 1,LENY
223 Y(IY) = ZERO
224 IY = IY + INCY
225 30 CONTINUE
226 ELSE
227 DO 40 I = 1,LENY
228 Y(IY) = BETA*Y(IY)
229 IY = IY + INCY
230 40 CONTINUE
231 END IF
232 END IF
233 END IF
234 IF (ALPHA.EQ.ZERO) RETURN
235 KUP1 = KU + 1
236 IF (LSAME(TRANS,'N')) THEN
237 *
238 * Form y := alpha*A*x + y.
239 *
240 JX = KX
241 IF (INCY.EQ.1) THEN
242 DO 60 J = 1,N
243 IF (X(JX).NE.ZERO) THEN
244 TEMP = ALPHA*X(JX)
245 K = KUP1 - J
246 DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
247 Y(I) = Y(I) + TEMP*A(K+I,J)
248 50 CONTINUE
249 END IF
250 JX = JX + INCX
251 60 CONTINUE
252 ELSE
253 DO 80 J = 1,N
254 IF (X(JX).NE.ZERO) THEN
255 TEMP = ALPHA*X(JX)
256 IY = KY
257 K = KUP1 - J
258 DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
259 Y(IY) = Y(IY) + TEMP*A(K+I,J)
260 IY = IY + INCY
261 70 CONTINUE
262 END IF
263 JX = JX + INCX
264 IF (J.GT.KU) KY = KY + INCY
265 80 CONTINUE
266 END IF
267 ELSE
268 *
269 * Form y := alpha*A**T*x + y.
270 *
271 JY = KY
272 IF (INCX.EQ.1) THEN
273 DO 100 J = 1,N
274 TEMP = ZERO
275 K = KUP1 - J
276 DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
277 TEMP = TEMP + A(K+I,J)*X(I)
278 90 CONTINUE
279 Y(JY) = Y(JY) + ALPHA*TEMP
280 JY = JY + INCY
281 100 CONTINUE
282 ELSE
283 DO 120 J = 1,N
284 TEMP = ZERO
285 IX = KX
286 K = KUP1 - J
287 DO 110 I = MAX(1,J-KU),MIN(M,J+KL)
288 TEMP = TEMP + A(K+I,J)*X(IX)
289 IX = IX + INCX
290 110 CONTINUE
291 Y(JY) = Y(JY) + ALPHA*TEMP
292 JY = JY + INCY
293 IF (J.GT.KU) KX = KX + INCX
294 120 CONTINUE
295 END IF
296 END IF
297 *
298 RETURN
299 *
300 * End of DGBMV .
301 *
302 END
2 * .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA,BETA
4 INTEGER INCX,INCY,KL,KU,LDA,M,N
5 CHARACTER TRANS
6 * ..
7 * .. Array Arguments ..
8 DOUBLE PRECISION A(LDA,*),X(*),Y(*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DGBMV performs one of the matrix-vector operations
15 *
16 * y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y,
17 *
18 * where alpha and beta are scalars, x and y are vectors and A is an
19 * m by n band matrix, with kl sub-diagonals and ku super-diagonals.
20 *
21 * Arguments
22 * ==========
23 *
24 * TRANS - CHARACTER*1.
25 * On entry, TRANS specifies the operation to be performed as
26 * follows:
27 *
28 * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
29 *
30 * TRANS = 'T' or 't' y := alpha*A**T*x + beta*y.
31 *
32 * TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y.
33 *
34 * Unchanged on exit.
35 *
36 * M - INTEGER.
37 * On entry, M specifies the number of rows of the matrix A.
38 * M must be at least zero.
39 * Unchanged on exit.
40 *
41 * N - INTEGER.
42 * On entry, N specifies the number of columns of the matrix A.
43 * N must be at least zero.
44 * Unchanged on exit.
45 *
46 * KL - INTEGER.
47 * On entry, KL specifies the number of sub-diagonals of the
48 * matrix A. KL must satisfy 0 .le. KL.
49 * Unchanged on exit.
50 *
51 * KU - INTEGER.
52 * On entry, KU specifies the number of super-diagonals of the
53 * matrix A. KU must satisfy 0 .le. KU.
54 * Unchanged on exit.
55 *
56 * ALPHA - DOUBLE PRECISION.
57 * On entry, ALPHA specifies the scalar alpha.
58 * Unchanged on exit.
59 *
60 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
61 * Before entry, the leading ( kl + ku + 1 ) by n part of the
62 * array A must contain the matrix of coefficients, supplied
63 * column by column, with the leading diagonal of the matrix in
64 * row ( ku + 1 ) of the array, the first super-diagonal
65 * starting at position 2 in row ku, the first sub-diagonal
66 * starting at position 1 in row ( ku + 2 ), and so on.
67 * Elements in the array A that do not correspond to elements
68 * in the band matrix (such as the top left ku by ku triangle)
69 * are not referenced.
70 * The following program segment will transfer a band matrix
71 * from conventional full matrix storage to band storage:
72 *
73 * DO 20, J = 1, N
74 * K = KU + 1 - J
75 * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
76 * A( K + I, J ) = matrix( I, J )
77 * 10 CONTINUE
78 * 20 CONTINUE
79 *
80 * Unchanged on exit.
81 *
82 * LDA - INTEGER.
83 * On entry, LDA specifies the first dimension of A as declared
84 * in the calling (sub) program. LDA must be at least
85 * ( kl + ku + 1 ).
86 * Unchanged on exit.
87 *
88 * X - DOUBLE PRECISION array of DIMENSION at least
89 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
90 * and at least
91 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
92 * Before entry, the incremented array X must contain the
93 * vector x.
94 * Unchanged on exit.
95 *
96 * INCX - INTEGER.
97 * On entry, INCX specifies the increment for the elements of
98 * X. INCX must not be zero.
99 * Unchanged on exit.
100 *
101 * BETA - DOUBLE PRECISION.
102 * On entry, BETA specifies the scalar beta. When BETA is
103 * supplied as zero then Y need not be set on input.
104 * Unchanged on exit.
105 *
106 * Y - DOUBLE PRECISION array of DIMENSION at least
107 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
108 * and at least
109 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
110 * Before entry, the incremented array Y must contain the
111 * vector y. On exit, Y is overwritten by the updated vector y.
112 *
113 * INCY - INTEGER.
114 * On entry, INCY specifies the increment for the elements of
115 * Y. INCY must not be zero.
116 * Unchanged on exit.
117 *
118 * Further Details
119 * ===============
120 *
121 * Level 2 Blas routine.
122 * The vector and matrix arguments are not referenced when N = 0, or M = 0
123 *
124 * -- Written on 22-October-1986.
125 * Jack Dongarra, Argonne National Lab.
126 * Jeremy Du Croz, Nag Central Office.
127 * Sven Hammarling, Nag Central Office.
128 * Richard Hanson, Sandia National Labs.
129 *
130 * =====================================================================
131 *
132 * .. Parameters ..
133 DOUBLE PRECISION ONE,ZERO
134 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
135 * ..
136 * .. Local Scalars ..
137 DOUBLE PRECISION TEMP
138 INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
139 * ..
140 * .. External Functions ..
141 LOGICAL LSAME
142 EXTERNAL LSAME
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC MAX,MIN
149 * ..
150 *
151 * Test the input parameters.
152 *
153 INFO = 0
154 IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
155 + .NOT.LSAME(TRANS,'C')) THEN
156 INFO = 1
157 ELSE IF (M.LT.0) THEN
158 INFO = 2
159 ELSE IF (N.LT.0) THEN
160 INFO = 3
161 ELSE IF (KL.LT.0) THEN
162 INFO = 4
163 ELSE IF (KU.LT.0) THEN
164 INFO = 5
165 ELSE IF (LDA.LT. (KL+KU+1)) THEN
166 INFO = 8
167 ELSE IF (INCX.EQ.0) THEN
168 INFO = 10
169 ELSE IF (INCY.EQ.0) THEN
170 INFO = 13
171 END IF
172 IF (INFO.NE.0) THEN
173 CALL XERBLA('DGBMV ',INFO)
174 RETURN
175 END IF
176 *
177 * Quick return if possible.
178 *
179 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
180 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
181 *
182 * Set LENX and LENY, the lengths of the vectors x and y, and set
183 * up the start points in X and Y.
184 *
185 IF (LSAME(TRANS,'N')) THEN
186 LENX = N
187 LENY = M
188 ELSE
189 LENX = M
190 LENY = N
191 END IF
192 IF (INCX.GT.0) THEN
193 KX = 1
194 ELSE
195 KX = 1 - (LENX-1)*INCX
196 END IF
197 IF (INCY.GT.0) THEN
198 KY = 1
199 ELSE
200 KY = 1 - (LENY-1)*INCY
201 END IF
202 *
203 * Start the operations. In this version the elements of A are
204 * accessed sequentially with one pass through the band part of A.
205 *
206 * First form y := beta*y.
207 *
208 IF (BETA.NE.ONE) THEN
209 IF (INCY.EQ.1) THEN
210 IF (BETA.EQ.ZERO) THEN
211 DO 10 I = 1,LENY
212 Y(I) = ZERO
213 10 CONTINUE
214 ELSE
215 DO 20 I = 1,LENY
216 Y(I) = BETA*Y(I)
217 20 CONTINUE
218 END IF
219 ELSE
220 IY = KY
221 IF (BETA.EQ.ZERO) THEN
222 DO 30 I = 1,LENY
223 Y(IY) = ZERO
224 IY = IY + INCY
225 30 CONTINUE
226 ELSE
227 DO 40 I = 1,LENY
228 Y(IY) = BETA*Y(IY)
229 IY = IY + INCY
230 40 CONTINUE
231 END IF
232 END IF
233 END IF
234 IF (ALPHA.EQ.ZERO) RETURN
235 KUP1 = KU + 1
236 IF (LSAME(TRANS,'N')) THEN
237 *
238 * Form y := alpha*A*x + y.
239 *
240 JX = KX
241 IF (INCY.EQ.1) THEN
242 DO 60 J = 1,N
243 IF (X(JX).NE.ZERO) THEN
244 TEMP = ALPHA*X(JX)
245 K = KUP1 - J
246 DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
247 Y(I) = Y(I) + TEMP*A(K+I,J)
248 50 CONTINUE
249 END IF
250 JX = JX + INCX
251 60 CONTINUE
252 ELSE
253 DO 80 J = 1,N
254 IF (X(JX).NE.ZERO) THEN
255 TEMP = ALPHA*X(JX)
256 IY = KY
257 K = KUP1 - J
258 DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
259 Y(IY) = Y(IY) + TEMP*A(K+I,J)
260 IY = IY + INCY
261 70 CONTINUE
262 END IF
263 JX = JX + INCX
264 IF (J.GT.KU) KY = KY + INCY
265 80 CONTINUE
266 END IF
267 ELSE
268 *
269 * Form y := alpha*A**T*x + y.
270 *
271 JY = KY
272 IF (INCX.EQ.1) THEN
273 DO 100 J = 1,N
274 TEMP = ZERO
275 K = KUP1 - J
276 DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
277 TEMP = TEMP + A(K+I,J)*X(I)
278 90 CONTINUE
279 Y(JY) = Y(JY) + ALPHA*TEMP
280 JY = JY + INCY
281 100 CONTINUE
282 ELSE
283 DO 120 J = 1,N
284 TEMP = ZERO
285 IX = KX
286 K = KUP1 - J
287 DO 110 I = MAX(1,J-KU),MIN(M,J+KL)
288 TEMP = TEMP + A(K+I,J)*X(IX)
289 IX = IX + INCX
290 110 CONTINUE
291 Y(JY) = Y(JY) + ALPHA*TEMP
292 JY = JY + INCY
293 IF (J.GT.KU) KX = KX + INCX
294 120 CONTINUE
295 END IF
296 END IF
297 *
298 RETURN
299 *
300 * End of DGBMV .
301 *
302 END