//
// String.cpp is designed to examine the oridinary differential equation problem 
// given by a vibrating string fixed at the ends.
// Using matrix algebra we can create a discrete matrix representing the finite
// difference approximation of the string's wave equation.  This method will 
// be used to to approximate the eigenvalues and eigenvectors of the wave
// equation for a non uniform string.  
//

#include <algorithm>
#include <vector>
#include <map>
#include "Function.h" 
#include <math.h>


typedef vector<vector<double> > ddvector ;

class Matrix;
void svdcmp(Matrix& a, int m, int n, vector<double>& w, Matrix&  v, int maxIterations);

class Matrix : public ddvector
{

public:
	Matrix(int a, int b) : vector<vector<double> >(a), nRows(a), nCols(b)
	{
		for (int i =0;i<a;i++)
			(*this)[i].resize(b);
	}

	Matrix(vector<double>& a) : vector<vector<double> > ((size_type)a.size()), nCols(1), nRows(a.size())
	{
		for (int i =0;i<a.size();i++)
			(*this)[i][0] = a[i];
	}

	Matrix(ddvector::const_iterator b, ddvector::const_iterator e, int firstColumn, int numberOfColumns)
		: ddvector(b,e), nCols(numberOfColumns)
	{//tbd
	}

	/*Matrix(double a)   : ddvector(1) , nRows(1), nCols(1)  // be wary of type conversions
	{
		(*this)[0].resize(1);
		(*this)[0][0] = a;
	}*/
	
	void subMatrix(Matrix m, int first_row, int numOfRows, int firstCol, int numOfCols)
	{//tbd
	}

	//overloaded function keeps track of size
	//probably should get size from vector class
	//but this causes trouble with instatiating matrices of unkown size
	virtual iterator erase(iterator _P) {nRows--;return ddvector::erase(_P);};  
	
	Matrix minor(int row, int col)
	{
		Matrix m(*this);
		ddvector::iterator row_iter = m.begin();
		for (int i = 0; i< m.numOfRows()-1;i++, row_iter++)
		{
			if (i== row)
				m.erase(row_iter);
			else
			{			
				vector<double>::iterator col_iter = (*row_iter).begin();
				for (int c=0;c<col;c++, col_iter++);
				(*row_iter).erase(col_iter);
			}
		};
		m.setsize(nRows-1,nCols-1);
		return m;
	    //copy(map.begin(), map2.end(), insert_iterator<Matrix> (map1, map1.begin())
	}

	//recursive alogrithm grinds to a halt for matrices with rank greater than 7
	double det_recursive(int row=0, int col=0)
	{
		if (nRows != nCols) cerr << "Determinant can not be found for non square matrix [" << nRows << "][" << nCols << "]" << endl;
		double sum = (nRows==1) ? (*this)[0][0] : 0;
		for (int r = 0;r < nRows;r++)
			if ((r)%2)
					sum += (*this)[row][col]*minor(row,col).det();
				else
					sum -= (*this)[row][col]*minor(row,col).det();
		return sum;
	}
	



	virtual double det()
	{
   		double d = 1;

		vector<double> w(nRows);
		Matrix v(nCols, nCols);

		svdcmp(*this, nRows, nCols, w, v, 30);

		for (int i = 0; i < nRows; i++)
		  d *= w[i];
		//cout << "The determinant was calculated using the vector \n";
		//for (i = 0; i < nRows; i++)
		//	cout <<i << "\t" <<  w[i] << endl;;
		//cout << "# The determinant is : " << d << endl;
		return d;
	}

	Matrix operator+(const Matrix& m)
	{
		Matrix ret(nRows, nCols);
		for (int i=0; i < m.size(); i++)
			for (int j=0;j < m[i].size(); j++)
				ret[i][j] = m[i][j]+(*this)[i][j];
		return ret;
	};

	Matrix operator~()  //transpose
	{
		Matrix ret(nCols,nRows);
		for (int i=0; i < nRows; i++)
			for (int j=0;j < nCols; j++)
				ret[j][i] = (*this)[i][j];
		return ret;
	};

