GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
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"
9void 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}
void qrecvc(double *ev, Cpx *evec, double *dp, int n)
Definition: qrecvc.c:9
Definition: ccmath.h:38
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38
#define x