programming6/qrevec_8c_source.html

00001 /*  qrevec.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 <math.h>
00009 int qrevec(double *ev, double *evec, double *dp, int n)
00010 {
00011     double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
00012
00013     int i, j, k, m, mqr = 8 * n;
00014
00015     double *p;
00016
00017     for (j = 0, m = n - 1;; ++j) {
00018         while (1) {
00019             if (m < 1)
00020                 return 0;
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];
00038                     p[0] = cc * h + sc * p[n];
00039                     p[n] = cc * p[n] - sc * h;
00040                 }
00041             }
00042         }
00043         if (j > mqr)
00044             return -1;
00045         if (x > 0.)
00046             d = ev[m] + x - h;
00047         else
00048             d = ev[m] + x + h;
00049         cc = 1.;
00050         y = 0.;
00051         ev[0] -= d;
00052         for (k = 0; k < m; ++k) {
00053             x = ev[k] * cc - y;
00054             y = dp[k] * cc;
00055             h = sqrt(x * x + dp[k] * dp[k]);
00056             if (k > 0)
00057                 dp[k - 1] = sc * h;
00058             ev[k] = cc * h;
00059             cc = x / h;
00060             sc = dp[k] / h;
00061             ev[k + 1] -= d;
00062             y *= sc;
00063             ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
00064             for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
00065                 h = p[0];
00066                 p[0] = cc * h + sc * p[n];
00067                 p[n] = cc * p[n] - sc * h;
00068             }
00069         }
00070         ev[k] = ev[k] * cc - y;
00071         dp[k - 1] = ev[k] * sc;
00072         ev[k] = ev[k] * cc + d;
00073     }
00074     return 0;
00075 }