GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
qreval.c
Go to the documentation of this file.
1 /* qreval.c CCMATH mathematics library source code.
2  *
3  * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4  * This code may be redistributed under the terms of the GNU library
5  * public license (LGPL). ( See the lgpl.license file for details.)
6  * ------------------------------------------------------------------------
7  */
8 #include "ccmath.h"
9 int qreval(double *ev, double *dp, int n)
10 {
11  double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
12 
13  int j, k, m, mqr = 8 * n;
14 
15  for (j = 0, m = n - 1;; ++j) {
16  while (1) {
17  if (m < 1)
18  return 0;
19  k = m - 1;
20  if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
21  --m;
22  else {
23  x = (ev[k] - ev[m]) / 2.;
24  h = sqrt(x * x + dp[k] * dp[k]);
25  if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
26  break;
27  x += ev[m];
28  ev[m--] = x - h;
29  ev[m--] = x + h;
30  }
31  }
32  if (j > mqr)
33  return -1;
34  if (x > 0.)
35  d = ev[m] + x - h;
36  else
37  d = ev[m] + x + h;
38  cc = 1.;
39  y = 0.;
40  ev[0] -= d;
41  for (k = 0; k < m; ++k) {
42  x = ev[k] * cc - y;
43  y = dp[k] * cc;
44  h = sqrt(x * x + dp[k] * dp[k]);
45  if (k > 0)
46  dp[k - 1] = sc * h;
47  ev[k] = cc * h;
48  cc = x / h;
49  sc = dp[k] / h;
50  ev[k + 1] -= d;
51  y *= sc;
52  ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
53  }
54  ev[k] = ev[k] * cc - y;
55  dp[k - 1] = ev[k] * sc;
56  ev[k] = ev[k] * cc + d;
57  }
58  return 0;
59 }
#define x
int qreval(double *ev, double *dp, int n)
Definition: qreval.c:9