GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-7d4cd0c030
qrbdu1.c
Go to the documentation of this file.
1 /* qrbdu1.c CCMATH mathematics library source code.
2  *
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 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 }
int qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)
Definition: qrbdu1.c:9
double b
Definition: r_raster.c:39
double t
Definition: r_raster.c:39
#define x