GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-381f661b3f
qrecvc.c
Go to the documentation of this file.
1 /* qrecvc.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
9 #include "ccmath.h"
10
11 void qrecvc(double *ev, Cpx *evec, double *dp, int n)
12 {
13  double cc, sc = 0.0, d, x = 0.0, y, h = 0.0, tzr = 1.e-15;
14
15  int i, j, k, m, nqr = 50 * n;
16
17  Cpx *p;
18
19  for (j = 0, m = n - 1; j < nqr; ++j) {
20  while (1) {
21  if (m < 1)
22  break;
23  k = m - 1;
24  if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
25  --m;
26  else {
27  x = (ev[k] - ev[m]) / 2.;
28  h = sqrt(x * x + dp[k] * dp[k]);
29  if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
30  break;
31  if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
32  sc = dp[k] / (2. * cc * h);
33  else
34  sc = 1.;
35  x += ev[m];
36  ev[m--] = x - h;
37  ev[m--] = x + h;
38  for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
39  h = p[0].re;
40  p[0].re = cc * h + sc * p[n].re;
41  p[n].re = cc * p[n].re - sc * h;
42  h = p[0].im;
43  p[0].im = cc * h + sc * p[n].im;
44  p[n].im = cc * p[n].im - sc * h;
45  }
46  }
47  }
48  if (x > 0.)
49  d = ev[m] + x - h;
50  else
51  d = ev[m] + x + h;
52  cc = 1.;
53  y = 0.;
54  ev[0] -= d;
55  for (k = 0; k < m; ++k) {
56  x = ev[k] * cc - y;
57  y = dp[k] * cc;
58  h = sqrt(x * x + dp[k] * dp[k]);
59  if (k > 0)
60  dp[k - 1] = sc * h;
61  ev[k] = cc * h;
62  cc = x / h;
63  sc = dp[k] / h;
64  ev[k + 1] -= d;
65  y *= sc;
66  ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
67  for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
68  h = p[0].re;
69  p[0].re = cc * h + sc * p[n].re;
70  p[n].re = cc * p[n].re - sc * h;
71  h = p[0].im;
72  p[0].im = cc * h + sc * p[n].im;
73  p[n].im = cc * p[n].im - sc * h;
74  }
75  }
76  ev[k] = ev[k] * cc - y;
77  dp[k - 1] = ev[k] * sc;
78  ev[k] = ev[k] * cc + d;
79  }
80 }
void qrecvc(double *ev, Cpx *evec, double *dp, int n)
Definition: qrecvc.c:11
Definition: la.h:54
double re
Definition: ccmath.h:39
double im
Definition: ccmath.h:39
#define x