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