Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
17 namespace ROOT {
18 
19  namespace Minuit2 {
20 
21 
22 int 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  }
65 L25:
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;
77 L30:
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  }
108 L47:
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__;
125 L65:
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  }
150 L100:
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  }
163 L110:
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 
197 L150:
198  if (m == l) {
199  goto L205;
200  }
201 
202 L160:
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;
244 L180:
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__;
250 L190:
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 
267 L205:
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];
284 L220:
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  }
300 L240:
301  ;
302  }
303  ifault = 0;
304 
305  return ifault;
306 } /* mneig_ */
307 
308 
309  } // namespace Minuit2
310 
311 } // namespace ROOT