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