evmax.c

Go to the documentation of this file.

/*  evmax.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 <stdlib.h>
#include "ccmath.h"
double evmax(double *a, double *u, int n)
{
    double *p, *q, *qm, *r, *s, *t;

    double ev, evm, c, h;

    int kc;

    q = (double *)calloc(n, sizeof(double));
    qm = q + n;
    *(qm - 1) = 1.;
    ev = 0.;
    for (kc = 0; kc < 200; ++kc) {
        h = c = 0.;
        evm = ev;
        for (p = u, r = a, s = q; s < qm;) {
            *p = 0.;
            for (t = q; t < qm;)
                *p += *r++ * *t++;
            c += *p * *p;
            h += *p++ * *s++;
        }
        ev = c / h;
        c = sqrt(c);
        for (p = u, s = q; s < qm;) {
            *p /= c;
            *s++ = *p++;
        }
        if (((c = ev - evm) < 0. ? -c : c) < 1.e-16 * (ev < 0. ? -ev : ev)) {
            free(q);
            return ev;
        }
    }
    free(q);
    for (kc = 0; kc < n;)
        u[kc++] = 0.;
    return 0.;
}