GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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"
9 int 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 }
#define x
double t
Definition: r_raster.c:39
double b
Definition: r_raster.c:39
int qrbdi(double *dm, double *em, int m)
Definition: qrbdi.c:9