Logo ROOT  
Reference Guide
mndspr.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 /* dspr.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 mndspr(const char *uplo, unsigned int n, double alpha, const double *x, int incx, double *ap)
23 {
24  /* System generated locals */
25  int i__1, i__2;
26 
27  /* Local variables */
28  int info;
29  double temp;
30  int i__, j, k;
31  int kk, ix, jx, kx = 0;
32 
33  /* .. Scalar Arguments .. */
34  /* .. Array Arguments .. */
35  /* .. */
36 
37  /* Purpose */
38  /* ======= */
39 
40  /* DSPR performs the symmetric rank 1 operation */
41 
42  /* A := alpha*x*x' + A, */
43 
44  /* where alpha is a real scalar, x is an n element vector and A is an */
45  /* n by n symmetric matrix, supplied in packed form. */
46 
47  /* Parameters */
48  /* ========== */
49 
50  /* UPLO - CHARACTER*1. */
51  /* On entry, UPLO specifies whether the Upper or Lower */
52  /* triangular part of the matrix A is supplied in the packed */
53  /* array AP as follows: */
54 
55  /* UPLO = 'U' or 'u' The Upper triangular part of A is */
56  /* supplied in AP. */
57 
58  /* UPLO = 'L' or 'l' The Lower triangular part of A is */
59  /* supplied in AP. */
60 
61  /* Unchanged on exit. */
62 
63  /* N - INTEGER. */
64  /* On entry, N specifies the order of the matrix A. */
65  /* N must be at least zero. */
66  /* Unchanged on exit. */
67 
68  /* ALPHA - DOUBLE PRECISION. */
69  /* On entry, ALPHA specifies the scalar alpha. */
70  /* Unchanged on exit. */
71 
72  /* X - DOUBLE PRECISION array of dimension at least */
73  /* ( 1 + ( n - 1 )*abs( INCX ) ). */
74  /* Before entry, the incremented array X must contain the n */
75  /* element vector x. */
76  /* Unchanged on exit. */
77 
78  /* INCX - INTEGER. */
79  /* On entry, INCX specifies the increment for the Elements of */
80  /* X. INCX must not be zero. */
81  /* Unchanged on exit. */
82 
83  /* AP - DOUBLE PRECISION array of DIMENSION at least */
84  /* ( ( n*( n + 1 ) )/2 ). */
85  /* Before entry with UPLO = 'U' or 'u', the array AP must */
86  /* contain the Upper triangular part of the symmetric matrix */
87  /* packed sequentially, column by column, so that AP( 1 ) */
88  /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
89  /* and a( 2, 2 ) respectively, and so on. On exit, the array */
90  /* AP is overwritten by the Upper triangular part of the */
91  /* updated matrix. */
92  /* Before entry with UPLO = 'L' or 'l', the array AP must */
93  /* contain the Lower triangular part of the symmetric matrix */
94  /* packed sequentially, column by column, so that AP( 1 ) */
95  /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
96  /* and a( 3, 1 ) respectively, and so on. On exit, the array */
97  /* AP is overwritten by the Lower triangular part of the */
98  /* updated matrix. */
99 
100  /* Level 2 Blas routine. */
101 
102  /* -- Written on 22-October-1986. */
103  /* Jack Dongarra, Argonne National Lab. */
104  /* Jeremy Du Croz, Nag Central Office. */
105  /* Sven Hammarling, Nag Central Office. */
106  /* Richard Hanson, Sandia National Labs. */
107 
108  /* .. Parameters .. */
109  /* .. Local Scalars .. */
110  /* .. External Functions .. */
111  /* .. External Subroutines .. */
112  /* .. */
113  /* .. Executable Statements .. */
114 
115  /* Test the input parameters. */
116 
117  /* Parameter adjustments */
118  --ap;
119  --x;
120 
121  /* Function Body */
122  info = 0;
123  if (!mnlsame(uplo, "U") && !mnlsame(uplo, "L")) {
124  info = 1;
125  }
126  // else if (n < 0) {
127  // info = 2;
128  // }
129  else if (incx == 0) {
130  info = 5;
131  }
132  if (info != 0) {
133  mnxerbla("DSPR ", info);
134  return 0;
135  }
136 
137  /* Quick return if possible. */
138 
139  if (n == 0 || alpha == 0.) {
140  return 0;
141  }
142 
143  /* Set the start point in X if the increment is not unity. */
144 
145  if (incx <= 0) {
146  kx = 1 - (n - 1) * incx;
147  } else if (incx != 1) {
148  kx = 1;
149  }
150 
151  /* Start the operations. In this version the Elements of the array AP */
152  /* are accessed sequentially with one pass through AP. */
153 
154  kk = 1;
155  if (mnlsame(uplo, "U")) {
156 
157  /* Form A when Upper triangle is stored in AP. */
158 
159  if (incx == 1) {
160  i__1 = n;
161  for (j = 1; j <= i__1; ++j) {
162  if (x[j] != 0.) {
163  temp = alpha * x[j];
164  k = kk;
165  i__2 = j;
166  for (i__ = 1; i__ <= i__2; ++i__) {
167  ap[k] += x[i__] * temp;
168  ++k;
169  /* L10: */
170  }
171  }
172  kk += j;
173  /* L20: */
174  }
175  } else {
176  jx = kx;
177  i__1 = n;
178  for (j = 1; j <= i__1; ++j) {
179  if (x[jx] != 0.) {
180  temp = alpha * x[jx];
181  ix = kx;
182  i__2 = kk + j - 1;
183  for (k = kk; k <= i__2; ++k) {
184  ap[k] += x[ix] * temp;
185  ix += incx;
186  /* L30: */
187  }
188  }
189  jx += incx;
190  kk += j;
191  /* L40: */
192  }
193  }
194  } else {
195 
196  /* Form A when Lower triangle is stored in AP. */
197 
198  if (incx == 1) {
199  i__1 = n;
200  for (j = 1; j <= i__1; ++j) {
201  if (x[j] != 0.) {
202  temp = alpha * x[j];
203  k = kk;
204  i__2 = n;
205  for (i__ = j; i__ <= i__2; ++i__) {
206  ap[k] += x[i__] * temp;
207  ++k;
208  /* L50: */
209  }
210  }
211  kk = kk + n - j + 1;
212  /* L60: */
213  }
214  } else {
215  jx = kx;
216  i__1 = n;
217  for (j = 1; j <= i__1; ++j) {
218  if (x[jx] != 0.) {
219  temp = alpha * x[jx];
220  ix = jx;
221  i__2 = kk + n - j;
222  for (k = kk; k <= i__2; ++k) {
223  ap[k] += x[ix] * temp;
224  ix += incx;
225  /* L70: */
226  }
227  }
228  jx += incx;
229  kk = kk + n - j + 1;
230  /* L80: */
231  }
232  }
233  }
234 
235  return 0;
236 
237  /* End of DSPR . */
238 
239 } /* dspr_ */
240 
241 } // namespace Minuit2
242 
243 } // 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::Minuit2::mndspr
int mndspr(const char *, unsigned int, double, const double *, int, double *)
Definition: mndspr.cxx:22
ROOT::Minuit2::mnlsame
bool mnlsame(const char *, const char *)
Definition: mnlsame.cxx:21
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4