/*  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 <math.h>
#include <stdlib.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.;
}