	Matrix operator*(const double& d)
	{
		Matrix ret(numOfRows(), numOfColumns());
		for (int i=0; i < ret.numOfRows(); i++)
		{
			for (int j=0;j < ret.numOfColumns(); j++)
			{
				ret[i][j] = d*(*this)[i][j];
			}
		}
		return ret;
	}
	Matrix operator*(const Matrix& m)
	{
		Matrix ret(nRows, m.numOfColumns());
		for (int i=0; i < ret.numOfRows(); i++)
		{
			for (int j=0;j < ret.numOfColumns(); j++)
			{
				double sum = 0;
				for (int k = 0; k < nCols;k++)
					sum += m[k][j]*(*this)[i][k];
				ret[i][j] = sum;
			}
		}
		return ret;
	}

	Matrix operator-(const Matrix& m)
	{
		Matrix ret(nRows, nCols);
		for (int i=0; i < ret.numOfRows(); i++)
		{
			for (int j=0;j < ret.numOfColumns(); j++)
			{
				ret[i][j] = (*this)[i][j] - m[i][j];
			}
		}
		return ret;
	}

	virtual void setsize(int a, int b)
	{
		resize(a);
		for (int i=0;i<a;i++)
			(*this)[i].resize(b);
		nRows = a;
		nCols = b;
	}

	vector<double> getColumn(int colNum)
	{
		vector<double> ret(nRows);
		for (int i=0;i<nRows;i++)
			ret[i] = (*this)[i][colNum];
		return ret;
	};
		
	void gaussJordanElimination(Matrix& multipliers, vector<int>& interchanges)
	{
		for (int i = 0 ; i < nCols;i++)
		{
			for (int j = i+1 ; j < nRows; j++)
			{
				//if ((*this)[i][i] == 0) //what to do.. what to do...
				pivot(multipliers, j , nRows);  //stolen from posted solution 
				multipliers[j][i] = (*this)[j][i]/(*this)[i][i]; //multipliers of ith row that is subtracted from the kth row during the ith step of the reduction process
				for (int k=i;k< nCols;k++)
					(*this)[j][k] -= multipliers[j][i]*(*this)[i][k];
			}
		};
		//cout << "After gauss-jordan elimination this is \n" << *this << endl;
		//cout << "The multipliers were \n" << multipliers << endl;
	}

	void gaussJordan(Matrix& solution)
	{
		
		Matrix multipliers(nRows, nCols);
		vector<int> interchanges;
		int i, j;

		gaussJordanElimination(multipliers, interchanges);

		//forward substitution
		for (i = 0;i<solution.numOfRows();i++)
			for (j=0;j<i;j++)
				solution[i][0] -= multipliers[i][j]*solution[j][0];

		//cout << "solution after forward substitution\n" << solution << endl;
		//back substitution
		for (i = solution.numOfRows()-1; i>=0;i--)
		{
			for (j = i+1;j< solution.numOfRows();j++)
				solution[i][0] -= (*this)[i][j]*solution[j][0];
			solution[i][0] /=(*this)[i][i];
		}
	};

	//from posted solution 
	void pivot(Matrix& solution,int level,int n)
	//this will find the largest pivot and rotate row
	//	to find the largest, it will only look downward
	{
		int i,maxi;
		double maxval;
		
		//find max pivot
		maxval=fabs((*this)[level][level]);
		maxi=level;
		for(i=level+1;i<n;i++)
			{
			if(fabs((*this)[i][level])>maxval)
				{
				maxval=(*this)[i][level];
				maxi=i;
				}
			}

		//rotate rows
		if(maxi!=level)
		{
			for(i=level;i<n;i++)
				{
				maxval=(*this)[level][i];
				(*this)[level][i]=(*this)[maxi][i];
				(*this)[maxi][i]=maxval;
				}
			maxval=solution[level][0];
			solution[level][0]=solution[maxi][0];
			solution[maxi][0]=maxval;
		}
	}


