#include <stdio.h>
#include <math.h>
#include "random.h"
#include "numcip.h"
#include "malgbr.h"

double T,m0,del,h,h2;
double *eignewt,*eigQR,**eigvect;
double **mat;
int N;


//Newton's Stuff
void geteignewt(double *eig);
double findbracket(double left);
void evalpn(double, double *);

//QR's Stuff
void geteigQR(double *eig);
void hqr(double **a,int n,double wr[],double wi[]);

//Inverse Power Method for eigenvects
//	we used ludcmp and lubksb here
void getvect(double *eig,double **eigvect);
double normalize(double *vec);

//Form Matrix
double mu(int where);
void fillin(double **mat);

main()
{
	int i,j;
	double scalefac;

	FILE *outfile;

	//system parameters
	T=1000.0;
	m0=0.954;
	del=0.5;
	h=0.05;

	N=(int)(1.0/h);//the grid will have N+1 points [0,N]
	h2=h*h;
	scalefac=T/h2;//actual freq w=sqrt(lambda*T/h2)

	eignewt=dvector(1,N-1);
	eigQR=dvector(1,N-1);
	eigvect=dmatrix(1,N-1,1,N-1);
	mat=dmatrix(1,N-1,1,N-1);//0 and N are at boundary
	for(i=1;i<=N-1;i++)
		eigvect[i][i]=1.0;

	fillin(mat);

	//use Newton's to find all the eigenvals
	geteignewt(eignewt);

	//use inverse power method to find eigenvectors
	getvect(eignewt,eigvect);
	
	//use QR to find all the eigenvals
	geteigQR(eigQR);


	outfile=fopen("eigval.dat","w");
	fprintf(outfile,"Eigenvalues from Newton\n");
	for(i=1;i<N;i++)
		fprintf(outfile,"\tEig%d: %20.15lg\n",i,eignewt[i]);
	fprintf(outfile,"Eigenvalues from QR\n");
	for(i=1;i<N;i++)
		fprintf(outfile,"\tEig%d: %20.15lg\n",i,eigQR[i]);
	fclose(outfile);

	outfile=fopen("eigvect.dat","w");

	for(i=1;i<N-1;i++)
		fprintf(outfile,"vec%d\t",i);
	fprintf(outfile,"vec%d\n",N-1);
	for(i=1;i<N-1;i++)
		fprintf(outfile,"%d\t",0);
	fprintf(outfile,"%d\n",0);
	for(i=1;i<N;i++)
		{
		for(j=1;j<N-1;j++)
			fprintf(outfile,"%20.15lg\t",eigvect[j][i]);
		fprintf(outfile,"%20.15lg\n",eigvect[N-1][i]);
		}
	for(i=1;i<N-1;i++)
		fprintf(outfile,"%d\t",0);
	fprintf(outfile,"%d\n",0);
	fclose(outfile);

	return 0;
}

void fillin(double **m)
//this is a fairly spare matrix.  In an efficient routine, we can only keep track of the nonzero
//	elements
{
	int i;
	double dummy;

	for(i=1;i<N;i++)
		{
		dummy=1.0/mu(i);
		if(i>1)
			m[i][i-1]=-dummy;
		m[i][i]=2.0*dummy;
		if(i<(N-1))
			m[i][i+1]=-dummy;
		}
}

double mu(int where)
{
	return m0+(where*h-0.5)*del;
}

void geteignewt(double *eig)
{
	double olambda,lambda;
	double startval;
	double xstep;
	int i;
	double pn[2];
	int n;
	

//this section was used to get a rough graph of the function
	FILE *outfile2;

	outfile2=fopen("pntest","w");
	startval=0.0;
	n=400;
	xstep=(6-startval)/n;
	for(i=1;i<=n;i++)
	{
		evalpn(startval,pn);
		fprintf(outfile2,"%lg\t%lg\n",startval,pn[0]);
		startval+=xstep;
	}
	fclose(outfile2);
//

	//initial guess at lambda=0
	startval=0.0;

	//find next newton guess
	for(n=1;n<N;n++)
	{
		olambda=findbracket(startval);
		
		do
		{
			lambda=olambda;
			evalpn(lambda,pn);//pn[0]=pn(lambda) and pn[1]=dpn(lambda)
			olambda=lambda-(pn[0]/pn[1]);

		}
		while(fabs((olambda-lambda)/lambda)>1e-7);
		eig[n]=olambda;

		startval=eig[n]+0.005;
	}
}

