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