ARTS  2.2.66
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 
18 
20 // File description
22 
36 // External declarations
39 
40 #include <iostream>
41 #include <cmath>
42 #include <stdexcept>
43 #include "poly_roots.h"
44 
45 
46 /* C-style matrix elements */
47 #define MAT(m,i,j,n) ((m)[(i)*(n) + (j)])
48 
49 /* Fortran-style matrix elements */
50 #define FMAT(m,i,j,n) ((m)[((i)-1)*(n) + ((j)-1)])
51 
52 
53 #define GSL_DBL_EPSILON 2.2204460492503131e-16
54 
55 #define RADIX 2
56 #define RADIX2 (RADIX*RADIX)
57 
58 #define GSL_SUCCESS 0
59 #define GSL_FAILURE -1
60 #define GSL_EINVAL 4
61 #define GSL_EFAILED 5
62 
63 #define GSL_SET_COMPLEX_PACKED(zp,n,x,y) do {*((zp)+2*(n))=(x); *((zp)+(2*(n)+1))=(y);} while(0)
64 
65 
66 typedef double *gsl_complex_packed_ptr ;
67 
68 typedef struct {
69  size_t nc;
70  double * matrix;
72 
73 
74 /* Begin Internal GSL function prototypes */
75 
77 gsl_poly_complex_workspace_alloc (size_t n);
78 
79 static void
80 gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w);
81 
82 static void
83 set_companion_matrix (const double *a, size_t n, double *m);
84 
85 static int
86 qr_companion (double *h, size_t nc, gsl_complex_packed_ptr z);
87 
88 static void
89 balance_companion_matrix (double *m, size_t n);
90 
91 static int
92 gsl_poly_complex_solve (const double *a, size_t n,
95 
96 /* End Internal GSL function prototypes */
97 
98 
99 int
100 poly_root_solve (Matrix& roots, Vector& coeffs)
101 {
102  Index a;
103  double *c;
104  double *s;
105 
106  a = coeffs.nelem();
107 
108  assert (roots.nrows() == a - 1);
109  assert (roots.ncols() == 2);
110  assert (coeffs[a - 1] != 0);
111 
112 #ifdef USE_DOUBLE
113  c = coeffs.mdata;
114  s = roots.mdata;
115 #else
116  c = (double *)malloc (a * sizeof (double));
117  for (Index i = 0; i < a; i++)
118  c[i] = (double)coeffs[i];
119 
120  s = (double *)malloc ((a-1) * 2 * sizeof (double));
121 #endif
122 
123  gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (a);
124 
125  int status = gsl_poly_complex_solve (c, a, w, s);
126 
127 #ifndef USE_DOUBLE
128  if (status == GSL_SUCCESS)
129  {
130  for (Index i = 0; i < a - 1; i++)
131  {
132  roots (i, 0) = (Numeric)s[i * 2];
133  roots (i, 1) = (Numeric)s[i * 2 + 1];
134  }
135  }
136 
137  free (c);
138  free (s);
139 #endif
140 
141  gsl_poly_complex_workspace_free (w);
142 
143  return status;
144 }
145 
146 
148 gsl_poly_complex_workspace_alloc (size_t n)
149 {
150  size_t nc ;
151 
153 
154  if (n == 0)
155  {
156  cerr << "matrix size n must be positive integer" << endl;
157 
158  return (NULL);
159  }
160 
162  malloc (sizeof(gsl_poly_complex_workspace));
163 
164  if (w == 0)
165  {
166  cerr << "failed to allocate space for struct" << endl;
167 
168  return (NULL);
169  }
170 
171  nc = n - 1;
172 
173  w->nc = nc;
174 
175  w->matrix = (double *) malloc (nc * nc * sizeof(double));
176 
177  if (w->matrix == 0)
178  {
179  free (w) ; /* error in constructor, avoid memory leak */
180 
181  cerr << "failed to allocate space for workspace matrix" << endl;
182 
183  return (NULL);
184  }
185 
186  return w ;
187 }
188 
189 
190 static void
191 gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w)
192 {
193  free(w->matrix) ;
194  free(w);
195 }
196 
197 
198 static void
199 balance_companion_matrix (double *m, size_t nc)
200 {
201  int not_converged = 1;
202 
203  double row_norm = 0;
204  double col_norm = 0;
205 
206  while (not_converged)
207  {
208  size_t i, j;
209  double g, f, s;
210 
211  not_converged = 0;
212 
213  for (i = 0; i < nc; i++)
214  {
215  /* column norm, excluding the diagonal */
216 
217  if (i != nc - 1)
218  {
219  col_norm = fabs (MAT (m, i + 1, i, nc));
220  }
221  else
222  {
223  col_norm = 0;
224 
225  for (j = 0; j < nc - 1; j++)
226  {
227  col_norm += fabs (MAT (m, j, nc - 1, nc));
228  }
229  }
230 
231  /* row norm, excluding the diagonal */
232 
233  if (i == 0)
234  {
235  row_norm = fabs (MAT (m, 0, nc - 1, nc));
236  }
237  else if (i == nc - 1)
238  {
239  row_norm = fabs (MAT (m, i, i - 1, nc));
240  }
241  else
242  {
243  row_norm = (fabs (MAT (m, i, i - 1, nc))
244  + fabs (MAT (m, i, nc - 1, nc)));
245  }
246 
247  if (col_norm == 0 || row_norm == 0)
248  {
249  continue;
250  }
251 
252  g = row_norm / RADIX;
253  f = 1;
254  s = col_norm + row_norm;
255 
256  while (col_norm < g)
257  {
258  f *= RADIX;
259  col_norm *= RADIX2;
260  }
261 
262  g = row_norm * RADIX;
263 
264  while (col_norm > g)
265  {
266  f /= RADIX;
267  col_norm /= RADIX2;
268  }
269 
270  if ((row_norm + col_norm) < 0.95 * s * f)
271  {
272  not_converged = 1;
273 
274  g = 1 / f;
275 
276  if (i == 0)
277  {
278  MAT (m, 0, nc - 1, nc) *= g;
279  }
280  else
281  {
282  MAT (m, i, i - 1, nc) *= g;
283  MAT (m, i, nc - 1, nc) *= g;
284  }
285 
286  if (i == nc - 1)
287  {
288  for (j = 0; j < nc; j++)
289  {
290  MAT (m, j, i, nc) *= f;
291  }
292  }
293  else
294  {
295  MAT (m, i + 1, i, nc) *= f;
296  }
297  }
298  }
299  }
300 }
301 
302 
303 static int
304 qr_companion (double *h, size_t nc, gsl_complex_packed_ptr zroot)
305 {
306  double t = 0.0;
307 
308  size_t iterations, e, i, j, k, m;
309 
310  double w, x, y, s, z;
311 
312  double p = 0, q = 0, r = 0;
313 
314  /* FIXME: if p,q,r, are not set to zero then the compiler complains
315  that they ``might be used uninitialized in this
316  function''. Looking at the code this does seem possible, so this
317  should be checked. */
318 
319  int notlast;
320 
321  size_t n = nc;
322 
323 next_root:
324 
325  if (n == 0)
326  return GSL_SUCCESS ;
327 
328  iterations = 0;
329 
330 next_iteration:
331 
332  for (e = n; e >= 2; e--)
333  {
334  double a1 = fabs (FMAT (h, e, e - 1, nc));
335  double a2 = fabs (FMAT (h, e - 1, e - 1, nc));
336  double a3 = fabs (FMAT (h, e, e, nc));
337 
338  if (a1 <= GSL_DBL_EPSILON * (a2 + a3))
339  break;
340  }
341 
342  x = FMAT (h, n, n, nc);
343 
344  if (e == n)
345  {
346  GSL_SET_COMPLEX_PACKED (zroot, n-1, x + t, 0); /* one real root */
347  n--;
348  goto next_root;
349  /*continue;*/
350  }
351 
352  y = FMAT (h, n - 1, n - 1, nc);
353  w = FMAT (h, n - 1, n, nc) * FMAT (h, n, n - 1, nc);
354 
355  if (e == n - 1)
356  {
357  p = (y - x) / 2;
358  q = p * p + w;
359  y = sqrt (fabs (q));
360 
361  x += t;
362 
363  if (q > 0) /* two real roots */
364  {
365  if (p < 0)
366  y = -y;
367  y += p;
368 
369  GSL_SET_COMPLEX_PACKED (zroot, n-1, x - w / y, 0);
370  GSL_SET_COMPLEX_PACKED (zroot, n-2, x + y, 0);
371  }
372  else
373  {
374  GSL_SET_COMPLEX_PACKED (zroot, n-1, x + p, -y);
375  GSL_SET_COMPLEX_PACKED (zroot, n-2, x + p, y);
376  }
377  n -= 2;
378 
379  goto next_root;
380  /*continue;*/
381  }
382 
383  /* No more roots found yet, do another iteration */
384 
385  if (iterations == 60) /* increased from 30 to 60 */
386  {
387  /* too many iterations - give up! */
388 
389  return GSL_FAILURE ;
390  }
391 
392  if (iterations % 10 == 0 && iterations > 0)
393  {
394  /* use an exceptional shift */
395 
396  t += x;
397  for (i = 1; i <= n; i++)
398  {
399  FMAT (h, i, i, nc) -= x;
400  }
401 
402  s = fabs (FMAT (h, n, n - 1, nc)) + fabs (FMAT (h, n - 1, n - 2, nc));
403  y = 0.75 * s;
404  x = y;
405  w = -0.4375 * s * s;
406  }
407 
408  iterations++;
409 
410  for (m = n - 2; m >= e; m--)
411  {
412  double a1, a2, a3;
413 
414  z = FMAT (h, m, m, nc);
415  r = x - z;
416  s = y - z;
417  p = FMAT (h, m, m + 1, nc) + (r * s - w) / FMAT (h, m + 1, m, nc);
418  q = FMAT (h, m + 1, m + 1, nc) - z - r - s;
419  r = FMAT (h, m + 2, m + 1, nc);
420  s = fabs (p) + fabs (q) + fabs (r);
421  p /= s;
422  q /= s;
423  r /= s;
424 
425  if (m == e)
426  break;
427 
428  a1 = fabs (FMAT (h, m, m - 1, nc));
429  a2 = fabs (FMAT (h, m - 1, m - 1, nc));
430  a3 = fabs (FMAT (h, m + 1, m + 1, nc));
431 
432  if (a1 * (fabs (q) + fabs (r)) <= GSL_DBL_EPSILON * fabs (p) * (a2 + a3))
433  break;
434  }
435 
436  for (i = m + 2; i <= n; i++)
437  {
438  FMAT (h, i, i - 2, nc) = 0;
439  }
440 
441  for (i = m + 3; i <= n; i++)
442  {
443  FMAT (h, i, i - 3, nc) = 0;
444  }
445 
446  /* double QR step */
447 
448  for (k = m; k <= n - 1; k++)
449  {
450  notlast = (k != n - 1);
451 
452  if (k != m)
453  {
454  p = FMAT (h, k, k - 1, nc);
455  q = FMAT (h, k + 1, k - 1, nc);
456  r = notlast ? FMAT (h, k + 2, k - 1, nc) : 0.0;
457 
458  x = fabs (p) + fabs (q) + fabs (r);
459 
460  if (x == 0)
461  continue; /* FIXME????? */
462 
463  p /= x;
464  q /= x;
465  r /= x;
466  }
467 
468  s = sqrt (p * p + q * q + r * r);
469 
470  if (p < 0)
471  s = -s;
472 
473  if (k != m)
474  {
475  FMAT (h, k, k - 1, nc) = -s * x;
476  }
477  else if (e != m)
478  {
479  FMAT (h, k, k - 1, nc) *= -1;
480  }
481 
482  p += s;
483  x = p / s;
484  y = q / s;
485  z = r / s;
486  q /= p;
487  r /= p;
488 
489  /* do row modifications */
490 
491  for (j = k; j <= n; j++)
492  {
493  p = FMAT (h, k, j, nc) + q * FMAT (h, k + 1, j, nc);
494 
495  if (notlast)
496  {
497  p += r * FMAT (h, k + 2, j, nc);
498  FMAT (h, k + 2, j, nc) -= p * z;
499  }
500 
501  FMAT (h, k + 1, j, nc) -= p * y;
502  FMAT (h, k, j, nc) -= p * x;
503  }
504 
505  j = (k + 3 < n) ? (k + 3) : n;
506 
507  /* do column modifications */
508 
509  for (i = e; i <= j; i++)
510  {
511  p = x * FMAT (h, i, k, nc) + y * FMAT (h, i, k + 1, nc);
512 
513  if (notlast)
514  {
515  p += z * FMAT (h, i, k + 2, nc);
516  FMAT (h, i, k + 2, nc) -= p * r;
517  }
518  FMAT (h, i, k + 1, nc) -= p * q;
519  FMAT (h, i, k, nc) -= p;
520  }
521  }
522 
523  goto next_iteration;
524 }
525 
526 
527 static void
528 set_companion_matrix (const double *a, size_t nc, double *m)
529 {
530  size_t i, j;
531 
532  for (i = 0; i < nc; i++)
533  for (j = 0; j < nc; j++)
534  MAT (m, i, j, nc) = 0.0;
535 
536  for (i = 1; i < nc; i++)
537  MAT (m, i, i - 1, nc) = 1.0;
538 
539  for (i = 0; i < nc; i++)
540  MAT (m, i, nc - 1, nc) = -a[i] / a[nc];
541 }
542 
543 
544 static int
545 gsl_poly_complex_solve (const double *a, size_t n,
548 {
549  int status;
550  double *m;
551 
552  if (n == 0)
553  {
554  cerr << "number of terms must be a positive integer" << endl;
555 
556  return (GSL_FAILURE);
557  }
558 
559  if (n == 1)
560  {
561  cerr << "cannot solve for only one term" << endl;
562 
563  return (GSL_FAILURE);
564  }
565 
566  if (a[n - 1] == 0)
567  {
568  cerr << "leading term of polynomial must be non-zero" << endl;
569 
570  return (GSL_FAILURE);
571  }
572 
573  if (w->nc != n - 1)
574  {
575  cerr << "size of workspace does not match polynomial" << endl;
576 
577  return (GSL_FAILURE);
578  }
579 
580  m = w->matrix;
581 
582 
583  set_companion_matrix (a, n - 1, m);
584 
585  balance_companion_matrix (m, n - 1);
586 
587  status = qr_companion (m, n - 1, z);
588 
589  if (status)
590  {
591  //cerr << "root solving qr method failed to converge" << endl;
592 
593  return (GSL_FAILURE);
594  }
595 
596  return GSL_SUCCESS;
597 }
598 
599 
600 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
#define q
Definition: continua.cc:21469
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:667
Contains the code to determine roots of polynomials.
#define GSL_SET_COMPLEX_PACKED(zp, n, x, y)
Definition: poly_roots.cc:63
#define FMAT(m, i, j, n)
Definition: poly_roots.cc:50
The Vector class.
Definition: matpackI.h:556
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
#define RADIX
Definition: poly_roots.cc:55
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
#define GSL_FAILURE
Definition: poly_roots.cc:59
double * gsl_complex_packed_ptr
Definition: poly_roots.cc:66
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define MAT(m, i, j, n)
Definition: poly_roots.cc:47
The Matrix class.
Definition: matpackI.h:788
#define GSL_DBL_EPSILON
Definition: poly_roots.cc:53
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:358
#define GSL_SUCCESS
Definition: poly_roots.cc:58
int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition: poly_roots.cc:100
#define RADIX2
Definition: poly_roots.cc:56
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832