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