double findbracket(double left)
//this will start with left value and move forward to look for f(x) changing sign
//	then it will return the right point after crossing zero as seed for Newton
//the increment steps are estimated from graph of f(x)
{
	double fleft[2],fright[2];
	double right;
	
	do
	{
		right=left+0.01;
		evalpn(left,fleft);
		evalpn(right,fright);
		left=right;

	}
	while(fleft[0]*fright[0]>0);

	return right;
}

void evalpn(double x,double *poly)
//this use the recursion relation to evaluate the characteristic polynomial 
//	and its derivative for the tridiagonl A-Ix
{
	int i;
	double ppn,pn1,pn2;
	double dpn,dpn1,dpn2;
	double bn,ac;

	//initialize
	pn2=1.0;
	pn1=mat[1][1]-x;
	dpn2=0.0;
	dpn1=-1.0;

	for(i=2;i<N;i++)
		{
		//coef for iteration
		bn=mat[i][i]-x;
		ac=mat[i][i-1]*mat[i-1][i];
		
		//iteration
		ppn=bn*pn1-ac*pn2;
		dpn=-pn1+bn*dpn1-ac*dpn2;

		//replace old values
		pn2=pn1;
		pn1=ppn;
		dpn2=dpn1;
		dpn1=dpn;
		}

	//passing values back
	poly[0]=ppn;
	poly[1]=dpn;
}

void geteigQR(double *eig)
//on output eig will contain the eigenvalues
//	and mat's columns will contain the eigenvectors
{

	double *eigi;
	double max;
	int imax;
	int i,j;

	eigi=dvector(1,N-1);

	hqr(mat,N-1,eig,eigi);

	//sort the eigs in accending order
	for(i=1;i<N;i++)
	{
		max=0.0;
		for(j=1;j<=N-i;j++)
			if(eig[j]>max)
			{
				imax=j;
				max=eig[j];
			}
		eig[imax]=eig[N-i];
		eig[N-i]=max;
	}

	free_dvector(eigi,1,N-1);
}

