Logo ROOT  
Reference Guide
mnteigen.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 /* mneig.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 #include <cmath>
16 
17 namespace ROOT {
18 
19 namespace Minuit2 {
20 
21 int mneigen(double *a, unsigned int ndima, unsigned int n, unsigned int mits, double *work, double precis)
22 {
23  // compute matrix eignevalues (transaltion from mneig.F of Minuit)
24 
25  /* System generated locals */
26  unsigned int a_dim1, a_offset, i__1, i__2, i__3;
27  double r__1, r__2;
28 
29  /* Local variables */
30  double b, c__, f, h__;
31  unsigned int i__, j, k, l, m = 0;
32  double r__, s;
33  unsigned int i0, i1, j1, m1, n1;
34  double hh, gl, pr, pt;
35 
36  /* PRECIS is the machine precision EPSMAC */
37  /* Parameter adjustments */
38  a_dim1 = ndima;
39  a_offset = 1 + a_dim1 * 1;
40  a -= a_offset;
41  --work;
42 
43  /* Function Body */
44  int ifault = 1;
45 
46  i__ = n;
47  i__1 = n;
48  for (i1 = 2; i1 <= i__1; ++i1) {
49  l = i__ - 2;
50  f = a[i__ + (i__ - 1) * a_dim1];
51  gl = (double)0.;
52 
53  if (l < 1) {
54  goto L25;
55  }
56 
57  i__2 = l;
58  for (k = 1; k <= i__2; ++k) {
59  /* Computing 2nd power */
60  r__1 = a[i__ + k * a_dim1];
61  gl += r__1 * r__1;
62  }
63  L25:
64  /* Computing 2nd power */
65  r__1 = f;
66  h__ = gl + r__1 * r__1;
67 
68  if (gl > (double)1e-35) {
69  goto L30;
70  }
71 
72  work[i__] = (double)0.;
73  work[n + i__] = f;
74  goto L65;
75  L30:
76  ++l;
77 
78  gl = std::sqrt(h__);
79 
80  if (f >= (double)0.) {
81  gl = -gl;
82  }
83 
84  work[n + i__] = gl;
85  h__ -= f * gl;
86  a[i__ + (i__ - 1) * a_dim1] = f - gl;
87  f = (double)0.;
88  i__2 = l;
89  for (j = 1; j <= i__2; ++j) {
90  a[j + i__ * a_dim1] = a[i__ + j * a_dim1] / h__;
91  gl = (double)0.;
92  i__3 = j;
93  for (k = 1; k <= i__3; ++k) {
94  gl += a[j + k * a_dim1] * a[i__ + k * a_dim1];
95  }
96 
97  if (j >= l) {
98  goto L47;
99  }
100 
101  j1 = j + 1;
102  i__3 = l;
103  for (k = j1; k <= i__3; ++k) {
104  gl += a[k + j * a_dim1] * a[i__ + k * a_dim1];
105  }
106  L47:
107  work[n + j] = gl / h__;
108  f += gl * a[j + i__ * a_dim1];
109  }
110  hh = f / (h__ + h__);
111  i__2 = l;
112  for (j = 1; j <= i__2; ++j) {
113  f = a[i__ + j * a_dim1];
114  gl = work[n + j] - hh * f;
115  work[n + j] = gl;
116  i__3 = j;
117  for (k = 1; k <= i__3; ++k) {
118  a[j + k * a_dim1] = a[j + k * a_dim1] - f * work[n + k] - gl * a[i__ + k * a_dim1];
119  }
120  }
121  work[i__] = h__;
122  L65:
123  --i__;
124  }
125  work[1] = (double)0.;
126  work[n + 1] = (double)0.;
127  i__1 = n;
128  for (i__ = 1; i__ <= i__1; ++i__) {
129  l = i__ - 1;
130 
131  if (work[i__] == (double)0. || l == 0) {
132  goto L100;
133  }
134 
135  i__3 = l;
136  for (j = 1; j <= i__3; ++j) {
137  gl = (double)0.;
138  i__2 = l;
139  for (k = 1; k <= i__2; ++k) {
140  gl += a[i__ + k * a_dim1] * a[k + j * a_dim1];
141  }
142  i__2 = l;
143  for (k = 1; k <= i__2; ++k) {
144  a[k + j * a_dim1] -= gl * a[k + i__ * a_dim1];
145  }
146  }
147  L100:
148  work[i__] = a[i__ + i__ * a_dim1];
149  a[i__ + i__ * a_dim1] = (double)1.;
150 
151  if (l == 0) {
152  goto L110;
153  }
154 
155  i__2 = l;
156  for (j = 1; j <= i__2; ++j) {
157  a[i__ + j * a_dim1] = (double)0.;
158  a[j + i__ * a_dim1] = (double)0.;
159  }
160  L110:;
161  }
162 
163  n1 = n - 1;
164  i__1 = n;
165  for (i__ = 2; i__ <= i__1; ++i__) {
166  i0 = n + i__ - 1;
167  work[i0] = work[i0 + 1];
168  }
169  work[n + n] = (double)0.;
170  b = (double)0.;
171  f = (double)0.;
172  i__1 = n;
173  for (l = 1; l <= i__1; ++l) {
174  j = 0;
175  h__ = precis * ((r__1 = work[l], std::fabs(r__1)) + (r__2 = work[n + l], std::fabs(r__2)));
176 
177  if (b < h__) {
178  b = h__;
179  }
180 
181  i__2 = n;
182  for (m1 = l; m1 <= i__2; ++m1) {
183  m = m1;
184 
185  if ((r__1 = work[n + m], std::fabs(r__1)) <= b) {
186  goto L150;
187  }
188  }
189 
190  L150:
191  if (m == l) {
192  goto L205;
193  }
194 
195  L160:
196  if (j == mits) {
197  return ifault;
198  }
199 
200  ++j;
201  pt = (work[l + 1] - work[l]) / (work[n + l] * (double)2.);
202  r__ = std::sqrt(pt * pt + (double)1.);
203  pr = pt + r__;
204 
205  if (pt < (double)0.) {
206  pr = pt - r__;
207  }
208 
209  h__ = work[l] - work[n + l] / pr;
210  i__2 = n;
211  for (i__ = l; i__ <= i__2; ++i__) {
212  work[i__] -= h__;
213  }
214  f += h__;
215  pt = work[m];
216  c__ = (double)1.;
217  s = (double)0.;
218  m1 = m - 1;
219  i__ = m;
220  i__2 = m1;
221  for (i1 = l; i1 <= i__2; ++i1) {
222  j = i__;
223  --i__;
224  gl = c__ * work[n + i__];
225  h__ = c__ * pt;
226 
227  if (std::fabs(pt) >= (r__1 = work[n + i__], std::fabs(r__1))) {
228  goto L180;
229  }
230 
231  c__ = pt / work[n + i__];
232  r__ = std::sqrt(c__ * c__ + (double)1.);
233  work[n + j] = s * work[n + i__] * r__;
234  s = (double)1. / r__;
235  c__ /= r__;
236  goto L190;
237  L180:
238  c__ = work[n + i__] / pt;
239  r__ = std::sqrt(c__ * c__ + (double)1.);
240  work[n + j] = s * pt * r__;
241  s = c__ / r__;
242  c__ = (double)1. / r__;
243  L190:
244  pt = c__ * work[i__] - s * gl;
245  work[j] = h__ + s * (c__ * gl + s * work[i__]);
246  i__3 = n;
247  for (k = 1; k <= i__3; ++k) {
248  h__ = a[k + j * a_dim1];
249  a[k + j * a_dim1] = s * a[k + i__ * a_dim1] + c__ * h__;
250  a[k + i__ * a_dim1] = c__ * a[k + i__ * a_dim1] - s * h__;
251  }
252  }
253  work[n + l] = s * pt;
254  work[l] = c__ * pt;
255 
256  if ((r__1 = work[n + l], std::fabs(r__1)) > b) {
257  goto L160;
258  }
259 
260  L205:
261  work[l] += f;
262  }
263  i__1 = n1;
264  for (i__ = 1; i__ <= i__1; ++i__) {
265  k = i__;
266  pt = work[i__];
267  i1 = i__ + 1;
268  i__3 = n;
269  for (j = i1; j <= i__3; ++j) {
270 
271  if (work[j] >= pt) {
272  goto L220;
273  }
274 
275  k = j;
276  pt = work[j];
277  L220:;
278  }
279 
280  if (k == i__) {
281  goto L240;
282  }
283 
284  work[k] = work[i__];
285  work[i__] = pt;
286  i__3 = n;
287  for (j = 1; j <= i__3; ++j) {
288  pt = a[j + i__ * a_dim1];
289  a[j + i__ * a_dim1] = a[j + k * a_dim1];
290  a[j + k * a_dim1] = pt;
291  }
292  L240:;
293  }
294  ifault = 0;
295 
296  return ifault;
297 } /* mneig_ */
298 
299 } // namespace Minuit2
300 
301 } // namespace ROOT
l
auto * l
Definition: textangle.C:4
m
auto * m
Definition: textangle.C:8
n
const Int_t n
Definition: legend1.C:16
e
#define e(i)
Definition: RSha256.hxx:103
f
#define f(i)
Definition: RSha256.hxx:104
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
b
#define b(i)
Definition: RSha256.hxx:100
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
a
auto * a
Definition: textangle.C:12
ROOT::Minuit2::mneigen
int mneigen(double *, unsigned int, unsigned int, unsigned int, double *, double)
Definition: mnteigen.cxx:21
double
double
Definition: Converters.cxx:921
sqrt
double sqrt(double)
pt
TPaveText * pt
Definition: entrylist_figure1.C:7
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4