qreval.c

Go to the documentation of this file.

/*  qreval.c    CCMATH mathematics library source code.
 *
 *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
 *  This code may be redistributed under the terms of the GNU library
 *  public license (LGPL). ( See the lgpl.license file for details.)
 * ------------------------------------------------------------------------
 */
#include "ccmath.h"
int qreval(double *ev, double *dp, int n)
{
    double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
 
    int j, k, m, mqr = 8 * n;
 
    for (j = 0, m = n - 1;; ++j) {
        while (1) {
            if (m < 1)
                return 0;
            k = m - 1;
            if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
                --m;
            else {
                x = (ev[k] - ev[m]) / 2.;
                h = sqrt(x * x + dp[k] * dp[k]);
                if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
                    break;
                x += ev[m];
                ev[m--] = x - h;
                ev[m--] = x + h;
            }
        }
        if (j > mqr)
            return -1;
        if (x > 0.)
            d = ev[m] + x - h;
        else
            d = ev[m] + x + h;
        cc = 1.;
        y = 0.;
        ev[0] -= d;
        for (k = 0; k < m; ++k) {
            x = ev[k] * cc - y;
            y = dp[k] * cc;
            h = sqrt(x * x + dp[k] * dp[k]);
            if (k > 0)
                dp[k - 1] = sc * h;
            ev[k] = cc * h;
            cc = x / h;
            sc = dp[k] / h;
            ev[k + 1] -= d;
            y *= sc;
            ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
        }
        ev[k] = ev[k] * cc - y;
        dp[k - 1] = ev[k] * sc;
        ev[k] = ev[k] * cc + d;
    }
    return 0;
}