GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
qrevec.c
Go to the documentation of this file.
1/* qrevec.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 <math.h>
9int qrevec(double *ev, double *evec, double *dp, int n)
10{
11 double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
12
13 int i, j, k, m, mqr = 8 * n;
14
15 double *p;
16
17 for (j = 0, m = n - 1;; ++j) {
18 while (1) {
19 if (m < 1)
20 return 0;
21 k = m - 1;
22 if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
23 --m;
24 else {
25 x = (ev[k] - ev[m]) / 2.;
26 h = sqrt(x * x + dp[k] * dp[k]);
27 if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
28 break;
29 if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
30 sc = dp[k] / (2. * cc * h);
31 else
32 sc = 1.;
33 x += ev[m];
34 ev[m--] = x - h;
35 ev[m--] = x + h;
36 for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
37 h = p[0];
38 p[0] = cc * h + sc * p[n];
39 p[n] = cc * p[n] - sc * h;
40 }
41 }
42 }
43 if (j > mqr)
44 return -1;
45 if (x > 0.)
46 d = ev[m] + x - h;
47 else
48 d = ev[m] + x + h;
49 cc = 1.;
50 y = 0.;
51 ev[0] -= d;
52 for (k = 0; k < m; ++k) {
53 x = ev[k] * cc - y;
54 y = dp[k] * cc;
55 h = sqrt(x * x + dp[k] * dp[k]);
56 if (k > 0)
57 dp[k - 1] = sc * h;
58 ev[k] = cc * h;
59 cc = x / h;
60 sc = dp[k] / h;
61 ev[k + 1] -= d;
62 y *= sc;
63 ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
64 for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
65 h = p[0];
66 p[0] = cc * h + sc * p[n];
67 p[n] = cc * p[n] - sc * h;
68 }
69 }
70 ev[k] = ev[k] * cc - y;
71 dp[k - 1] = ev[k] * sc;
72 ev[k] = ev[k] * cc + d;
73 }
74 return 0;
75}
int qrevec(double *ev, double *evec, double *dp, int n)
Definition qrevec.c:9
#define x