GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 #include "ccmath.h"
9 void qrecvc(double *ev, Cpx * evec, double *dp, int n)
10 {
11  double cc, sc = 0.0, d, x = 0.0, y, h = 0.0, tzr = 1.e-15;
12 
13  int i, j, k, m, nqr = 50 * n;
14 
15  Cpx *p;
16 
17  for (j = 0, m = n - 1; j < nqr; ++j) {
18  while (1) {
19  if (m < 1)
20  break;
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].re;
38  p[0].re = cc * h + sc * p[n].re;
39  p[n].re = cc * p[n].re - sc * h;
40  h = p[0].im;
41  p[0].im = cc * h + sc * p[n].im;
42  p[n].im = cc * p[n].im - sc * h;
43  }
44  }
45  }
46  if (x > 0.)
47  d = ev[m] + x - h;
48  else
49  d = ev[m] + x + h;
50  cc = 1.;
51  y = 0.;
52  ev[0] -= d;
53  for (k = 0; k < m; ++k) {
54  x = ev[k] * cc - y;
55  y = dp[k] * cc;
56  h = sqrt(x * x + dp[k] * dp[k]);
57  if (k > 0)
58  dp[k - 1] = sc * h;
59  ev[k] = cc * h;
60  cc = x / h;
61  sc = dp[k] / h;
62  ev[k + 1] -= d;
63  y *= sc;
64  ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
65  for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
66  h = p[0].re;
67  p[0].re = cc * h + sc * p[n].re;
68  p[n].re = cc * p[n].re - sc * h;
69  h = p[0].im;
70  p[0].im = cc * h + sc * p[n].im;
71  p[n].im = cc * p[n].im - sc * h;
72  }
73  }
74  ev[k] = ev[k] * cc - y;
75  dp[k - 1] = ev[k] * sc;
76  ev[k] = ev[k] * cc + d;
77  }
78 }
#define x
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38
Definition: la.h:54
void qrecvc(double *ev, Cpx *evec, double *dp, int n)
Definition: qrecvc.c:9