programming6/qrecvc_8c_source.html

00001 /*  qrecvc.c    CCMATH mathematics library source code.
00002  *
00003  *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
00004  *  This code may be redistributed under the terms of the GNU library
00005  *  public license (LGPL). ( See the lgpl.license file for details.)
00006  * ------------------------------------------------------------------------
00007  */
00008 #include "ccmath.h"
00009 void qrecvc(double *ev, Cpx * evec, double *dp, int n)
00010 {
00011     double cc, sc = 0.0, d, x = 0.0, y, h = 0.0, tzr = 1.e-15;
00012
00013     int i, j, k, m, nqr = 50 * n;
00014
00015     Cpx *p;
00016
00017     for (j = 0, m = n - 1; j < nqr; ++j) {
00018         while (1) {
00019             if (m < 1)
00020                 break;
00021             k = m - 1;
00022             if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
00023                 --m;
00024             else {
00025                 x = (ev[k] - ev[m]) / 2.;
00026                 h = sqrt(x * x + dp[k] * dp[k]);
00027                 if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
00028                     break;
00029                 if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
00030                     sc = dp[k] / (2. * cc * h);
00031                 else
00032                     sc = 1.;
00033                 x += ev[m];
00034                 ev[m--] = x - h;
00035                 ev[m--] = x + h;
00036                 for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
00037                     h = p[0].re;
00038                     p[0].re = cc * h + sc * p[n].re;
00039                     p[n].re = cc * p[n].re - sc * h;
00040                     h = p[0].im;
00041                     p[0].im = cc * h + sc * p[n].im;
00042                     p[n].im = cc * p[n].im - sc * h;
00043                 }
00044             }
00045         }
00046         if (x > 0.)
00047             d = ev[m] + x - h;
00048         else
00049             d = ev[m] + x + h;
00050         cc = 1.;
00051         y = 0.;
00052         ev[0] -= d;
00053         for (k = 0; k < m; ++k) {
00054             x = ev[k] * cc - y;
00055             y = dp[k] * cc;
00056             h = sqrt(x * x + dp[k] * dp[k]);
00057             if (k > 0)
00058                 dp[k - 1] = sc * h;
00059             ev[k] = cc * h;
00060             cc = x / h;
00061             sc = dp[k] / h;
00062             ev[k + 1] -= d;
00063             y *= sc;
00064             ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
00065             for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
00066                 h = p[0].re;
00067                 p[0].re = cc * h + sc * p[n].re;
00068                 p[n].re = cc * p[n].re - sc * h;
00069                 h = p[0].im;
00070                 p[0].im = cc * h + sc * p[n].im;
00071                 p[n].im = cc * p[n].im - sc * h;
00072             }
00073         }
00074         ev[k] = ev[k] * cc - y;
00075         dp[k - 1] = ev[k] * sc;
00076         ev[k] = ev[k] * cc + d;
00077     }
00078 }