dspmv.c
Go to the documentation of this file.
1 /* dspmv.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 dspmv_(char *uplo, integer *n, doublereal *alpha,
16  doublereal *ap, doublereal *x, integer *incx, doublereal *beta,
17  doublereal *y, integer *incy, ftnlen uplo_len)
18 {
19  /* System generated locals */
20  integer i__1, i__2;
21 
22  /* Local variables */
23  integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
24  doublereal temp1, temp2;
25  extern logical lsame_(char *, char *, ftnlen, ftnlen);
26  extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
27 
28 /* .. Scalar Arguments .. */
29 /* .. */
30 /* .. Array Arguments .. */
31 /* .. */
32 
33 /* Purpose */
34 /* ======= */
35 
36 /* DSPMV performs the matrix-vector operation */
37 
38 /* y := alpha*A*x + beta*y, */
39 
40 /* where alpha and beta are scalars, x and y are n element vectors and */
41 /* A is an n by n symmetric matrix, supplied in packed form. */
42 
43 /* Arguments */
44 /* ========== */
45 
46 /* UPLO - CHARACTER*1. */
47 /* On entry, UPLO specifies whether the upper or lower */
48 /* triangular part of the matrix A is supplied in the packed */
49 /* array AP as follows: */
50 
51 /* UPLO = 'U' or 'u' The upper triangular part of A is */
52 /* supplied in AP. */
53 
54 /* UPLO = 'L' or 'l' The lower triangular part of A is */
55 /* supplied in AP. */
56 
57 /* Unchanged on exit. */
58 
59 /* N - INTEGER. */
60 /* On entry, N specifies the order of the matrix A. */
61 /* N must be at least zero. */
62 /* Unchanged on exit. */
63 
64 /* ALPHA - DOUBLE PRECISION. */
65 /* On entry, ALPHA specifies the scalar alpha. */
66 /* Unchanged on exit. */
67 
68 /* AP - DOUBLE PRECISION array of DIMENSION at least */
69 /* ( ( n*( n + 1 ) )/2 ). */
70 /* Before entry with UPLO = 'U' or 'u', the array AP must */
71 /* contain the upper triangular part of the symmetric matrix */
72 /* packed sequentially, column by column, so that AP( 1 ) */
73 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
74 /* and a( 2, 2 ) respectively, and so on. */
75 /* Before entry with UPLO = 'L' or 'l', the array AP must */
76 /* contain the lower triangular part of the symmetric matrix */
77 /* packed sequentially, column by column, so that AP( 1 ) */
78 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
79 /* and a( 3, 1 ) respectively, and so on. */
80 /* Unchanged on exit. */
81 
82 /* X - DOUBLE PRECISION array of dimension at least */
83 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
84 /* Before entry, the incremented array X must contain the n */
85 /* element vector x. */
86 /* Unchanged on exit. */
87 
88 /* INCX - INTEGER. */
89 /* On entry, INCX specifies the increment for the elements of */
90 /* X. INCX must not be zero. */
91 /* Unchanged on exit. */
92 
93 /* BETA - DOUBLE PRECISION. */
94 /* On entry, BETA specifies the scalar beta. When BETA is */
95 /* supplied as zero then Y need not be set on input. */
96 /* Unchanged on exit. */
97 
98 /* Y - DOUBLE PRECISION array of dimension at least */
99 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
100 /* Before entry, the incremented array Y must contain the n */
101 /* element vector y. On exit, Y is overwritten by the updated */
102 /* vector y. */
103 
104 /* INCY - INTEGER. */
105 /* On entry, INCY specifies the increment for the elements of */
106 /* Y. INCY must not be zero. */
107 /* Unchanged on exit. */
108 
109 /* Further Details */
110 /* =============== */
111 
112 /* Level 2 Blas routine. */
113 
114 /* -- Written on 22-October-1986. */
115 /* Jack Dongarra, Argonne National Lab. */
116 /* Jeremy Du Croz, Nag Central Office. */
117 /* Sven Hammarling, Nag Central Office. */
118 /* Richard Hanson, Sandia National Labs. */
119 
120 /* ===================================================================== */
121 
122 /* .. Parameters .. */
123 /* .. */
124 /* .. Local Scalars .. */
125 /* .. */
126 /* .. External Functions .. */
127 /* .. */
128 /* .. External Subroutines .. */
129 /* .. */
130 
131 /* Test the input parameters. */
132 
133  /* Parameter adjustments */
134  --y;
135  --x;
136  --ap;
137 
138  /* Function Body */
139  info = 0;
140  if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
141  ftnlen)1, (ftnlen)1)) {
142  info = 1;
143  } else if (*n < 0) {
144  info = 2;
145  } else if (*incx == 0) {
146  info = 6;
147  } else if (*incy == 0) {
148  info = 9;
149  }
150  if (info != 0) {
151  xerbla_("DSPMV ", &info, (ftnlen)6);
152  return 0;
153  }
154 
155 /* Quick return if possible. */
156 
157  if (*n == 0 || (*alpha == 0. && *beta == 1.)) {
158  return 0;
159  }
160 
161 /* Set up the start points in X and Y. */
162 
163  if (*incx > 0) {
164  kx = 1;
165  } else {
166  kx = 1 - (*n - 1) * *incx;
167  }
168  if (*incy > 0) {
169  ky = 1;
170  } else {
171  ky = 1 - (*n - 1) * *incy;
172  }
173 
174 /* Start the operations. In this version the elements of the array AP */
175 /* are accessed sequentially with one pass through AP. */
176 
177 /* First form y := beta*y. */
178 
179  if (*beta != 1.) {
180  if (*incy == 1) {
181  if (*beta == 0.) {
182  i__1 = *n;
183  for (i__ = 1; i__ <= i__1; ++i__) {
184  y[i__] = 0.;
185 /* L10: */
186  }
187  } else {
188  i__1 = *n;
189  for (i__ = 1; i__ <= i__1; ++i__) {
190  y[i__] = *beta * y[i__];
191 /* L20: */
192  }
193  }
194  } else {
195  iy = ky;
196  if (*beta == 0.) {
197  i__1 = *n;
198  for (i__ = 1; i__ <= i__1; ++i__) {
199  y[iy] = 0.;
200  iy += *incy;
201 /* L30: */
202  }
203  } else {
204  i__1 = *n;
205  for (i__ = 1; i__ <= i__1; ++i__) {
206  y[iy] = *beta * y[iy];
207  iy += *incy;
208 /* L40: */
209  }
210  }
211  }
212  }
213  if (*alpha == 0.) {
214  return 0;
215  }
216  kk = 1;
217  if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
218 
219 /* Form y when AP contains the upper triangle. */
220 
221  if (*incx == 1 && *incy == 1) {
222  i__1 = *n;
223  for (j = 1; j <= i__1; ++j) {
224  temp1 = *alpha * x[j];
225  temp2 = 0.;
226  k = kk;
227  i__2 = j - 1;
228  for (i__ = 1; i__ <= i__2; ++i__) {
229  y[i__] += temp1 * ap[k];
230  temp2 += ap[k] * x[i__];
231  ++k;
232 /* L50: */
233  }
234  y[j] = y[j] + temp1 * ap[kk + j - 1] + *alpha * temp2;
235  kk += j;
236 /* L60: */
237  }
238  } else {
239  jx = kx;
240  jy = ky;
241  i__1 = *n;
242  for (j = 1; j <= i__1; ++j) {
243  temp1 = *alpha * x[jx];
244  temp2 = 0.;
245  ix = kx;
246  iy = ky;
247  i__2 = kk + j - 2;
248  for (k = kk; k <= i__2; ++k) {
249  y[iy] += temp1 * ap[k];
250  temp2 += ap[k] * x[ix];
251  ix += *incx;
252  iy += *incy;
253 /* L70: */
254  }
255  y[jy] = y[jy] + temp1 * ap[kk + j - 1] + *alpha * temp2;
256  jx += *incx;
257  jy += *incy;
258  kk += j;
259 /* L80: */
260  }
261  }
262  } else {
263 
264 /* Form y when AP contains the lower triangle. */
265 
266  if (*incx == 1 && *incy == 1) {
267  i__1 = *n;
268  for (j = 1; j <= i__1; ++j) {
269  temp1 = *alpha * x[j];
270  temp2 = 0.;
271  y[j] += temp1 * ap[kk];
272  k = kk + 1;
273  i__2 = *n;
274  for (i__ = j + 1; i__ <= i__2; ++i__) {
275  y[i__] += temp1 * ap[k];
276  temp2 += ap[k] * x[i__];
277  ++k;
278 /* L90: */
279  }
280  y[j] += *alpha * temp2;
281  kk += *n - j + 1;
282 /* L100: */
283  }
284  } else {
285  jx = kx;
286  jy = ky;
287  i__1 = *n;
288  for (j = 1; j <= i__1; ++j) {
289  temp1 = *alpha * x[jx];
290  temp2 = 0.;
291  y[jy] += temp1 * ap[kk];
292  ix = jx;
293  iy = jy;
294  i__2 = kk + *n - j;
295  for (k = kk + 1; k <= i__2; ++k) {
296  ix += *incx;
297  iy += *incy;
298  y[iy] += temp1 * ap[k];
299  temp2 += ap[k] * x[ix];
300 /* L110: */
301  }
302  y[jy] += *alpha * temp2;
303  jx += *incx;
304  jy += *incy;
305  kk += *n - j + 1;
306 /* L120: */
307  }
308  }
309  }
310 
311  return 0;
312 
313 /* End of DSPMV . */
314 
315 } /* dspmv_ */
316 
dspmv_
int dspmv_(char *uplo, integer *n, doublereal *alpha, doublereal *ap, doublereal *x, integer *incx, doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len)
Definition: dspmv.c:15
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
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
doublereal
double doublereal
Definition: datatypes.h:11
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
y
Scalar * y
Definition: level1_cplx_impl.h:124
integer
int integer
Definition: datatypes.h:8
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
logical
int logical
Definition: datatypes.h:15

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