dsbmv.c
Go to the documentation of this file.
1 /* dsbmv.f -- translated by f2c (version 20100827).
2  You must link the resulting object file with libf2c:
3  on Microsoft Windows system, link with libf2c.lib;
4  on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5  or, if you install libf2c.a in a standard place, with -lf2c -lm
6  -- in that order, at the end of the command line, as in
7  cc *.o -lf2c -lm
8  Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9 
10  http://www.netlib.org/f2c/libf2c.zip
11 */
12 
13 #include "datatypes.h"
14 
15 /* Subroutine */ int dsbmv_(char *uplo, integer *n, integer *k, doublereal *
16  alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
17  doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len)
18 {
19  /* System generated locals */
20  integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
21 
22  /* Local variables */
23  integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
24  doublereal temp1, temp2;
25  extern logical lsame_(char *, char *, ftnlen, ftnlen);
26  integer kplus1;
27  extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
28 
29 /* .. Scalar Arguments .. */
30 /* .. */
31 /* .. Array Arguments .. */
32 /* .. */
33 
34 /* Purpose */
35 /* ======= */
36 
37 /* DSBMV performs the matrix-vector operation */
38 
39 /* y := alpha*A*x + beta*y, */
40 
41 /* where alpha and beta are scalars, x and y are n element vectors and */
42 /* A is an n by n symmetric band matrix, with k super-diagonals. */
43 
44 /* Arguments */
45 /* ========== */
46 
47 /* UPLO - CHARACTER*1. */
48 /* On entry, UPLO specifies whether the upper or lower */
49 /* triangular part of the band matrix A is being supplied as */
50 /* follows: */
51 
52 /* UPLO = 'U' or 'u' The upper triangular part of A is */
53 /* being supplied. */
54 
55 /* UPLO = 'L' or 'l' The lower triangular part of A is */
56 /* being supplied. */
57 
58 /* Unchanged on exit. */
59 
60 /* N - INTEGER. */
61 /* On entry, N specifies the order of the matrix A. */
62 /* N must be at least zero. */
63 /* Unchanged on exit. */
64 
65 /* K - INTEGER. */
66 /* On entry, K specifies the number of super-diagonals of the */
67 /* matrix A. K must satisfy 0 .le. K. */
68 /* Unchanged on exit. */
69 
70 /* ALPHA - DOUBLE PRECISION. */
71 /* On entry, ALPHA specifies the scalar alpha. */
72 /* Unchanged on exit. */
73 
74 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
75 /* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
76 /* by n part of the array A must contain the upper triangular */
77 /* band part of the symmetric matrix, supplied column by */
78 /* column, with the leading diagonal of the matrix in row */
79 /* ( k + 1 ) of the array, the first super-diagonal starting at */
80 /* position 2 in row k, and so on. The top left k by k triangle */
81 /* of the array A is not referenced. */
82 /* The following program segment will transfer the upper */
83 /* triangular part of a symmetric band matrix from conventional */
84 /* full matrix storage to band storage: */
85 
86 /* DO 20, J = 1, N */
87 /* M = K + 1 - J */
88 /* DO 10, I = MAX( 1, J - K ), J */
89 /* A( M + I, J ) = matrix( I, J ) */
90 /* 10 CONTINUE */
91 /* 20 CONTINUE */
92 
93 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
94 /* by n part of the array A must contain the lower triangular */
95 /* band part of the symmetric matrix, supplied column by */
96 /* column, with the leading diagonal of the matrix in row 1 of */
97 /* the array, the first sub-diagonal starting at position 1 in */
98 /* row 2, and so on. The bottom right k by k triangle of the */
99 /* array A is not referenced. */
100 /* The following program segment will transfer the lower */
101 /* triangular part of a symmetric band matrix from conventional */
102 /* full matrix storage to band storage: */
103 
104 /* DO 20, J = 1, N */
105 /* M = 1 - J */
106 /* DO 10, I = J, MIN( N, J + K ) */
107 /* A( M + I, J ) = matrix( I, J ) */
108 /* 10 CONTINUE */
109 /* 20 CONTINUE */
110 
111 /* Unchanged on exit. */
112 
113 /* LDA - INTEGER. */
114 /* On entry, LDA specifies the first dimension of A as declared */
115 /* in the calling (sub) program. LDA must be at least */
116 /* ( k + 1 ). */
117 /* Unchanged on exit. */
118 
119 /* X - DOUBLE PRECISION array of DIMENSION at least */
120 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
121 /* Before entry, the incremented array X must contain the */
122 /* vector x. */
123 /* Unchanged on exit. */
124 
125 /* INCX - INTEGER. */
126 /* On entry, INCX specifies the increment for the elements of */
127 /* X. INCX must not be zero. */
128 /* Unchanged on exit. */
129 
130 /* BETA - DOUBLE PRECISION. */
131 /* On entry, BETA specifies the scalar beta. */
132 /* Unchanged on exit. */
133 
134 /* Y - DOUBLE PRECISION array of DIMENSION at least */
135 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
136 /* Before entry, the incremented array Y must contain the */
137 /* vector y. On exit, Y is overwritten by the updated vector y. */
138 
139 /* INCY - INTEGER. */
140 /* On entry, INCY specifies the increment for the elements of */
141 /* Y. INCY must not be zero. */
142 /* Unchanged on exit. */
143 
144 
145 /* Level 2 Blas routine. */
146 
147 /* -- Written on 22-October-1986. */
148 /* Jack Dongarra, Argonne National Lab. */
149 /* Jeremy Du Croz, Nag Central Office. */
150 /* Sven Hammarling, Nag Central Office. */
151 /* Richard Hanson, Sandia National Labs. */
152 
153 /* ===================================================================== */
154 
155 /* .. Parameters .. */
156 /* .. */
157 /* .. Local Scalars .. */
158 /* .. */
159 /* .. External Functions .. */
160 /* .. */
161 /* .. External Subroutines .. */
162 /* .. */
163 /* .. Intrinsic Functions .. */
164 /* .. */
165 
166 /* Test the input parameters. */
167 
168  /* Parameter adjustments */
169  a_dim1 = *lda;
170  a_offset = 1 + a_dim1;
171  a -= a_offset;
172  --x;
173  --y;
174 
175  /* Function Body */
176  info = 0;
177  if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
178  ftnlen)1, (ftnlen)1)) {
179  info = 1;
180  } else if (*n < 0) {
181  info = 2;
182  } else if (*k < 0) {
183  info = 3;
184  } else if (*lda < *k + 1) {
185  info = 6;
186  } else if (*incx == 0) {
187  info = 8;
188  } else if (*incy == 0) {
189  info = 11;
190  }
191  if (info != 0) {
192  xerbla_("DSBMV ", &info, (ftnlen)6);
193  return 0;
194  }
195 
196 /* Quick return if possible. */
197 
198  if (*n == 0 || (*alpha == 0. && *beta == 1.)) {
199  return 0;
200  }
201 
202 /* Set up the start points in X and Y. */
203 
204  if (*incx > 0) {
205  kx = 1;
206  } else {
207  kx = 1 - (*n - 1) * *incx;
208  }
209  if (*incy > 0) {
210  ky = 1;
211  } else {
212  ky = 1 - (*n - 1) * *incy;
213  }
214 
215 /* Start the operations. In this version the elements of the array A */
216 /* are accessed sequentially with one pass through A. */
217 
218 /* First form y := beta*y. */
219 
220  if (*beta != 1.) {
221  if (*incy == 1) {
222  if (*beta == 0.) {
223  i__1 = *n;
224  for (i__ = 1; i__ <= i__1; ++i__) {
225  y[i__] = 0.;
226 /* L10: */
227  }
228  } else {
229  i__1 = *n;
230  for (i__ = 1; i__ <= i__1; ++i__) {
231  y[i__] = *beta * y[i__];
232 /* L20: */
233  }
234  }
235  } else {
236  iy = ky;
237  if (*beta == 0.) {
238  i__1 = *n;
239  for (i__ = 1; i__ <= i__1; ++i__) {
240  y[iy] = 0.;
241  iy += *incy;
242 /* L30: */
243  }
244  } else {
245  i__1 = *n;
246  for (i__ = 1; i__ <= i__1; ++i__) {
247  y[iy] = *beta * y[iy];
248  iy += *incy;
249 /* L40: */
250  }
251  }
252  }
253  }
254  if (*alpha == 0.) {
255  return 0;
256  }
257  if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
258 
259 /* Form y when upper triangle of A is stored. */
260 
261  kplus1 = *k + 1;
262  if (*incx == 1 && *incy == 1) {
263  i__1 = *n;
264  for (j = 1; j <= i__1; ++j) {
265  temp1 = *alpha * x[j];
266  temp2 = 0.;
267  l = kplus1 - j;
268 /* Computing MAX */
269  i__2 = 1, i__3 = j - *k;
270  i__4 = j - 1;
271  for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
272  y[i__] += temp1 * a[l + i__ + j * a_dim1];
273  temp2 += a[l + i__ + j * a_dim1] * x[i__];
274 /* L50: */
275  }
276  y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
277 /* L60: */
278  }
279  } else {
280  jx = kx;
281  jy = ky;
282  i__1 = *n;
283  for (j = 1; j <= i__1; ++j) {
284  temp1 = *alpha * x[jx];
285  temp2 = 0.;
286  ix = kx;
287  iy = ky;
288  l = kplus1 - j;
289 /* Computing MAX */
290  i__4 = 1, i__2 = j - *k;
291  i__3 = j - 1;
292  for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
293  y[iy] += temp1 * a[l + i__ + j * a_dim1];
294  temp2 += a[l + i__ + j * a_dim1] * x[ix];
295  ix += *incx;
296  iy += *incy;
297 /* L70: */
298  }
299  y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha *
300  temp2;
301  jx += *incx;
302  jy += *incy;
303  if (j > *k) {
304  kx += *incx;
305  ky += *incy;
306  }
307 /* L80: */
308  }
309  }
310  } else {
311 
312 /* Form y when lower triangle of A is stored. */
313 
314  if (*incx == 1 && *incy == 1) {
315  i__1 = *n;
316  for (j = 1; j <= i__1; ++j) {
317  temp1 = *alpha * x[j];
318  temp2 = 0.;
319  y[j] += temp1 * a[j * a_dim1 + 1];
320  l = 1 - j;
321 /* Computing MIN */
322  i__4 = *n, i__2 = j + *k;
323  i__3 = min(i__4,i__2);
324  for (i__ = j + 1; i__ <= i__3; ++i__) {
325  y[i__] += temp1 * a[l + i__ + j * a_dim1];
326  temp2 += a[l + i__ + j * a_dim1] * x[i__];
327 /* L90: */
328  }
329  y[j] += *alpha * temp2;
330 /* L100: */
331  }
332  } else {
333  jx = kx;
334  jy = ky;
335  i__1 = *n;
336  for (j = 1; j <= i__1; ++j) {
337  temp1 = *alpha * x[jx];
338  temp2 = 0.;
339  y[jy] += temp1 * a[j * a_dim1 + 1];
340  l = 1 - j;
341  ix = jx;
342  iy = jy;
343 /* Computing MIN */
344  i__4 = *n, i__2 = j + *k;
345  i__3 = min(i__4,i__2);
346  for (i__ = j + 1; i__ <= i__3; ++i__) {
347  ix += *incx;
348  iy += *incy;
349  y[iy] += temp1 * a[l + i__ + j * a_dim1];
350  temp2 += a[l + i__ + j * a_dim1] * x[ix];
351 /* L110: */
352  }
353  y[jy] += *alpha * temp2;
354  jx += *incx;
355  jy += *incy;
356 /* L120: */
357  }
358  }
359  }
360 
361  return 0;
362 
363 /* End of DSBMV . */
364 
365 } /* dsbmv_ */
366 
max
#define max(a, b)
Definition: datatypes.h:20
y
Scalar * y
Definition: level1_cplx_impl.h:124
min
#define min(a, b)
Definition: datatypes.h:19
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
n
int n
Definition: BiCGSTAB_simple.cpp:1
ftnlen
int ftnlen
Definition: datatypes.h:14
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
l
static const Line3 l(Rot3(), 1, 1)
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
lda
* lda
Definition: eigenvalues.cpp:59
logical
int logical
Definition: datatypes.h:15
lsame_
logical lsame_(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: lsame.c:15
doublereal
double doublereal
Definition: datatypes.h:11
dsbmv_
int dsbmv_(char *uplo, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len)
Definition: dsbmv.c:15
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
integer
int integer
Definition: datatypes.h:8
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
incy
int RealScalar int RealScalar int * incy
Definition: level1_cplx_impl.h:119

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:11