1 /* sgbmv.f -- translated by f2c (version 19991025). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */
2 #include "FLA_f2c.h"
3 /* Subroutine */
sgbmv_(char * trans,integer * m,integer * n,integer * kl,integer * ku,real * alpha,real * a,integer * lda,real * x,integer * incx,real * beta,real * y,integer * incy)4 int sgbmv_(char *trans, integer *m, integer *n, integer *kl, integer *ku, real *alpha, real *a, integer *lda, real *x, integer * incx, real *beta, real *y, integer *incy)
5 {
6 /* System generated locals */
7 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
8 /* Local variables */
9 integer info;
10 real temp;
11 integer lenx, leny, i__, j, k;
12 extern logical lsame_(char *, char *);
13 integer ix, iy, jx, jy, kx, ky;
14 extern /* Subroutine */
15 int xerbla_(char *, integer *);
16 integer kup1;
17 /* .. Scalar Arguments .. */
18 /* .. Array Arguments .. */
19 /* .. */
20 /* Purpose */
21 /* ======= */
22 /* SGBMV performs one of the matrix-vector operations */
23 /* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
24 /* where alpha and beta are scalars, x and y are vectors and A is an */
25 /* m by n band matrix, with kl sub-diagonals and ku super-diagonals. */
26 /* Parameters */
27 /* ========== */
28 /* TRANS - CHARACTER*1. */
29 /* On entry, TRANS specifies the operation to be performed as */
30 /* follows: */
31 /* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
32 /* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
33 /* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
34 /* Unchanged on exit. */
35 /* M - INTEGER. */
36 /* On entry, M specifies the number of rows of the matrix A. */
37 /* M must be at least zero. */
38 /* Unchanged on exit. */
39 /* N - INTEGER. */
40 /* On entry, N specifies the number of columns of the matrix A. */
41 /* N must be at least zero. */
42 /* Unchanged on exit. */
43 /* KL - INTEGER. */
44 /* On entry, KL specifies the number of sub-diagonals of the */
45 /* matrix A. KL must satisfy 0 .le. KL. */
46 /* Unchanged on exit. */
47 /* KU - INTEGER. */
48 /* On entry, KU specifies the number of super-diagonals of the */
49 /* matrix A. KU must satisfy 0 .le. KU. */
50 /* Unchanged on exit. */
51 /* ALPHA - REAL . */
52 /* On entry, ALPHA specifies the scalar alpha. */
53 /* Unchanged on exit. */
54 /* A - REAL array of DIMENSION ( LDA, n ). */
55 /* Before entry, the leading ( kl + ku + 1 ) by n part of the */
56 /* array A must contain the matrix of coefficients, supplied */
57 /* column by column, with the leading diagonal of the matrix in */
58 /* row ( ku + 1 ) of the array, the first super-diagonal */
59 /* starting at position 2 in row ku, the first sub-diagonal */
60 /* starting at position 1 in row ( ku + 2 ), and so on. */
61 /* Elements in the array A that do not correspond to elements */
62 /* in the band matrix (such as the top left ku by ku triangle) */
63 /* are not referenced. */
64 /* The following program segment will transfer a band matrix */
65 /* from conventional full matrix storage to band storage: */
66 /* DO 20, J = 1, N */
67 /* K = KU + 1 - J */
68 /* DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) */
69 /* A( K + I, J ) = matrix( I, J ) */
70 /* 10 CONTINUE */
71 /* 20 CONTINUE */
72 /* Unchanged on exit. */
73 /* LDA - INTEGER. */
74 /* On entry, LDA specifies the first dimension of A as declared */
75 /* in the calling (sub) program. LDA must be at least */
76 /* ( kl + ku + 1 ). */
77 /* Unchanged on exit. */
78 /* X - REAL array of DIMENSION at least */
79 /* ( 1 + ( n - 1 )*f2c_abs( INCX ) ) when TRANS = 'N' or 'n' */
80 /* and at least */
81 /* ( 1 + ( m - 1 )*f2c_abs( INCX ) ) otherwise. */
82 /* Before entry, the incremented array X must contain the */
83 /* vector x. */
84 /* Unchanged on exit. */
85 /* INCX - INTEGER. */
86 /* On entry, INCX specifies the increment for the elements of */
87 /* X. INCX must not be zero. */
88 /* Unchanged on exit. */
89 /* BETA - REAL . */
90 /* On entry, BETA specifies the scalar beta. When BETA is */
91 /* supplied as zero then Y need not be set on input. */
92 /* Unchanged on exit. */
93 /* Y - REAL array of DIMENSION at least */
94 /* ( 1 + ( m - 1 )*f2c_abs( INCY ) ) when TRANS = 'N' or 'n' */
95 /* and at least */
96 /* ( 1 + ( n - 1 )*f2c_abs( INCY ) ) otherwise. */
97 /* Before entry, the incremented array Y must contain the */
98 /* vector y. On exit, Y is overwritten by the updated vector y. */
99 /* INCY - INTEGER. */
100 /* On entry, INCY specifies the increment for the elements of */
101 /* Y. INCY must not be zero. */
102 /* Unchanged on exit. */
103 /* Level 2 Blas routine. */
104 /* -- Written on 22-October-1986. */
105 /* Jack Dongarra, Argonne National Lab. */
106 /* Jeremy Du Croz, Nag Central Office. */
107 /* Sven Hammarling, Nag Central Office. */
108 /* Richard Hanson, Sandia National Labs. */
109 /* .. Parameters .. */
110 /* .. Local Scalars .. */
111 /* .. External Functions .. */
112 /* .. External Subroutines .. */
113 /* .. Intrinsic Functions .. */
114 /* .. */
115 /* .. Executable Statements .. */
116 /* Test the input parameters. */
117 /* Parameter adjustments */
118 a_dim1 = *lda;
119 a_offset = 1 + a_dim1 * 1;
120 a -= a_offset;
121 --x;
122 --y;
123 /* Function Body */
124 info = 0;
125 if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C") )
126 {
127 info = 1;
128 }
129 else if (*m < 0)
130 {
131 info = 2;
132 }
133 else if (*n < 0)
134 {
135 info = 3;
136 }
137 else if (*kl < 0)
138 {
139 info = 4;
140 }
141 else if (*ku < 0)
142 {
143 info = 5;
144 }
145 else if (*lda < *kl + *ku + 1)
146 {
147 info = 8;
148 }
149 else if (*incx == 0)
150 {
151 info = 10;
152 }
153 else if (*incy == 0)
154 {
155 info = 13;
156 }
157 if (info != 0)
158 {
159 xerbla_("SGBMV ", &info);
160 return 0;
161 }
162 /* Quick return if possible. */
163 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f)
164 {
165 return 0;
166 }
167 /* Set LENX and LENY, the lengths of the vectors x and y, and set */
168 /* up the start points in X and Y. */
169 if (lsame_(trans, "N"))
170 {
171 lenx = *n;
172 leny = *m;
173 }
174 else
175 {
176 lenx = *m;
177 leny = *n;
178 }
179 if (*incx > 0)
180 {
181 kx = 1;
182 }
183 else
184 {
185 kx = 1 - (lenx - 1) * *incx;
186 }
187 if (*incy > 0)
188 {
189 ky = 1;
190 }
191 else
192 {
193 ky = 1 - (leny - 1) * *incy;
194 }
195 /* Start the operations. In this version the elements of A are */
196 /* accessed sequentially with one pass through the band part of A. */
197 /* First form y := beta*y. */
198 if (*beta != 1.f)
199 {
200 if (*incy == 1)
201 {
202 if (*beta == 0.f)
203 {
204 i__1 = leny;
205 for (i__ = 1;
206 i__ <= i__1;
207 ++i__)
208 {
209 y[i__] = 0.f;
210 /* L10: */
211 }
212 }
213 else
214 {
215 i__1 = leny;
216 for (i__ = 1;
217 i__ <= i__1;
218 ++i__)
219 {
220 y[i__] = *beta * y[i__];
221 /* L20: */
222 }
223 }
224 }
225 else
226 {
227 iy = ky;
228 if (*beta == 0.f)
229 {
230 i__1 = leny;
231 for (i__ = 1;
232 i__ <= i__1;
233 ++i__)
234 {
235 y[iy] = 0.f;
236 iy += *incy;
237 /* L30: */
238 }
239 }
240 else
241 {
242 i__1 = leny;
243 for (i__ = 1;
244 i__ <= i__1;
245 ++i__)
246 {
247 y[iy] = *beta * y[iy];
248 iy += *incy;
249 /* L40: */
250 }
251 }
252 }
253 }
254 if (*alpha == 0.f)
255 {
256 return 0;
257 }
258 kup1 = *ku + 1;
259 if (lsame_(trans, "N"))
260 {
261 /* Form y := alpha*A*x + y. */
262 jx = kx;
263 if (*incy == 1)
264 {
265 i__1 = *n;
266 for (j = 1;
267 j <= i__1;
268 ++j)
269 {
270 if (x[jx] != 0.f)
271 {
272 temp = *alpha * x[jx];
273 k = kup1 - j;
274 /* Computing MAX */
275 i__2 = 1, i__3 = j - *ku;
276 /* Computing MIN */
277 i__5 = *m, i__6 = j + *kl;
278 i__4 = min(i__5,i__6);
279 for (i__ = max(i__2,i__3);
280 i__ <= i__4;
281 ++i__)
282 {
283 y[i__] += temp * a[k + i__ + j * a_dim1];
284 /* L50: */
285 }
286 }
287 jx += *incx;
288 /* L60: */
289 }
290 }
291 else
292 {
293 i__1 = *n;
294 for (j = 1;
295 j <= i__1;
296 ++j)
297 {
298 if (x[jx] != 0.f)
299 {
300 temp = *alpha * x[jx];
301 iy = ky;
302 k = kup1 - j;
303 /* Computing MAX */
304 i__4 = 1, i__2 = j - *ku;
305 /* Computing MIN */
306 i__5 = *m, i__6 = j + *kl;
307 i__3 = min(i__5,i__6);
308 for (i__ = max(i__4,i__2);
309 i__ <= i__3;
310 ++i__)
311 {
312 y[iy] += temp * a[k + i__ + j * a_dim1];
313 iy += *incy;
314 /* L70: */
315 }
316 }
317 jx += *incx;
318 if (j > *ku)
319 {
320 ky += *incy;
321 }
322 /* L80: */
323 }
324 }
325 }
326 else
327 {
328 /* Form y := alpha*A'*x + y. */
329 jy = ky;
330 if (*incx == 1)
331 {
332 i__1 = *n;
333 for (j = 1;
334 j <= i__1;
335 ++j)
336 {
337 temp = 0.f;
338 k = kup1 - j;
339 /* Computing MAX */
340 i__3 = 1, i__4 = j - *ku;
341 /* Computing MIN */
342 i__5 = *m, i__6 = j + *kl;
343 i__2 = min(i__5,i__6);
344 for (i__ = max(i__3,i__4);
345 i__ <= i__2;
346 ++i__)
347 {
348 temp += a[k + i__ + j * a_dim1] * x[i__];
349 /* L90: */
350 }
351 y[jy] += *alpha * temp;
352 jy += *incy;
353 /* L100: */
354 }
355 }
356 else
357 {
358 i__1 = *n;
359 for (j = 1;
360 j <= i__1;
361 ++j)
362 {
363 temp = 0.f;
364 ix = kx;
365 k = kup1 - j;
366 /* Computing MAX */
367 i__2 = 1, i__3 = j - *ku;
368 /* Computing MIN */
369 i__5 = *m, i__6 = j + *kl;
370 i__4 = min(i__5,i__6);
371 for (i__ = max(i__2,i__3);
372 i__ <= i__4;
373 ++i__)
374 {
375 temp += a[k + i__ + j * a_dim1] * x[ix];
376 ix += *incx;
377 /* L110: */
378 }
379 y[jy] += *alpha * temp;
380 jy += *incy;
381 if (j > *ku)
382 {
383 kx += *incx;
384 }
385 /* L120: */
386 }
387 }
388 }
389 return 0;
390 /* End of SGBMV . */
391 }
392 /* sgbmv_ */
393
394