GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
qrbdu1.c
Go to the documentation of this file.
1/* qrbdu1.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 qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)
10{
11 int i, j, k, n, jj, nm;
12
13 double u, x, y, a, b, c, s, t, w, *p, *q;
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 nm = m;
21 for (j = 0; m > 1 && j < n; ++j) {
22 for (k = m - 1; k > 0; --k) {
23 if (fabs(em[k - 1]) < t)
24 break;
25 if (fabs(dm[k - 1]) < t) {
26 for (i = k, s = 1., c = 0.; i < m; ++i) {
27 a = s * em[i - 1];
28 b = dm[i];
29 em[i - 1] *= c;
30 dm[i] = u = sqrt(a * a + b * b);
31 s = -a / u;
32 c = b / u;
33 for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += nm) {
34 q = p + i - k + 1;
35 w = c * *p + s * *q;
36 *q = c * *q - s * *p;
37 *p = w;
38 }
39 }
40 break;
41 }
42 }
43 y = dm[k];
44 x = dm[m - 1];
45 u = em[m - 2];
46 a = (y + x) * (y - x) - u * u;
47 s = y * em[k];
48 b = s + s;
49 u = sqrt(a * a + b * b);
50 if (u > 0.) {
51 c = sqrt((u + a) / (u + u));
52 if (c != 0.)
53 s /= (c * u);
54 else
55 s = 1.;
56 for (i = k; i < m - 1; ++i) {
57 b = em[i];
58 if (i > k) {
59 a = s * em[i];
60 b *= c;
61 em[i - 1] = u = sqrt(x * x + a * a);
62 c = x / u;
63 s = a / u;
64 }
65 a = c * y + s * b;
66 b = c * b - s * y;
67 for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
68 w = c * *p + s * *(p + 1);
69 *(p + 1) = c * *(p + 1) - s * *p;
70 *p = w;
71 }
72 s *= dm[i + 1];
73 dm[i] = u = sqrt(a * a + s * s);
74 y = c * dm[i + 1];
75 c = a / u;
76 s /= u;
77 x = c * b + s * y;
78 y = c * y - s * b;
79 for (jj = 0, p = um + i; jj < mm; ++jj, p += nm) {
80 w = c * *p + s * *(p + 1);
81 *(p + 1) = c * *(p + 1) - s * *p;
82 *p = w;
83 }
84 }
85 }
86 em[m - 2] = x;
87 dm[m - 1] = y;
88 if (fabs(x) < t)
89 --m;
90 if (m == k + 1)
91 --m;
92 }
93 return j;
94}
double b
double t
int qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)
Definition qrbdu1.c:9
#define x