GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71915
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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>
9 int 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 }
#define x
int qrevec(double *eval, double *evec, double *dp, int n)
Definition: qrevec.c:9