ssbmv.c
Go to the documentation of this file.
1 /* ssbmv.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 ssbmv_(char *uplo, integer *n, integer *k, real *alpha,
16  real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
17  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  real 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 /* SSBMV 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 - REAL . */
71 /* On entry, ALPHA specifies the scalar alpha. */
72 /* Unchanged on exit. */
73 
74 /* A - REAL 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 - REAL 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 - REAL . */
131 /* On entry, BETA specifies the scalar beta. */
132 /* Unchanged on exit. */
133 
134 /* Y - REAL 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 /* Further Details */
145 /* =============== */
146 
147 /* Level 2 Blas routine. */
148 
149 /* -- Written on 22-October-1986. */
150 /* Jack Dongarra, Argonne National Lab. */
151 /* Jeremy Du Croz, Nag Central Office. */
152 /* Sven Hammarling, Nag Central Office. */
153 /* Richard Hanson, Sandia National Labs. */
154 
155 /* ===================================================================== */
156 
157 /* .. Parameters .. */
158 /* .. */
159 /* .. Local Scalars .. */
160 /* .. */
161 /* .. External Functions .. */
162 /* .. */
163 /* .. External Subroutines .. */
164 /* .. */
165 /* .. Intrinsic Functions .. */
166 /* .. */
167 
168 /* Test the input parameters. */
169 
170  /* Parameter adjustments */
171  a_dim1 = *lda;
172  a_offset = 1 + a_dim1;
173  a -= a_offset;
174  --x;
175  --y;
176 
177  /* Function Body */
178  info = 0;
179  if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
180  ftnlen)1, (ftnlen)1)) {
181  info = 1;
182  } else if (*n < 0) {
183  info = 2;
184  } else if (*k < 0) {
185  info = 3;
186  } else if (*lda < *k + 1) {
187  info = 6;
188  } else if (*incx == 0) {
189  info = 8;
190  } else if (*incy == 0) {
191  info = 11;
192  }
193  if (info != 0) {
194  xerbla_("SSBMV ", &info, (ftnlen)6);
195  return 0;
196  }
197 
198 /* Quick return if possible. */
199 
200  if (*n == 0 || (*alpha == 0.f && *beta == 1.f)) {
201  return 0;
202  }
203 
204 /* Set up the start points in X and Y. */
205 
206  if (*incx > 0) {
207  kx = 1;
208  } else {
209  kx = 1 - (*n - 1) * *incx;
210  }
211  if (*incy > 0) {
212  ky = 1;
213  } else {
214  ky = 1 - (*n - 1) * *incy;
215  }
216 
217 /* Start the operations. In this version the elements of the array A */
218 /* are accessed sequentially with one pass through A. */
219 
220 /* First form y := beta*y. */
221 
222  if (*beta != 1.f) {
223  if (*incy == 1) {
224  if (*beta == 0.f) {
225  i__1 = *n;
226  for (i__ = 1; i__ <= i__1; ++i__) {
227  y[i__] = 0.f;
228 /* L10: */
229  }
230  } else {
231  i__1 = *n;
232  for (i__ = 1; i__ <= i__1; ++i__) {
233  y[i__] = *beta * y[i__];
234 /* L20: */
235  }
236  }
237  } else {
238  iy = ky;
239  if (*beta == 0.f) {
240  i__1 = *n;
241  for (i__ = 1; i__ <= i__1; ++i__) {
242  y[iy] = 0.f;
243  iy += *incy;
244 /* L30: */
245  }
246  } else {
247  i__1 = *n;
248  for (i__ = 1; i__ <= i__1; ++i__) {
249  y[iy] = *beta * y[iy];
250  iy += *incy;
251 /* L40: */
252  }
253  }
254  }
255  }
256  if (*alpha == 0.f) {
257  return 0;
258  }
259  if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
260 
261 /* Form y when upper triangle of A is stored. */
262 
263  kplus1 = *k + 1;
264  if (*incx == 1 && *incy == 1) {
265  i__1 = *n;
266  for (j = 1; j <= i__1; ++j) {
267  temp1 = *alpha * x[j];
268  temp2 = 0.f;
269  l = kplus1 - j;
270 /* Computing MAX */
271  i__2 = 1, i__3 = j - *k;
272  i__4 = j - 1;
273  for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
274  y[i__] += temp1 * a[l + i__ + j * a_dim1];
275  temp2 += a[l + i__ + j * a_dim1] * x[i__];
276 /* L50: */
277  }
278  y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
279 /* L60: */
280  }
281  } else {
282  jx = kx;
283  jy = ky;
284  i__1 = *n;
285  for (j = 1; j <= i__1; ++j) {
286  temp1 = *alpha * x[jx];
287  temp2 = 0.f;
288  ix = kx;
289  iy = ky;
290  l = kplus1 - j;
291 /* Computing MAX */
292  i__4 = 1, i__2 = j - *k;
293  i__3 = j - 1;
294  for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
295  y[iy] += temp1 * a[l + i__ + j * a_dim1];
296  temp2 += a[l + i__ + j * a_dim1] * x[ix];
297  ix += *incx;
298  iy += *incy;
299 /* L70: */
300  }
301  y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha *
302  temp2;
303  jx += *incx;
304  jy += *incy;
305  if (j > *k) {
306  kx += *incx;
307  ky += *incy;
308  }
309 /* L80: */
310  }
311  }
312  } else {
313 
314 /* Form y when lower triangle of A is stored. */
315 
316  if (*incx == 1 && *incy == 1) {
317  i__1 = *n;
318  for (j = 1; j <= i__1; ++j) {
319  temp1 = *alpha * x[j];
320  temp2 = 0.f;
321  y[j] += temp1 * a[j * a_dim1 + 1];
322  l = 1 - j;
323 /* Computing MIN */
324  i__4 = *n, i__2 = j + *k;
325  i__3 = min(i__4,i__2);
326  for (i__ = j + 1; i__ <= i__3; ++i__) {
327  y[i__] += temp1 * a[l + i__ + j * a_dim1];
328  temp2 += a[l + i__ + j * a_dim1] * x[i__];
329 /* L90: */
330  }
331  y[j] += *alpha * temp2;
332 /* L100: */
333  }
334  } else {
335  jx = kx;
336  jy = ky;
337  i__1 = *n;
338  for (j = 1; j <= i__1; ++j) {
339  temp1 = *alpha * x[jx];
340  temp2 = 0.f;
341  y[jy] += temp1 * a[j * a_dim1 + 1];
342  l = 1 - j;
343  ix = jx;
344  iy = jy;
345 /* Computing MIN */
346  i__4 = *n, i__2 = j + *k;
347  i__3 = min(i__4,i__2);
348  for (i__ = j + 1; i__ <= i__3; ++i__) {
349  ix += *incx;
350  iy += *incy;
351  y[iy] += temp1 * a[l + i__ + j * a_dim1];
352  temp2 += a[l + i__ + j * a_dim1] * x[ix];
353 /* L110: */
354  }
355  y[jy] += *alpha * temp2;
356  jx += *incx;
357  jy += *incy;
358 /* L120: */
359  }
360  }
361  }
362 
363  return 0;
364 
365 /* End of SSBMV . */
366 
367 } /* ssbmv_ */
368 
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
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
ssbmv_
int ssbmv_(char *uplo, integer *n, integer *k, real *alpha, real *a, integer *lda, real *x, integer *incx, real *beta, real *y, integer *incy, ftnlen uplo_len)
Definition: ssbmv.c:15
ftnlen
int ftnlen
Definition: datatypes.h:14
beta
double beta(double a, double b)
Definition: beta.c:61
datatypes.h
n
int n
Definition: BiCGSTAB_simple.cpp:1
incy
int RealScalar int RealScalar int * incy
Definition: level1_cplx_impl.h:119
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
lsame_
logical lsame_(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: lsame.c:15
l
static const Line3 l(Rot3(), 1, 1)
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
lda
* lda
Definition: eigenvalues.cpp:59
y
Scalar * y
Definition: level1_cplx_impl.h:124
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
integer
int integer
Definition: datatypes.h:8
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
min
#define min(a, b)
Definition: datatypes.h:19
logical
int logical
Definition: datatypes.h:15
max
#define max(a, b)
Definition: datatypes.h:20
real
Definition: main.h:100

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:04:59