	//householder reduction 
	void householder()
	{
		;
	};

	//uses a similarity transform (such as a householder reduction) 
	//to turn A into an upper Hessenberg Matrix
	void upperHessenberg();

	//returns the eigenvalues of A
	//A should be an upper Hessenberg Matrix
	vector<double> qr()
	{
		vector<double> ret(nRows);

		return ret;
	}

	//nr function
	double SIGN(double a, double b)
	{
		return (b >=0 ? (a >= 0 ? a : -a)  : (a >= 0 ? -a : a));
	}
	//adapted from the solution
	void hqr(vector<double>& wr,vector<double>& wi)
	{
		Matrix& a = (*this);
		int n = a.numOfRows();
		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[0][0]);
		for (i=1;i<n;i++)
			for (j=(i-1);j<n;j++)
				anorm += fabs(a[i][j]);
		nn=n-1;
		t=0.0;
		while (nn >= 0) {
			its=0;
			do {
				for (l=nn;l>0;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) cerr << ("Too many iterations in hqr");
						if (its == 10 || its == 20) {
							t += x;
							for (i=0;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 ludcmp(Matrix &a, vector<int>& indx, double& d)
	{
		double Tiny=1e-20;
		int i, imax, j, k;
		double big, dum, sum, temp;

		int n=a. numOfRows ( );
		vector<double> vv(n);
		d=1.0;
		for (i=0; i<n;i++) 
		{
		   big=0.0;
			for (j=0; j<n; j++)
				if ((temp=fabs(a[i] [j])) > big) big=temp;
			if (big == 0.0) cerr << "Singular matrix in routine ludcmp";
			vv [i]=1.0/big;
		}
					
		for (j=0; j<n; j++) 
		{
		   for(i=0; i<j; i++) 
		   {
			   sum=a[i] [j];
			   for (k=0; k<i; k++) sum -= a [i] [k]*a [k] [j];
			   a [i] [j]=sum;
		   }
		   big=0.0;
		   for (i=j; i<n;i++) 
		   {
			   sum=a [i] [j];
			   for (k=0; k<j; k++) sum -= a [i] [k] *a [k] [j];
			   a [i] [j]=sum;
			   if (( dum=vv [i]*fabs (sum)) >=big)
			   {
				   //Is the figure of merit for the pivot better than the best so far?
				   big=dum;
				   imax=i;
			   }
		   }
		   if (j  != imax) 
		   { 
			   //Do we need to change rows?
			   for (k=0; k<n; k++) 
			   {
				   dum =a [imax] [k];
				   a [imax] [k] =a [j] [k];
				   a [j] [k] = dum;
			   }
			   d=-d;
			   vv[imax]=vv [j];
		   }
		   indx[j] = imax;
		   if (a[j][j] == 0.0) a[j][j] = Tiny;
		   if (j != n-1) 
		   {
			   dum = 1.0/a[j][j];
			   for (i=j+1;i<n;i++) a[i][j] *= dum;
		   }
	   }
	}
	
	void lubksb(Matrix& a, const vector<int>& indx, vector<double>& b)
	{
		int i , ii=0, ip, j;
		double sum;

		int n = a.numOfRows();
		for (i=0;i<n;i++)
		{
			ip = indx[i];
			sum = b[ip];
			b[ip]=b[i];
			if (ii != 0)
				for (j=ii=1;j<i;j++) sum -= a[i][j]*b[j];
			else if (sum != 0.0)
				ii=i+1;
			b[i]=sum;
		}
		for (i=n-1;i>=0;i--)
		{
			sum = b[i];
			for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
			b[i] = sum/a[i][i];
		}
	}

//	}


   //returns the eigenvectors
   	//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

	Matrix eigenVectors(vector<double>& eig)
	{
		//int RAND_MAX = (~0);
		int seed=-1;//seed for random number generator
		int maxIterations = 4;
		Matrix& mat = (*this);
		Matrix aq(nRows,nCols);
		vector<int> luindx(nRows);
		double lud;
		int n=nRows;
		int i,j,k;
		Matrix eigvect(nRows, nCols);
		
		for(k=0;k<n;k++)
		{
			for(i=0;i<n;i++)
			{
				for(j=0;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=0;i<n;i++)
				eigvect[k][i]=((double) rand()/(double)RAND_MAX)*2 - 1;

			//cout << "# The random eigenvector is given by : \n" << eigvect << endl;
			//begin inverse power method iteration
			ludcmp(aq,luindx,lud);
			//cout << "# The eigenvector after LU decomposition is : " << aq << endl;
			
			lubksb(aq,luindx,eigvect[k]);
			//cout << "# The eigenvector after " << 1 << " backsubstitutions is given by : \n" << eigvect << endl;
			normalize(eigvect[k]);

			for (i=0;i<maxIterations;i++)
			{
				//cout << "# The normalized eigenvector after " << i+1 << " backsubstitutions is given by : \n" << eigvect << endl;
				lubksb(aq,luindx,eigvect[k]);
				//cout << "# The eigenvector after " << i+1 << " backsubstitutions is given by : \n" << eigvect << endl;
				normalize(eigvect[k]);
			}
		}

		return eigvect;
	}


	double normalize(vector<double>& vec)
	{
		int i;
		double temp;

		temp=vec[0]*vec[0];
		for(i=1;i< vec.size();i++)
			temp+=vec[i]*vec[i];
		temp=sqrt(temp);
		
		for(i=0;i< vec.size();i++)
			vec[i]=vec[i]/temp;

		return temp;
	}


	friend ostream& operator<<(ostream& out, const Matrix& m)
	{
		for (int i =0; i < m.size();i++)
		{
			out << "#" <<  i << "\t";
			for (int j = 0;j < m[i].size();j++)
				out << ((double)m[i][j]) << "\t";
			out  << endl;
		};
		return out;
	};
	int numOfRows() const {return nRows;}
	int numOfColumns() const {return nCols;};
protected:
	int nRows;
	int nCols;
};

//taken from Press (;-)shhh
void svdbksb(Matrix& u, vector<double>& w, Matrix& v, Matrix& b, vector<double>& x)
{
	int jj, j, i;
	double s;

	int m = u.numOfRows();
	int n = u.numOfColumns();

	vector<double> tmp(n);
	for (j=0;j<n;j++)
	{
		s=0.0;
		if (w[j])
		{
			for (i=0;i<m;i++) s +=- u[i][j]*b[i][0];
			s /=w[j];
		}
		tmp[j] = s;
	}
	for (j=0;j<n;j++)
	{
		s=0;
		for (jj=0;jj<n;jj++) s += v[j][jj]*tmp[jj];
		x[j] = s;
	}
};


//SVD algorithms taken from sourceforge's linalg, a Simple Linear Algebra C++ Library 

/*
PYTHAG computes sqrt(a^{2} + b^{2}) without destructive overflow or underflow.
*/
static double at, bt, ct;
#define PYTHAG(a, b) ((at = fabs(a)) > (bt = fabs(b)) ? \
(ct = bt/at, at*sqrt(1.0+ct*ct)): (bt ? (ct = at/bt, bt*sqrt(1.0+ct*ct)): 0.0))

static double maxarg1, maxarg2;
#define MAX(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? \
(maxarg1) : (maxarg2))

#define SIGN(a, b) ((b) < 0.0 ? -fabs(a): fabs(a))

void error(const char error_text[])
  /* Numerical Recipes standard error handler.			*/
{
  cerr << "Numerical Recipes run-time error..." << endl;
  cerr << error_text << endl;
  throw runtime_error(error_text);
}

/*
Given a matrix a[m][n], this routine computes its singular value
decomposition, A = U*W*V^{T}.  The matrix U replaces a on output.
The diagonal matrix of singular values W is output as a vector w[n].
The matrix V (not the transpose V^{T}) is output as v[n][n].
m must be greater or equal to n;  if it is smaller, then a should be
filled up to square with zero rows.
*/
void svdcmp(Matrix& a, int m, int n, vector<double>& w, Matrix&  v, int maxIterations)
  {
    int flag, i, its, j, jj, k;
    int l  = 0;
    int nm = 0;
    double c, f, h, s, x, y, z;
    double anorm = 0.0, g = 0.0, scale = 0.0;

    if (m < n)
      cerr << "SVDCMP: Matrix A must be augmented with extra rows of zeros.";
    double* rv1 = new double [n];

    /* Householder reduction to bidiagonal form.			*/
    for (i = 0; i < n; i++)
      {
	l = i + 1;
        rv1[i] = scale*g;
        g = s = scale = 0.0;
        if (i < m)
	  {
	    for (k = i; k < m; k++)
	      scale += fabs(a[k][i]);
	    if (scale)
	      {
	        for (k = i; k < m; k++)
	          {
		    a[k][i] /= scale;
		    s += a[k][i]*a[k][i];
	          };
	        f = a[i][i];
	        g = -SIGN(sqrt(s), f);
	        h = f*g - s;
	        a[i][i] = f - g;
	        if (i != n - 1)
	          {
		    for (j = l; j < n; j++)
		      {
		        for (s  = 0.0, k = i; k < m; k++)
		          s += a[k][i]*a[k][j];
		        f = s/h;
		        for ( k = i; k < m; k++)
			  a[k][j] += f*a[k][i];
		      };
	          };
	        for (k = i; k < m; k++)
	          a[k][i] *= scale;
	      };
          };
        w[i] = scale*g;
        g = s= scale = 0.0;
        if (i < m && i != n - 1)
          {
	    for (k = l; k < n; k++)
	      scale += fabs(a[i][k]);
	    if (scale)
	      {
	        for (k = l; k < n; k++)
	          {
		    a[i][k] /= scale;
		    s += a[i][k]*a[i][k];
	          };
	        f = a[i][l];
	        g = -SIGN(sqrt(s), f);
	        h = f*g - s;
	        a[i][l] = f - g;
	        for (k = l; k < n; k++)
	          rv1[k] = a[i][k]/h;
	        if (i != m - 1)
	          {
		    for (j = l; j < m; j++)
		      {
		        for (s = 0.0, k = l; k < n; k++)
		          s += a[j][k]*a[i][k];
		        for (k = l; k < n; k++)
		          a[j][k] += s*rv1[k];
		      };
	          };
	        for (k = l; k < n; k++)
	          a[i][k] *= scale;
	      };
	  };
        anorm = MAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
      };
    /* Accumulation of right-hand transformations.			*/
    for (i = n - 1; 0 <= i; i--)
      {
	if (i < n - 1)
	  {
	    if (g)
	      {
		for (j = l; j < n; j++)
		  v[j][i] = (a[i][j]/a[i][l])/g;
		  /* Double division to avoid possible underflow:	*/
		for (j = l; j < n; j++)
		  {
		    for (s = 0.0, k = l; k < n; k++)
		      s += a[i][k]*v[k][j];
		    for (k = l; k < n; k++)
		      v[k][j] += s*v[k][i];
		  };
	      };
	    for (j = l; j < n; j++)
	      v[i][j] = v[j][i] = 0.0;
	  };
	v[i][i] = 1.0;
	g = rv1[i];
	l = i;
      };
    /* Accumulation of left-hand transformations.			*/
    for (i = n - 1; 0 <= i; i--)
      {
	l = i + 1;
	g = w[i];
	if (i < n - 1)
	  for (j = l; j < n; j++)
	    a[i][j] = 0.0;
	if (g)
	  {
	    g = 1.0/g;
	    if (i != n - 1)
	      {
		for (j = l; j < n; j++)
		  {
		    for (s = 0.0, k = l; k < m; k++)
		      s += a[k][i]*a[k][j];
		    f = (s/a[i][i])*g;
		    for (k = i; k < m; k++)
		      a[k][j] += f*a[k][i];
		  };
	       };
	    for (j = i; j < m; j++)
	      a[j][i] *= g;
	  }
	else
	  for (j = i; j < m; j++)
	    a[j][i] = 0.0;
	++a[i][i];
      };
    /* Diagonalization of the bidiagonal form.				*/
    for (k = n - 1; 0 <= k; k--)	/* Loop over singular values.	*/
      {
	for (its = 0; its < maxIterations; its++)	/* Loop over allowed iterations.*/
	  {
	    flag = 1;
	    for (l = k; 0 <= l; l--)	/* Test for splitting:		*/
	      {
		nm = l - 1;		/* Note that rv1[0] is always zero.*/
		if (fabs(rv1[l]) + anorm == anorm)
		  {
		    flag = 0;
		    break;
		  };
		if (fabs(w[nm]) + anorm == anorm)
		  break;
	      };
	    if (flag)
	      {
		c = 0.0;		/* Cancellation of rv1[l], if l>0:*/
		s = 1.0;
		for (i = l; i <= k; i++) {
		    f = s*rv1[i];
		    if (fabs(f) + anorm != anorm)
		      {
			g = w[i];
			h = PYTHAG(f, g);
			w[i] = h;
			h = 1.0/h;
			c = g*h;
			s = (-f*h);
			for (j = 0; j < m; j++)
			  {
			    y = a[j][nm];
			    z = a[j][i];
			    a[j][nm] = y*c + z*s;
			    a[j][i]  = z*c - y*s;
			  };
		      };
		  };
	      };
	    z = w[k];
	    if (l == k)		/* Convergence.				*/
	      {
		if (z < 0.0)	/* Singular value is made non-negative.	*/
		  {
		    w[k] = -z;
		    for (j = 0; j < n; j++)
		      v[j][k] = (-v[j][k]);
		  };
		break;
	      };
	    if (its == (maxIterations - 1))
	      error("No convergence in SVDCMP iterations.");
	    x = w[l];		/* Shift from bottom 2-by-2 minor.	*/
	    nm = k - 1;
	    y = w[nm];
	    g = rv1[nm];
	    h = rv1[k];
	    f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
	    g = PYTHAG(f, 1.0);
	    f = ((x - z)*(x + z) + h*((y/(f + SIGN(g, f))) - h))/x;
	    /* Next QR transformation:					*/
	    c = s = 1.0;
	    for (j = l; j <= nm; j++)
	      {
		i = j + 1;
		g = rv1[i];
		y = w[i];
		h = s*g;
		g = c*g;
		z = PYTHAG(f, h);
		rv1[j] = z;
		c = f/z;
		s = h/z;
		f = x*c + g*s;
		g = g*c - x*s;
		h = y*s;
		y = y*c;
		for (jj = 0; jj < n;  jj++)
		  {
		    x = v[jj][j];
		    z = v[jj][i];
		    v[jj][j] = x*c + z*s;
		    v[jj][i] = z*c - x*s;
		  };
		z = PYTHAG(f, h);
		w[j] = z;	/* Rotation can be arbitrary if z = 0.	*/
		if (z)
		  {
		    z = 1.0/z;
		    c = f*z;
		    s = h*z;
		  };
		f = (c*g) + (s*y);
		x = (c*y) - (s*g);
		for (jj = 0; jj < m; jj++)
		  {
		    y = a[jj][j];
		    z = a[jj][i];
		    a[jj][j] = y*c + z*s;
		    a[jj][i] = z*c - y*s;
		  };
	      };
	    rv1[l] = 0.0;
	    rv1[k] = f;
	    w[k] = x;
	  };
      };
    delete [] rv1;
  }  


//typedef MatrixFunction class Function<Matrix>;

class ScalarMatrix : public Matrix
{
public:
	ScalarMatrix() : Matrix(1,1) {};
	ScalarMatrix(double value) : Matrix(1,1) { setValue(value);};
	void setValue(double value) { (*this)[0][0]=value;};
	//operator double() { return (*this)[0][0];};
};

class SquareMatrix : public Matrix
{
public:
	SquareMatrix(int n) : Matrix(n,n) {};
	SquareMatrix(SquareMatrix& m) : Matrix(m) {};
	int getSize() {return nRows;};
	
};

class SecularEquation : public Function
{
public:
	SecularEquation(SquareMatrix& A, Matrix& Y) :a(A), y(Y), identity(a.getSize())
	{
		for (int i=0;i<a.numOfRows();i++)
			identity[i][i]=1;
	};

	Matrix operator()(Matrix& ev)
	{
		double eigenvalue = ev[0][0];
		Matrix ret = (a-(identity*eigenvalue));
		return ret;
	};
protected:
	SquareMatrix a;
	Matrix y;
	SquareMatrix identity;
};

class TridiagonalMatrix : public Matrix
{
public:
	TridiagonalMatrix(Matrix& m) : Matrix(m) {};
	virtual double det() 
	{
		double p0 = 1, p1 = (*this)[0][0];
		double pn = p1;
		for (int i = 0; i < nRows-1;i++)
		{
			pn = (*this)[i+1][i+1]*p1 - (*this)[i][i+1]*(*this)[i+1][i]*p0;
			p0 = p1;
			p1 = pn;
			//cout << "Polynomial Approximation: " << pn  << endl;
		};
		return pn;
	}
};

class Approximation
{
public:
	Approximation(double lowerBound, double upperBound, Function& function)
		:a(lowerBound), b(upperBound), f(function)
	{
		if (a > b) cerr << "Approximation lower bound is greater than upper bound" << endl;
	};

	virtual double calculate(double precision, int maxTries = 3000)
	{
		cout << "Default approximation not implemented" << endl;
		return -1;
	}

protected:
	double a;
	double b;
	Function& f;
};

//Bisection Approximation for the degenerate case
//
// bisects twice and chooses 3/4 of the interval around the least value
//
// has not been proven
//
class BisectionDegenerateApproximation : public Approximation
{	
public:
	BisectionDegenerateApproximation(double lower, double upper, Function& f)
		: Approximation(lower,upper,f)
	{};
	double calculate(double precision=0.0, int maxTries=2000)
	{
		double p=precision+1;
		while ((p > precision) && (--maxTries))
		{
			if ((f((a+(b-a)/4)))< (f(a+3*(b-a)/4)))
				b = a+3*(b-a)/4;
			else 
				a = a+(b-a)/4;
			p = fabs(a-b);
		}
		return (a+b)/2;
	};
};

//Bisection Approximation
class Bisection : public Approximation
{	
public:
	Bisection(double lower, double upper, Function& f)
		: Approximation(lower,upper,f)
	{};
	double calculate(double precision=0.0, int maxTries=2000)
	{
		double p=precision+1;
		double x = a, y=b;
		double z=(y-x)/2.0;
		while ((z > precision) && (--maxTries))
		{
			if ((f(x+z)*f(x)) > 0)
				x += z;
			else
				y = x + z;
			z /= 2.0;
		}
		if (!maxTries) cerr << "# maxTries reached in Bisection Approximation" << endl;
		return (x+y)/2.0;
	};
};

class SecularDeterminant : public Function //secular determinant
{
public:
	SecularDeterminant(SecularEquation& eq) : secular(eq) {};
	virtual double operator()(double ev) 
	{
		ScalarMatrix eigenvalue(ev);
		TridiagonalMatrix secular_matrix = secular(eigenvalue);
		return secular_matrix.det();
		//cout << ev << "\t" << secular_matrix.det() << "\t" << sm.det() << endl;
	};	
protected:
	SecularEquation& secular;
};

int main()
{
	double tension = 1000; //Newtons
	double length = 1;          //meters
	double mue0 = 0.954;   //grams/meter (average mass density)
	double delta = 0.5;         //grams/meter^2  (mass density variations per unit length)
	int gridsize = 20;             // number of points on the discretized grid
	int n = gridsize/length;  // number of rows and columns in the design matrix
	double h = 1./(double)n;						 // length of the increments on the grid
	int i,j,k;							 //counters
	SquareMatrix a(n-1);					//tridiagonal matrix containing difference equations
	Matrix x( n-1,1);                  //solutions vector of the difference equations (eigenvector)

	double ev;                      //eigenvalue of matrix equation (negative of angular frequencey squared)
	double ev_lower=0, ev_upper=5, ev_inc=.01;

	Function self;
	Function& mue = ((self-length/2)*delta + mue0);

	//construct tridiagonal design matrix 
	double c = tension/(h*h);		//coefficient on every term factored out 
	double xi = length/n;
	for (i = 0; i < n - 2; i++,xi += length/n)
	{
		
		cout << "# mue(" << xi << ") = " << ((self-length/2)*delta + mue0)(xi) << " or " <<  mue(xi) << endl;;
		a[i][i] = 2/mue(xi);
		a[i+1][i] = -1/mue(xi + length/n);
		a[i][i+1] = -1/mue(xi);
	};
	a[i][i] = 2/mue((i+1)*length/n);


	cout << "# The design matrix is given by :\n" << a << endl;
	Matrix designMatrix = a;
	Matrix designMatrix2 = a;


	//plot the det(A-lambda(1))
	SecularEquation secular(a, x);
	cout << "\n# Plot of possible eigenvalue vs. the determinant of the secular equation \n";
	for (ev = ev_lower+ev_inc; ev < ev_upper;ev += ev_inc)
	{
		ScalarMatrix eigenvalue(ev);
   		TridiagonalMatrix secular_matrix = secular(eigenvalue);
		Matrix sm = secular(eigenvalue);
		//cout << "The secular_matrix is given by\n" << secular_matrix << endl;
		cout << ev << "\t" << secular_matrix.det() << "\t" << sm.det() << endl;
	};

	//find the zero's of the plot to determine the eigenvalues
	SecularDeterminant sD(secular);
	Bisection eigenvalue0(0,0.05, sD);
	Bisection eigenvalue1(0.05,0.2, sD);
	Bisection eigenvalue2(0.2,.4, sD);
	double precision = 0.00000001;
	double e0,e1,e2;
	cout.precision(14);
	cout << "# The eigenvalues without the common factor were found to be " << endl;
	cout << "# e0 = " << (e0 = eigenvalue0.calculate(precision, 20)) << endl;;
	cout << "# e1 = " << (e1 = eigenvalue1.calculate(precision, 20)) << endl;;
	cout << "# e2 = " << (e2 = eigenvalue2.calculate(precision, 20)) << endl;;
	cout << "# The eigenfrequencies are given by " << endl;
	cout << "# w0 = " << sqrt(e0*c) << endl;
	cout << "# w1 = " << sqrt(e1*c) << endl;
	cout << "# w2 = " << sqrt(e2*c) << endl;

	//use QR to determine all the eignenvalues
	vector<double> w1(n-1), w2(n-1);
	cout << "# The Eigenvalues given by QR decomposition algorithm are : " << endl;
	designMatrix.hqr(w1, w2);
	sort(w1.begin(), w1.end());
	for (i = 0; i< w1.size();i++)
		cout << "# " << w1[i] << endl;

	Matrix eigenVectors = designMatrix2.eigenVectors(w1);
	cout << "\n# The eigenvectors are give by the matrix " << endl;
	cout << eigenVectors;

	cout << "# The plot of the matrix is given by " << endl;
	for (i = 0; i < eigenVectors.numOfRows(); i++)
	{
		cout << i  << "\t";
		for (j=0;j< eigenVectors.numOfColumns(); j++)
			cout << eigenVectors[j][i] << "\t";
		cout << endl;
	}
   
	return 0;
}