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
15namespace ROOT {
16
17 namespace Minuit2 {
18
19
20bool mnlsame(const char*, const char*);
21int mnxerbla(const char*, int);
22
23int 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
double beta(double x, double y)
Calculates the beta function.
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
bool mnlsame(const char *, const char *)
Definition: mnlsame.cxx:22
int mnxerbla(const char *, int)
Definition: mnxerbla.cxx:27
int Mndspmv(const char *, unsigned int, double, const double *, const double *, int, double, double *, int)
Definition: mndspmv.cxx:23
VSD Structures.
Definition: StringConv.hxx:21