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
17namespace Minuit2 {
18
19bool mnlsame(const char *, const char *);
20int mnxerbla(const char *, int);
21
22int 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
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:21
int mnxerbla(const char *, int)
Definition: mnxerbla.cxx:26
int Mndspmv(const char *, unsigned int, double, const double *, const double *, int, double, double *, int)
Definition: mndspmv.cxx:22
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.