ARTS  2.4.0(git:4fb77825)
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 
33 // 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 
64 typedef double *gsl_complex_packed_ptr;
65 
66 typedef struct {
67  size_t nc;
68  double *matrix;
70 
71 /* Begin Internal GSL function prototypes */
72 
73 static gsl_poly_complex_workspace *gsl_poly_complex_workspace_alloc(size_t n);
74 
75 static void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w);
76 
77 static void set_companion_matrix(const double *a, size_t n, double *m);
78 
79 static int qr_companion(double *h, size_t nc, gsl_complex_packed_ptr z);
80 
81 static void balance_companion_matrix(double *m, size_t n);
82 
83 static int gsl_poly_complex_solve(const double *a,
84  size_t n,
87 
88 /* End Internal GSL function prototypes */
89 
90 int poly_root_solve(Matrix &roots, Vector &coeffs) {
91  Index a;
92  double *c;
93  double *s;
94 
95  a = coeffs.nelem();
96 
97  assert(roots.nrows() == a - 1);
98  assert(roots.ncols() == 2);
99  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 
132 static 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 
168 static void gsl_poly_complex_workspace_free(gsl_poly_complex_workspace *w) {
169  free(w->matrix);
170  free(w);
171 }
172 
173 static 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 
252 static 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 
270 next_root:
271 
272  if (n == 0) return GSL_SUCCESS;
273 
274  iterations = 0;
275 
276 next_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 
448 static 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 
459 static 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 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
gsl_poly_complex_workspace::matrix
double * matrix
Definition: poly_roots.cc:68
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
GSL_SUCCESS
#define GSL_SUCCESS
Definition: poly_roots.cc:53
GSL_DBL_EPSILON
#define GSL_DBL_EPSILON
Definition: poly_roots.cc:48
ConstMatrixView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1081
poly_roots.h
Contains the code to determine roots of polynomials.
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
RADIX2
#define RADIX2
Definition: poly_roots.cc:51
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
poly_root_solve
int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition: poly_roots.cc:90
GSL_SET_COMPLEX_PACKED
#define GSL_SET_COMPLEX_PACKED(zp, n, x, y)
Definition: poly_roots.cc:58
ConstVectorView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:596
a1
#define a1
Definition: complex.h:56
FMAT
#define FMAT(m, i, j, n)
Definition: poly_roots.cc:46
RADIX
#define RADIX
Definition: poly_roots.cc:50
q
#define q
Definition: legacy_continua.cc:21712
GSL_FAILURE
#define GSL_FAILURE
Definition: poly_roots.cc:54
gsl_poly_complex_workspace
Definition: poly_roots.cc:66
a2
#define a2
Definition: complex.h:58
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
gsl_poly_complex_workspace::nc
size_t nc
Definition: poly_roots.cc:67
Vector
The Vector class.
Definition: matpackI.h:860
gsl_complex_packed_ptr
double * gsl_complex_packed_ptr
Definition: poly_roots.cc:64
MAT
#define MAT(m, i, j, n)
Definition: poly_roots.cc:43