void hqr(double **a,int n,double wr[],double wi[])
{
	int nn,m,l,k,j,its,i,mmin;
	double z,y,x,w,v,u,t,s,r,q,p,anorm;

	anorm=fabs(a[1][1]);
	for (i=2;i<=n;i++)
		for (j=(i-1);j<=n;j++)
			anorm += fabs(a[i][j]);
	nn=n;
	t=0.0;
	while (nn >= 1) {
		its=0;
		do {
			for (l=nn;l>=2;l--) {
				s=fabs(a[l-1][l-1])+fabs(a[l][l]);
				if (s == 0.0) s=anorm;
				if ((double)(fabs(a[l][l-1]) + s) == s) break;
			}
			x=a[nn][nn];
			if (l == nn) {
				wr[nn]=x+t;
				wi[nn--]=0.0;
			} else {
				y=a[nn-1][nn-1];
				w=a[nn][nn-1]*a[nn-1][nn];
				if (l == (nn-1)) {
					p=0.5*(y-x);
					q=p*p+w;
					z=sqrt(fabs(q));
					x += t;
					if (q >= 0.0) {
						z=p+SIGN(z,p);
						wr[nn-1]=wr[nn]=x+z;
						if (z) wr[nn]=x-w/z;
						wi[nn-1]=wi[nn]=0.0;
					} else {
						wr[nn-1]=wr[nn]=x+p;
						wi[nn-1]= -(wi[nn]=z);
					}
					nn -= 2;
				} else {
					if (its == 30) nrerror("Too many iterations in hqr");
					if (its == 10 || its == 20) {
						t += x;
						for (i=1;i<=nn;i++) a[i][i] -= x;
						s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
						y=x=0.75*s;
						w = -0.4375*s*s;
					}
					++its;
					for (m=(nn-2);m>=l;m--) {
						z=a[m][m];
						r=x-z;
						s=y-z;
						p=(r*s-w)/a[m+1][m]+a[m][m+1];
						q=a[m+1][m+1]-z-r-s;
						r=a[m+2][m+1];
						s=fabs(p)+fabs(q)+fabs(r);
						p /= s;
						q /= s;
						r /= s;
						if (m == l) break;
						u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
						v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
						if ((double)(u+v) == v) break;
					}
					for (i=m+2;i<=nn;i++) {
						a[i][i-2]=0.0;
						if (i != (m+2)) a[i][i-3]=0.0;
					}
					for (k=m;k<=nn-1;k++) {
						if (k != m) {
							p=a[k][k-1];
							q=a[k+1][k-1];
							r=0.0;
							if (k != (nn-1)) r=a[k+2][k-1];
							if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) {
								p /= x;
								q /= x;
								r /= x;
							}
						}
						if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {
							if (k == m) {
								if (l != m)
								a[k][k-1] = -a[k][k-1];
							} else
								a[k][k-1] = -s*x;
							p += s;
							x=p/s;
							y=q/s;
							z=r/s;
							q /= p;
							r /= p;
							for (j=k;j<=nn;j++) {
								p=a[k][j]+q*a[k+1][j];
								if (k != (nn-1)) {
									p += r*a[k+2][j];
									a[k+2][j] -= p*z;
								}
								a[k+1][j] -= p*y;
								a[k][j] -= p*x;
							}
							mmin = nn<k+3 ? nn : k+3;
							for (i=l;i<=mmin;i++) {
								p=x*a[i][k]+y*a[i][k+1];
								if (k != (nn-1)) {
									p += z*a[i][k+2];
									a[i][k+2] -= p*r;
								}
								a[i][k+1] -= p*q;
								a[i][k] -= p;
							}
						}
					}
				}
			}
		} while (l < nn-1);
	}
}

void getvect(double *eig,double **eigvect)
//will use power method to find the eigenmode correspondng eignval x
//(A-Ix)y_n+1=y_n can be solved by ludcmp/lubksb without explicitly using (A-Ix)^-1
{
	int seed=-1;//seed for random number generator
	double **aq;
	int *luindx;
	double lud;
	int n=N-1;
	int i,j,k;
	int M;

	aq=dmatrix(1,n,1,n);
	luindx=ivector(1,n);
	
	for(k=1;k<=n;k++)
	{
		for(i=1;i<=n;i++)
		{
			for(j=1;j<=i;j++)
				aq[i][j]=mat[i][j];
			aq[i][i]=mat[i][i]-eig[k];
			for(j=i+1;j<=n;j++)
				aq[i][j]=mat[i][j];
		}
		//initialize eigenvects in random directions
		for(i=1;i<N;i++)
			eigvect[k][i]=ran1(&seed);

		//begin inverse power method iteration
		ludcmp(aq,n,luindx,&lud);
		lubksb(aq,n,luindx,eigvect[k]);
	
		normalize(eigvect[k]);
		M=0;
		do
		{
			lubksb(aq,n,luindx,eigvect[k]);
			normalize(eigvect[k]);
			M++;
		}
		while(M<3);
	}

	free_dmatrix(aq,1,N-1,1,N-1);
	free_ivector(luindx,1,N-1);
}


double normalize(double *vec)
{
	int i;
	double temp;

	temp=vec[1]*vec[1];
	for(i=2;i<=N-1;i++)
		temp+=vec[i]*vec[i];
	temp=sqrt(temp);
	
	for(i=1;i<=N-1;i++)
		vec[i]=vec[i]/temp;

	return temp;
}