1 /* eigen/qrstep.c
2 *
3 * Copyright (C) 2007, 2010 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* remove off-diagonal elements which are neglegible compared with the
21 neighboring diagonal elements */
22
23 static void
chop_small_elements(const size_t N,const double d[],double sd[])24 chop_small_elements (const size_t N, const double d[], double sd[])
25 {
26 double d_i = d[0];
27
28 size_t i;
29
30 for (i = 0; i < N - 1; i++)
31 {
32 double sd_i = sd[i];
33 double d_ip1 = d[i + 1];
34
35 if (fabs (sd_i) < GSL_DBL_EPSILON * (fabs (d_i) + fabs (d_ip1)))
36 {
37 sd[i] = 0.0;
38 }
39 d_i = d_ip1;
40 }
41 }
42
43 /* Generate a Givens rotation (cos,sin) which takes v=(x,y) to (|v|,0)
44
45 From Golub and Van Loan, "Matrix Computations", Section 5.1.8 */
46
47 inline static void
create_givens(const double a,const double b,double * c,double * s)48 create_givens (const double a, const double b, double *c, double *s)
49 {
50 if (b == 0)
51 {
52 *c = 1;
53 *s = 0;
54 }
55 else if (fabs (b) > fabs (a))
56 {
57 double t = -a / b;
58 double s1 = 1.0 / sqrt (1 + t * t);
59 *s = s1;
60 *c = s1 * t;
61 }
62 else
63 {
64 double t = -b / a;
65 double c1 = 1.0 / sqrt (1 + t * t);
66 *c = c1;
67 *s = c1 * t;
68 }
69 }
70
71 inline static double
trailing_eigenvalue(const size_t n,const double d[],const double sd[])72 trailing_eigenvalue (const size_t n, const double d[], const double sd[])
73 {
74 double ta = d[n - 2];
75 double tb = d[n - 1];
76 double tab = sd[n - 2];
77
78 double dt = (ta - tb) / 2.0;
79
80 double mu;
81
82 if (dt > 0)
83 {
84 mu = tb - tab * (tab / (dt + hypot (dt, tab)));
85 }
86 else if (dt == 0)
87 {
88 mu = tb - fabs(tab);
89 }
90 else
91 {
92 mu = tb + tab * (tab / ((-dt) + hypot (dt, tab)));
93 }
94
95 return mu;
96 }
97
98 static void
qrstep(const size_t n,double d[],double sd[],double gc[],double gs[])99 qrstep (const size_t n, double d[], double sd[], double gc[], double gs[])
100 {
101 double x, z;
102 double ak, bk, zk, ap, bp, aq, bq;
103 size_t k;
104
105 double mu = trailing_eigenvalue (n, d, sd);
106
107 /* If mu is large relative to d_0 and sd_0 then the Givens rotation
108 will have no effect, leading to an infinite loop.
109
110 We set mu to zero in this case, which at least diagonalises the
111 submatrix [d_0, sd_0 ; sd_0, d_0] and allows further progress. */
112
113 if (GSL_DBL_EPSILON * fabs(mu) > (fabs(d[0]) + fabs(sd[0]))) {
114 mu = 0;
115 }
116
117 x = d[0] - mu;
118 z = sd[0];
119
120 ak = 0;
121 bk = 0;
122 zk = 0;
123
124 ap = d[0];
125 bp = sd[0];
126
127 aq = d[1];
128
129 if (n == 2)
130 {
131 double c, s;
132 create_givens (x, z, &c, &s);
133
134 if (gc != NULL)
135 gc[0] = c;
136 if (gs != NULL)
137 gs[0] = s;
138
139 {
140 double ap1 = c * (c * ap - s * bp) + s * (s * aq - c * bp);
141 double bp1 = c * (s * ap + c * bp) - s * (s * bp + c * aq);
142
143 double aq1 = s * (s * ap + c * bp) + c * (s * bp + c * aq);
144
145 ak = ap1;
146 bk = bp1;
147
148 ap = aq1;
149 }
150
151 d[0] = ak;
152 sd[0] = bk;
153 d[1] = ap;
154
155 return;
156 }
157
158 bq = sd[1];
159
160 for (k = 0; k < n - 1; k++)
161 {
162 double c, s;
163 create_givens (x, z, &c, &s);
164
165 /* store Givens rotation */
166 if (gc != NULL)
167 gc[k] = c;
168 if (gs != NULL)
169 gs[k] = s;
170
171 /* compute G' T G */
172
173 {
174 double bk1 = c * bk - s * zk;
175
176 double ap1 = c * (c * ap - s * bp) + s * (s * aq - c * bp);
177 double bp1 = c * (s * ap + c * bp) - s * (s * bp + c * aq);
178 double zp1 = -s * bq;
179
180 double aq1 = s * (s * ap + c * bp) + c * (s * bp + c * aq);
181 double bq1 = c * bq;
182
183 ak = ap1;
184 bk = bp1;
185 zk = zp1;
186
187 ap = aq1;
188 bp = bq1;
189
190 if (k < n - 2)
191 aq = d[k + 2];
192 if (k < n - 3)
193 bq = sd[k + 2];
194
195 d[k] = ak;
196
197 if (k > 0)
198 sd[k - 1] = bk1;
199
200 if (k < n - 2)
201 sd[k + 1] = bp;
202
203 x = bk;
204 z = zk;
205 }
206 }
207
208 /* k = n - 1 */
209 d[k] = ap;
210 sd[k - 1] = bk;
211 }
212