GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
11void 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}
struct complex Cpx
Definition ccmath.h:41
void qrecvc(double *ev, Cpx *evec, double *dp, int n)
Definition qrecvc.c:11
double re
Definition ccmath.h:39
double im
Definition ccmath.h:39
#define x