diff options
Diffstat (limited to 'gsl-1.9/poly/qr.c')
| -rw-r--r-- | gsl-1.9/poly/qr.c | 245 |
1 files changed, 245 insertions, 0 deletions
diff --git a/gsl-1.9/poly/qr.c b/gsl-1.9/poly/qr.c new file mode 100644 index 0000000..23e7bfb --- /dev/null +++ b/gsl-1.9/poly/qr.c @@ -0,0 +1,245 @@ +/* poly/qr.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +static int qr_companion (double *h, size_t nc, gsl_complex_packed_ptr z); + +static int +qr_companion (double *h, size_t nc, gsl_complex_packed_ptr zroot) +{ + double t = 0.0; + + size_t iterations, e, i, j, k, m; + + double w, x, y, s, z; + + double p = 0, q = 0, r = 0; + + /* FIXME: if p,q,r, are not set to zero then the compiler complains + that they ``might be used uninitialized in this + function''. Looking at the code this does seem possible, so this + should be checked. */ + + int notlast; + + size_t n = nc; + +next_root: + + if (n == 0) + return GSL_SUCCESS ; + + iterations = 0; + +next_iteration: + + for (e = n; e >= 2; e--) + { + double a1 = fabs (FMAT (h, e, e - 1, nc)); + double a2 = fabs (FMAT (h, e - 1, e - 1, nc)); + double a3 = fabs (FMAT (h, e, e, nc)); + + if (a1 <= GSL_DBL_EPSILON * (a2 + a3)) + break; + } + + x = FMAT (h, n, n, nc); + + if (e == n) + { + GSL_SET_COMPLEX_PACKED (zroot, n-1, x + t, 0); /* one real root */ + n--; + goto next_root; + /*continue;*/ + } + + y = FMAT (h, n - 1, n - 1, nc); + w = FMAT (h, n - 1, n, nc) * FMAT (h, n, n - 1, nc); + + if (e == n - 1) + { + p = (y - x) / 2; + q = p * p + w; + y = sqrt (fabs (q)); + + x += t; + + if (q > 0) /* two real roots */ + { + if (p < 0) + y = -y; + y += p; + + GSL_SET_COMPLEX_PACKED (zroot, n-1, x - w / y, 0); + GSL_SET_COMPLEX_PACKED (zroot, n-2, x + y, 0); + } + else + { + GSL_SET_COMPLEX_PACKED (zroot, n-1, x + p, -y); + GSL_SET_COMPLEX_PACKED (zroot, n-2, x + p, y); + } + n -= 2; + + goto next_root; + /*continue;*/ + } + + /* No more roots found yet, do another iteration */ + + if (iterations == 60) /* increased from 30 to 60 */ + { + /* too many iterations - give up! */ + + return GSL_FAILURE ; + } + + if (iterations % 10 == 0 && iterations > 0) + { + /* use an exceptional shift */ + + t += x; + + for (i = 1; i <= n; i++) + { + FMAT (h, i, i, nc) -= x; + } + + s = fabs (FMAT (h, n, n - 1, nc)) + fabs (FMAT (h, n - 1, n - 2, nc)); + y = 0.75 * s; + x = y; + w = -0.4375 * s * s; + } + + iterations++; + + for (m = n - 2; m >= e; m--) + { + double a1, a2, a3; + + z = FMAT (h, m, m, nc); + r = x - z; + s = y - z; + p = FMAT (h, m, m + 1, nc) + (r * s - w) / FMAT (h, m + 1, m, nc); + q = FMAT (h, m + 1, m + 1, nc) - z - r - s; + r = FMAT (h, m + 2, m + 1, nc); + s = fabs (p) + fabs (q) + fabs (r); + p /= s; + q /= s; + r /= s; + + if (m == e) + break; + + a1 = fabs (FMAT (h, m, m - 1, nc)); + a2 = fabs (FMAT (h, m - 1, m - 1, nc)); + a3 = fabs (FMAT (h, m + 1, m + 1, nc)); + + if (a1 * (fabs (q) + fabs (r)) <= GSL_DBL_EPSILON * fabs (p) * (a2 + a3)) + break; + } + + for (i = m + 2; i <= n; i++) + { + FMAT (h, i, i - 2, nc) = 0; + } + + for (i = m + 3; i <= n; i++) + { + FMAT (h, i, i - 3, nc) = 0; + } + + /* double QR step */ + + for (k = m; k <= n - 1; k++) + { + notlast = (k != n - 1); + + if (k != m) + { + p = FMAT (h, k, k - 1, nc); + q = FMAT (h, k + 1, k - 1, nc); + r = notlast ? FMAT (h, k + 2, k - 1, nc) : 0.0; + + x = fabs (p) + fabs (q) + fabs (r); + + if (x == 0) + continue; /* FIXME????? */ + + p /= x; + q /= x; + r /= x; + } + + s = sqrt (p * p + q * q + r * r); + + if (p < 0) + s = -s; + + if (k != m) + { + FMAT (h, k, k - 1, nc) = -s * x; + } + else if (e != m) + { + FMAT (h, k, k - 1, nc) *= -1; + } + + p += s; + x = p / s; + y = q / s; + z = r / s; + q /= p; + r /= p; + + /* do row modifications */ + + for (j = k; j <= n; j++) + { + p = FMAT (h, k, j, nc) + q * FMAT (h, k + 1, j, nc); + + if (notlast) + { + p += r * FMAT (h, k + 2, j, nc); + FMAT (h, k + 2, j, nc) -= p * z; + } + + FMAT (h, k + 1, j, nc) -= p * y; + FMAT (h, k, j, nc) -= p * x; + } + + j = (k + 3 < n) ? (k + 3) : n; + + /* do column modifications */ + + for (i = e; i <= j; i++) + { + p = x * FMAT (h, i, k, nc) + y * FMAT (h, i, k + 1, nc); + + if (notlast) + { + p += z * FMAT (h, i, k + 2, nc); + FMAT (h, i, k + 2, nc) -= p * r; + } + FMAT (h, i, k + 1, nc) -= p * q; + FMAT (h, i, k, nc) -= p; + } + } + + goto next_iteration; +} + |