ARTS 2.5.4 (git: 31ce4f0e)
poly_roots.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012 Oliver Lemke <olemke@core-dump.info>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
19// File description
21
34// External declarations
36
37#include "poly_roots.h"
38#include <cmath>
39#include <iostream>
40#include <stdexcept>
41
42/* C-style matrix elements */
43#define MAT(m, i, j, n) ((m)[(i) * (n) + (j)])
44
45/* Fortran-style matrix elements */
46#define FMAT(m, i, j, n) ((m)[((i)-1) * (n) + ((j)-1)])
47
48#define GSL_DBL_EPSILON 2.2204460492503131e-16
49
50#define RADIX 2
51#define RADIX2 (RADIX * RADIX)
52
53#define GSL_SUCCESS 0
54#define GSL_FAILURE -1
55#define GSL_EINVAL 4
56#define GSL_EFAILED 5
57
58#define GSL_SET_COMPLEX_PACKED(zp, n, x, y) \
59 do { \
60 *((zp) + 2 * (n)) = (x); \
61 *((zp) + (2 * (n) + 1)) = (y); \
62 } while (0)
63
64typedef double *gsl_complex_packed_ptr;
65
66typedef struct {
67 size_t nc;
68 double *matrix;
70
71/* Begin Internal GSL function prototypes */
72
73static gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n);
74
75static void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w);
76
77static void set_companion_matrix(const double *a, size_t n, double *m);
78
79static int qr_companion(double *h, size_t nc, gsl_complex_packed_ptr z);
80
81static void balance_companion_matrix(double *m, size_t n);
82
83static int gsl_poly_complex_solve(const double *a,
84 size_t n,
87
88/* End Internal GSL function prototypes */
89
90int poly_root_solve(Matrix &roots, Vector &coeffs) {
91 Index a;
92 double *c;
93 double *s;
94
95 a = coeffs.nelem();
96
97 ARTS_ASSERT(roots.nrows() == a - 1);
98 ARTS_ASSERT(roots.ncols() == 2);
99 ARTS_ASSERT(coeffs[a - 1] != 0);
100
101#ifdef USE_DOUBLE
102 c = coeffs.mdata;
103 s = roots.mdata;
104#else
105 c = (double *)malloc(a * sizeof(double));
106 for (Index i = 0; i < a; i++) c[i] = (double)coeffs[i];
107
108 s = (double *)malloc((a - 1) * 2 * sizeof(double));
109#endif
110
111 gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(a);
112
113 int status = gsl_poly_complex_solve(c, a, w, s);
114
115#ifndef USE_DOUBLE
116 if (status == GSL_SUCCESS) {
117 for (Index i = 0; i < a - 1; i++) {
118 roots(i, 0) = (Numeric)s[i * 2];
119 roots(i, 1) = (Numeric)s[i * 2 + 1];
120 }
121 }
122
123 free(c);
124 free(s);
125#endif
126
127 gsl_poly_complex_workspace_free(w);
128
129 return status;
130}
131
132static gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n) {
133 size_t nc;
134
136
137 if (n == 0) {
138 cerr << "matrix size n must be positive integer" << endl;
139
140 return (NULL);
141 }
142
144
145 if (w == 0) {
146 cerr << "failed to allocate space for struct" << endl;
147
148 return (NULL);
149 }
150
151 nc = n - 1;
152
153 w->nc = nc;
154
155 w->matrix = (double *)malloc(nc * nc * sizeof(double));
156
157 if (w->matrix == 0) {
158 free(w); /* error in constructor, avoid memory leak */
159
160 cerr << "failed to allocate space for workspace matrix" << endl;
161
162 return (NULL);
163 }
164
165 return w;
166}
167
168static void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w) {
169 free(w->matrix);
170 free(w);
171}
172
173static void balance_companion_matrix(double *m, size_t nc) {
174 int not_converged = 1;
175
176 double row_norm = 0;
177 double col_norm = 0;
178
179 while (not_converged) {
180 size_t i, j;
181 double g, f, s;
182
183 not_converged = 0;
184
185 for (i = 0; i < nc; i++) {
186 /* column norm, excluding the diagonal */
187
188 if (i != nc - 1) {
189 col_norm = fabs(MAT(m, i + 1, i, nc));
190 } else {
191 col_norm = 0;
192
193 for (j = 0; j < nc - 1; j++) {
194 col_norm += fabs(MAT(m, j, nc - 1, nc));
195 }
196 }
197
198 /* row norm, excluding the diagonal */
199
200 if (i == 0) {
201 row_norm = fabs(MAT(m, 0, nc - 1, nc));
202 } else if (i == nc - 1) {
203 row_norm = fabs(MAT(m, i, i - 1, nc));
204 } else {
205 row_norm = (fabs(MAT(m, i, i - 1, nc)) + fabs(MAT(m, i, nc - 1, nc)));
206 }
207
208 if (col_norm == 0 || row_norm == 0) {
209 continue;
210 }
211
212 g = row_norm / RADIX;
213 f = 1;
214 s = col_norm + row_norm;
215
216 while (col_norm < g) {
217 f *= RADIX;
218 col_norm *= RADIX2;
219 }
220
221 g = row_norm * RADIX;
222
223 while (col_norm > g) {
224 f /= RADIX;
225 col_norm /= RADIX2;
226 }
227
228 if ((row_norm + col_norm) < 0.95 * s * f) {
229 not_converged = 1;
230
231 g = 1 / f;
232
233 if (i == 0) {
234 MAT(m, 0, nc - 1, nc) *= g;
235 } else {
236 MAT(m, i, i - 1, nc) *= g;
237 MAT(m, i, nc - 1, nc) *= g;
238 }
239
240 if (i == nc - 1) {
241 for (j = 0; j < nc; j++) {
242 MAT(m, j, i, nc) *= f;
243 }
244 } else {
245 MAT(m, i + 1, i, nc) *= f;
246 }
247 }
248 }
249 }
250}
251
252static int qr_companion(double *h, size_t nc, gsl_complex_packed_ptr zroot) {
253 double t = 0.0;
254
255 size_t iterations, e, i, j, k, m;
256
257 double w, x, y, s, z;
258
259 double p = 0, q = 0, r = 0;
260
261 /* FIXME: if p,q,r, are not set to zero then the compiler complains
262 that they ``might be used uninitialized in this
263 function''. Looking at the code this does seem possible, so this
264 should be checked. */
265
266 int notlast;
267
268 size_t n = nc;
269
270next_root:
271
272 if (n == 0) return GSL_SUCCESS;
273
274 iterations = 0;
275
276next_iteration:
277
278 for (e = n; e >= 2; e--) {
279 double a1 = fabs(FMAT(h, e, e - 1, nc));
280 double a2 = fabs(FMAT(h, e - 1, e - 1, nc));
281 double a3 = fabs(FMAT(h, e, e, nc));
282
283 if (a1 <= GSL_DBL_EPSILON * (a2 + a3)) break;
284 }
285
286 x = FMAT(h, n, n, nc);
287
288 if (e == n) {
289 GSL_SET_COMPLEX_PACKED(zroot, n - 1, x + t, 0); /* one real root */
290 n--;
291 goto next_root;
292 /*continue;*/
293 }
294
295 y = FMAT(h, n - 1, n - 1, nc);
296 w = FMAT(h, n - 1, n, nc) * FMAT(h, n, n - 1, nc);
297
298 if (e == n - 1) {
299 p = (y - x) / 2;
300 q = p * p + w;
301 y = sqrt(fabs(q));
302
303 x += t;
304
305 if (q > 0) /* two real roots */
306 {
307 if (p < 0) y = -y;
308 y += p;
309
310 GSL_SET_COMPLEX_PACKED(zroot, n - 1, x - w / y, 0);
311 GSL_SET_COMPLEX_PACKED(zroot, n - 2, x + y, 0);
312 } else {
313 GSL_SET_COMPLEX_PACKED(zroot, n - 1, x + p, -y);
314 GSL_SET_COMPLEX_PACKED(zroot, n - 2, x + p, y);
315 }
316 n -= 2;
317
318 goto next_root;
319 /*continue;*/
320 }
321
322 /* No more roots found yet, do another iteration */
323
324 if (iterations == 60) /* increased from 30 to 60 */
325 {
326 /* too many iterations - give up! */
327
328 return GSL_FAILURE;
329 }
330
331 if (iterations % 10 == 0 && iterations > 0) {
332 /* use an exceptional shift */
333
334 t += x;
335 for (i = 1; i <= n; i++) {
336 FMAT(h, i, i, nc) -= x;
337 }
338
339 s = fabs(FMAT(h, n, n - 1, nc)) + fabs(FMAT(h, n - 1, n - 2, nc));
340 y = 0.75 * s;
341 x = y;
342 w = -0.4375 * s * s;
343 }
344
345 iterations++;
346
347 for (m = n - 2; m >= e; m--) {
348 double a1, a2, a3;
349
350 z = FMAT(h, m, m, nc);
351 r = x - z;
352 s = y - z;
353 p = FMAT(h, m, m + 1, nc) + (r * s - w) / FMAT(h, m + 1, m, nc);
354 q = FMAT(h, m + 1, m + 1, nc) - z - r - s;
355 r = FMAT(h, m + 2, m + 1, nc);
356 s = fabs(p) + fabs(q) + fabs(r);
357 p /= s;
358 q /= s;
359 r /= s;
360
361 if (m == e) break;
362
363 a1 = fabs(FMAT(h, m, m - 1, nc));
364 a2 = fabs(FMAT(h, m - 1, m - 1, nc));
365 a3 = fabs(FMAT(h, m + 1, m + 1, nc));
366
367 if (a1 * (fabs(q) + fabs(r)) <= GSL_DBL_EPSILON * fabs(p) * (a2 + a3))
368 break;
369 }
370
371 for (i = m + 2; i <= n; i++) {
372 FMAT(h, i, i - 2, nc) = 0;
373 }
374
375 for (i = m + 3; i <= n; i++) {
376 FMAT(h, i, i - 3, nc) = 0;
377 }
378
379 /* double QR step */
380
381 for (k = m; k <= n - 1; k++) {
382 notlast = (k != n - 1);
383
384 if (k != m) {
385 p = FMAT(h, k, k - 1, nc);
386 q = FMAT(h, k + 1, k - 1, nc);
387 r = notlast ? FMAT(h, k + 2, k - 1, nc) : 0.0;
388
389 x = fabs(p) + fabs(q) + fabs(r);
390
391 if (x == 0) continue; /* FIXME????? */
392
393 p /= x;
394 q /= x;
395 r /= x;
396 }
397
398 s = sqrt(p * p + q * q + r * r);
399
400 if (p < 0) s = -s;
401
402 if (k != m) {
403 FMAT(h, k, k - 1, nc) = -s * x;
404 } else if (e != m) {
405 FMAT(h, k, k - 1, nc) *= -1;
406 }
407
408 p += s;
409 x = p / s;
410 y = q / s;
411 z = r / s;
412 q /= p;
413 r /= p;
414
415 /* do row modifications */
416
417 for (j = k; j <= n; j++) {
418 p = FMAT(h, k, j, nc) + q * FMAT(h, k + 1, j, nc);
419
420 if (notlast) {
421 p += r * FMAT(h, k + 2, j, nc);
422 FMAT(h, k + 2, j, nc) -= p * z;
423 }
424
425 FMAT(h, k + 1, j, nc) -= p * y;
426 FMAT(h, k, j, nc) -= p * x;
427 }
428
429 j = (k + 3 < n) ? (k + 3) : n;
430
431 /* do column modifications */
432
433 for (i = e; i <= j; i++) {
434 p = x * FMAT(h, i, k, nc) + y * FMAT(h, i, k + 1, nc);
435
436 if (notlast) {
437 p += z * FMAT(h, i, k + 2, nc);
438 FMAT(h, i, k + 2, nc) -= p * r;
439 }
440 FMAT(h, i, k + 1, nc) -= p * q;
441 FMAT(h, i, k, nc) -= p;
442 }
443 }
444
445 goto next_iteration;
446}
447
448static void set_companion_matrix(const double *a, size_t nc, double *m) {
449 size_t i, j;
450
451 for (i = 0; i < nc; i++)
452 for (j = 0; j < nc; j++) MAT(m, i, j, nc) = 0.0;
453
454 for (i = 1; i < nc; i++) MAT(m, i, i - 1, nc) = 1.0;
455
456 for (i = 0; i < nc; i++) MAT(m, i, nc - 1, nc) = -a[i] / a[nc];
457}
458
459static int gsl_poly_complex_solve(const double *a,
460 size_t n,
463 int status;
464 double *m;
465
466 if (n == 0) {
467 cerr << "number of terms must be a positive integer" << endl;
468
469 return (GSL_FAILURE);
470 }
471
472 if (n == 1) {
473 cerr << "cannot solve for only one term" << endl;
474
475 return (GSL_FAILURE);
476 }
477
478 if (a[n - 1] == 0) {
479 cerr << "leading term of polynomial must be non-zero" << endl;
480
481 return (GSL_FAILURE);
482 }
483
484 if (w->nc != n - 1) {
485 cerr << "size of workspace does not match polynomial" << endl;
486
487 return (GSL_FAILURE);
488 }
489
490 m = w->matrix;
491
492 set_companion_matrix(a, n - 1, m);
493
494 balance_companion_matrix(m, n - 1);
495
496 status = qr_companion(m, n - 1, z);
497
498 if (status) {
499 //cerr << "root solving qr method failed to converge" << endl;
500
501 return (GSL_FAILURE);
502 }
503
504 return GSL_SUCCESS;
505}
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1152
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:649
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Matrix class.
Definition: matpackI.h:1261
The Vector class.
Definition: matpackI.h:899
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define q
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
#define a2
#define a1
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric e
Elementary charge convenience name [C].
constexpr Numeric h
Planck constant convenience name [J s].
#define FMAT(m, i, j, n)
Definition: poly_roots.cc:46
#define GSL_DBL_EPSILON
Definition: poly_roots.cc:48
#define GSL_SET_COMPLEX_PACKED(zp, n, x, y)
Definition: poly_roots.cc:58
int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition: poly_roots.cc:90
#define RADIX2
Definition: poly_roots.cc:51
#define MAT(m, i, j, n)
Definition: poly_roots.cc:43
#define RADIX
Definition: poly_roots.cc:50
double * gsl_complex_packed_ptr
Definition: poly_roots.cc:64
#define GSL_SUCCESS
Definition: poly_roots.cc:53
#define GSL_FAILURE
Definition: poly_roots.cc:54
Contains the code to determine roots of polynomials.
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
#define w
#define a
#define c