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
17namespace ROOT {
18
19namespace Minuit2 {
20
21int 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
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
TPaveText * pt
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
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:21
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double s
TMarker m
Definition: textangle.C:8
TLine l
Definition: textangle.C:4
TArc a
Definition: textangle.C:12