GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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"
9int 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}
int qreval(double *ev, double *dp, int n)
Definition qreval.c:9
#define x