$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 1998-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

/*

*/


#ifdef _WIN32
  #include <windows.h>
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <math.h>
#include <iostream>

#include "matrixutil.h"

using namespace std;
using namespace fq;

//----- Local functions ( optional declaration ) 
void balanc(double** a, int n);
void elmhes0(double** a, int n);		// TODO - not yet used 
void elmhes(double** a, int n);
void hmges(int n, double* w, double** v, double* x);
BOOL hqr(double** a, int n, double* wr, double* wi);
double pythag(double a, double b);
void swap(double *a, double *b, int iItems);


// ----------------------------------------------------------
// Swap two doubles - swap(&n,&n) function or SWAP(n,n) macro
// you must define locally:
// double swap_nTemp;
#define SWAP(g,h) {swap_nTemp=(g);(g)=(h);(h)=swap_nTemp;}

#define sign(a,b)   ((b)>=0.0 ? fabs(a) : -fabs(a))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

void swap(double *a, double *b, int iItems){
	double tmp;
	for( int i=0; i<iItems; ++i ) {
		tmp=a[i];
		a[i]=b[i];
		b[i]=tmp;
		}
	}

/*	---------------------------------------------------------------
	Gauss-Jordan linear equation solver
	2 versions - 0 based arrays and 1 based arrays
	a[1..n][1..n] is the input matrix
	b[1..n][1..m] is input containing the m right-hand side vectors
	on output, a is replaced by its matrix inverse,
	and b is replaced by the corresponding set of solution vectors.
*/
/*
BOOL gaussj(double** a, int n, double **b, int m, ostream &msgOut){
	int i, icol=-1, irow=-1, j, k, l, ll;
	double big,dum,pivinv;
	int * indxc = int_alloc1D(n+1);
	int * indxr = int_alloc1D(n+1);
	int * ipiv  = int_alloc1D(n+1);
	
	for (j=1; j<=n; j++) 
		ipiv[j]=0;

	for (i=1; i<=n; i++) {
		irow = icol = -1; // should be inside i ?
		big=0.0;
		for (j=1; j<=n; j++) {
			if (ipiv[j] != 1) {
				for (k=1; k<=n; k++) {
					if (ipiv[k] == 0) {
						if (fabs(a[j][k]) >= big) {
							big=fabs(a[j][k]);
							irow=j;
							icol=k;
							}
						} 
					else if (ipiv[k]>1) {
						msgOut << "Fatal err: Singular matrix in gaussj" << endl;
						free(ipiv);
						free(indxr);
						free(indxc);
						return FALSE;
						}
					}
				}
			}
		
		if ( irow == -1 || icol == -1 ) {
			msgOut << "gaussj found no pivot" << endl;
			free(ipiv);
			free(indxr);
			free(indxc);
			return FALSE;
			}
		
		++(ipiv[icol]);
		if (irow != icol) {
			swap(a[irow], a[icol], n+1);	// SWAP ROWS
			if( b != NULL ) 
				swap(b[irow], b[icol], m+1);
			}

		indxr[i]=irow;
		indxc[i]=icol;
		if (a[icol][icol] == 0.0) {
			msgOut << "Fatal err: Singular matrix in gaussj" << endl;
			free(ipiv);
			free(indxr);
			free(indxc);
			return FALSE;
			}
		
		pivinv=1.0/a[icol][icol];
		a[icol][icol]=1.0;
		for (l=1; l<=n; l++) 
			a[icol][l] *= pivinv;
		for (l=1; l<=m; l++) 
			b[icol][l] *= pivinv;
		for (ll=1; ll<=n; ll++) {
			if (ll != icol) {
				dum=a[ll][icol];
				a[ll][icol]=0.0;
				for (l=1; l<=n; l++) 
					a[ll][l] -= a[icol][l]*dum;
				for (l=1; l<=m; l++) 
					b[ll][l] -= b[icol][l]*dum;
				}
			}
		}
	
	for (l=n; l>=1; l--)
		if (indxr[l] != indxc[l])
			for (k=1; k<=n; k++)
				SWAP( (a[k][indxr[l]]), (a[k][indxc[l]]) );	//SWAP COLUMNS
			
	free(ipiv);
	free(indxr);
	free(indxc);
	return TRUE;
	}
*/

//--  Similar to Gauss-J above.  Removes memory leaks.  Uses 0 based arrays.
//--  Added comments.  More meaningful variable names.  Better Swap usage.
extern "C" QUBFAST_API BOOL gaussj0(double** a, int n, double **b, int m, ostream &msgOut){
	int iIter, iRow, iCol, iMaxCol, iMaxRow;	
	double 	big, dum, pivinv;
	int * indxc = int_alloc1D(n);	// [i] is the column of the i'th pivot element
	int * indxr = int_alloc1D(n);	// [i] is the original row for the i'th pivot 
	int * ipiv = int_alloc1D(n);
	double swap_nTemp;
	
	//-- Set all pivots to 0 ( false ) 
	for( iRow=0; iRow<n; iRow++) 
		ipiv[iRow]=0;

	iMaxRow = iMaxCol = -1;

	//-- For each iteration, find max and {divide, swap}.
	for (iIter=0; iIter<n; iIter++) {
		//-- For i,j not yet pivoted, find where max abs(a[i][j]) occurs ==> Pivot 
		big=0.0;
		for (iRow=0; iRow<n; iRow++) {
			if (ipiv[iRow] != 1) {
				for (iCol=0; iCol<n; iCol++) {
					if (ipiv[iCol] == 0) { 
						if (fabs(a[iRow][iCol]) > big)
							big=fabs(a[iMaxRow=iRow][iMaxCol=iCol]);
						}
					else if (ipiv[iCol]>1) {
						msgOut << "Fatal err: Singular matrix in gauss-j" << endl;
						free(ipiv);
						free(indxr);
						free(indxc);
						return FALSE;
						}
					}
				}
			}
		if ( iMaxRow == -1 || iMaxCol == -1 ) {
		  // msgOut << "gauss-j found no pivot" << endl;
			free(ipiv);
			free(indxr);
			free(indxc);
			return FALSE;
			}

		//-- Got the pivot element, interchange rows to force the pivot element to 
		//-- be on the diagonal.  
		++(ipiv[iMaxCol]);
		if (iMaxRow != iMaxCol) {
			swap( a[iMaxRow], a[iMaxCol], n );
			if( b != NULL ) 
				swap( b[iMaxRow], b[iMaxCol], m );
			}

		indxr[iIter]=iMaxRow;
		indxc[iIter]=iMaxCol;
		pivinv = a[iMaxCol][iMaxCol];
		a[iMaxCol][iMaxCol]=1.0;		//-- this is unusual, but is fixed by the next if stmt.
		//-- divide row iCol by big in a and b
		for (iCol=0; iCol<n; iCol++) 
			a[iMaxCol][iCol] /= pivinv;
		for (iCol=0; iCol<m; iCol++) 
			b[iMaxCol][iCol] /= pivinv;

		//-- For each row, subtract a multiple of the maxrow which makes maxcol element 0
		for (iRow=0; iRow<n; iRow++) {
			if (iRow != iMaxCol) {
				dum=a[iRow][iMaxCol];
				a[iRow][iMaxCol]=0.0;
				for (iCol=0; iCol<n; iCol++) 
					a[iRow][iCol] -= a[iMaxCol][iCol]*dum;
				for (iCol=0; iCol<m; iCol++) 
					b[iRow][iCol] -= b[iMaxCol][iCol]*dum;
				}
			}
		}
	
	// unscramble the solution w.r.t. the column interchanges
	// by interchanging columns in the reverse order 
	for( iIter=n-1; iIter>=0; iIter--)
		if (indxr[iIter] != indxc[iIter])
			for (iRow=0; iRow<n; iRow++)
				SWAP( (a[iRow][indxr[iIter]]), (a[iRow][indxc[iIter]]) );

	free(ipiv);
	free(indxr);
	free(indxc);
	return TRUE;
	}

extern "C" QUBFAST_API BOOL gaussj_invert(matrix<double>& mat, int n, std::ostream& msgOut){
	return gaussj0(mat, n, NULL, 0, msgOut );
	}

/*	---------------------------------------------------------------
	Pythagoras theorem without destructive underflow or overflow.
	Used by svdecomp function.
*/
double pythag(double a, double b){
	a=fabs(a);
	b=fabs(b);
	if( a>b )
		return a*sqrt(1.0+DSQR(b/a));
	if ( b==0.0 )		
		return 0.0;
	return b*sqrt(1.0+DSQR(a/b));
	}

/*	---------------------------------------------------------------
	Singular-value decomposition routine
*/

extern "C" QUBFAST_API int svdecomp(double **a, int m, int n, double *w, double **v, ostream &msgOut){
	int     flag,i,its,j,jj,k,l,nm;
	double  anorm,c,f,g,h,s,scale,x,y,z,*rv1;
	int rtn = 0;

	rv1 = double_alloc1D(n+1);
	for (i=0; i<=n; ++i)
	  rv1[i] = 0.0;
	
	// householder reduction to bidiagonal form 
	g=scale=anorm=0.0;
	for (i=1; 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;
				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) {
			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;
				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=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
		}

	//-- accumulation of right hand transformations 
	for (i=n; i>=1; i--) {
		if (i < n) {
			if (g) {
				for (j=l; j<=n; j++) v[j][i]=(a[i][j]/a[i][l])/g;
				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 right hand transformations
	for (i=IMIN(m,n); i>=1; i--) {
		l=i+1;
		g=w[i];
		for (j=l; j<=n; j++) a[i][j]=0.0;
		if (g) {
			g=1.0/g;
			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: 
	//   loop over singular values and allowed iterations
	for (k=n; k>=1; k--) {
		for (its=1; its<=30; its++) {
			flag=1;
			for (l=k; l>=1; l--) {		// test for splitting
				nm=l-1;					// note that rv1[1] is always zero...
				if (fabs(rv1[l])+anorm == anorm) {
					flag=0;
					break;
					}
				if (fabs(w[nm])+anorm == anorm)	// ...w[nm] could err if rv1[1]!=0
					break;
				}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l; i<=k; i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if (fabs(f)+anorm == anorm) 
						break;
					g=w[i];
					h=pythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=1; 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.  Singular value is made non negative.
				if (z < 0.0) {
					w[k] = -z;
					for (j=1; j<=n; j++) 
						v[j][k] = -v[j][k];
					}
				break;
				}
			if(its==50) {
				rtn = 1;
				msgOut << "Err: no convergence in svdecomp" << endl;
				}
			
			// shift from bottom 2-by-2 minor
			x=w[l];
			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;
			c=s=1.0;					// Next QR transformation
			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 *= c;
				for (jj=1; 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=1; 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;
			}
		}

	// Reorders columns of a,v matrices and w vector ordered by abs(w) descending.
	// This appears to be an inefficient swap sort but check more thoroughly.
	// More efficient : Quick sort the W array and a parallel vector of column 
	// destinations then move the a and v columns efficiently.
	double *aa,*vv,ww;
	aa  = double_alloc1D(m+1);
	vv  = double_alloc1D(n+1);
	for (j=2; j<=n; j++) {
	    ww=w[j];
	    for (i=1; i<=m; i++) aa[i]=a[i][j];		// aa = col j of a
	    for (i=1; i<=n; i++) vv[i]=v[i][j];		// vv = col j of v 
	    k=j-1;
	    while ( k > 0 && -fabs(w[k]) > -fabs(ww) ) {
	        w[k+1]=w[k];
	        for (i=1; i<=m; i++) a[i][k+1]=a[i][k];	// copy col k of a to k+1
	        for (i=1; i<=n; i++) v[i][k+1]=v[i][k];	// copy col k of v to k+1
	        k--;
			}
	    w[k+1]=ww;
	    for (i=1; i<=m; i++) a[i][k+1]=aa[i];
	    for (i=1; i<=n; i++) v[i][k+1]=vv[i];
		}
	free(aa);
	free(vv);
	
	free(rv1);

	return rtn;
	}


/*
 * eigen utility - Reduction to hessenberg form by the elimination method.
 Replace with elmhes0 from qublib ( same ) 
 */
void elmhes0(double **a, int n){
	int m,j,i;
	double y,x;
	double swap_nTemp;
	
	for (m=2; m<n; m++) {
		x=0.0;
		i=m;
		for (j=m; j<=n; j++) {			// find the pivot
			y= fabs(a[j-1][m-1-1]);
			if ( y > fabs(x) ) {
				x=a[j-1][m-1-1];
				i=j;
				}             
			}

		if (i != m) {
			for (j=m-1; j<=n; j++) 
				SWAP(a[i-1][j-1],a[m-1][j-1])
			for (j=1; j<=n; j++) 
				SWAP(a[j-1][i-1],a[j-1][m-1])
			}

		if (x) {
			for (i=m+1; i<=n; i++) {
				if ((y=a[i-1][m-1-1]) != 0.0) {
					y /= x;
					a[i-1][m-1-1]=y;
					for (j=m; j<=n; j++) 
						a[i-1][j-1] -= y*a[m-1][j-1];
					for (j=1; j<=n; j++) 
						a[j-1][m-1] += y*a[j-1][i-1];
					}
				}
			}
		}
	}

void elmhes(double **a, int n){			
	int m,j,i;
	double y,x;
	double swap_nTemp;
	
	for (m=2; m<n; m++) {
		x=0.0;
		i=m;
		for (j=m; j<=n; j++) {			// find the pivot
			y= fabs(a[j][m-1]);
			if ( y > fabs(x) ) {
				x=a[j][m-1];
				i=j;
				}             
			}

		if (i != m) {
			for (j=m-1; j<=n; j++) 
				SWAP(a[i][j],a[m][j])
			for (j=1; j<=n; j++) 
				SWAP(a[j][i],a[j][m])
			}

		if (x) {
			for (i=m+1; i<=n; i++) {
				if ((y=a[i][m-1]) != 0.0) {
					y /= x;
					a[i][m-1]=y;
					for (j=m; j<=n; j++) 
						a[i][j] -= y*a[m][j];
					for (j=1; j<=n; j++) 
						a[j][m] += y*a[j][i];
					}
				}
			}
		}
	}

void hmges(int n, double* w, double** v, double* x){
	 const double TOL = 1.0e-8;
	 int	i,j;
	 double	wmx,s,*y;
	 
	 y=double_alloc1D(n+1);
	 wmx=0.0;
	 for (i=1; i<=n; i++)
		 if(fabs(w[i]) > wmx) 
			 wmx=fabs(w[i]);
		 
	 for (i=1; i<=n; i++)
		 if(fabs(w[i])<TOL*wmx) 
			 w[i]=0.0;
		 
	 /* y=V'x  */
	 for (i=1; i<=n; i++)
		 y[i]= (w[i]==0.0 ? 1.0 : 0.0 );   /* + ncall + i;*/

	 for (i=1; i<=n; i++) {
		 s=0.0;
		 for (j=1; j<=n; j++) 
			 s=s+v[i][j]*y[j];
		 x[i]=s;
		 }
	 
	 s=0.0;
	 for (i=1; i<=n; i++) 
		 s=s+x[i]*x[i];
	 if( s > 0.0 ) { 
		 s=sqrt(s);
		 for (i=1; i<=n; i++) 
			 x[i]=x[i]/s;
		 free(y);
		 }

	 }
 
// WARNING : a, v are 0 based arrays. wr and wi are 1 based arrays.
extern "C" QUBFAST_API BOOL eigen(double** a, int n, double* wr, double* wi, double** v, ostream &msgOut){
	const double TINY = 1.0e-5;

	int	i,j;
	BOOL lRet=TRUE;
	
     /* find eigenvalues */
	double **u=double_alloc2D(n+1,n+1);
	for (i=1; i<=n; i++) 
	   for (j=1; j<=n; j++) 
         u[i][j]=a[i-1][j-1];

	balanc(u,n); 

	elmhes(u,n);

	if ( ! hqr(u,n,wr,wi) )    {
	    free_2D( (char **)u );
	    return FALSE;
		}	
    free_2D( (char **)u );

    

	double **u1=double_alloc2D(n+1,n+1);
	double **v1=double_alloc2D(n+1,n+1);
	double *w1=double_alloc1D(n+1);
	double *x1=double_alloc1D(n+1);
	double **u2=double_alloc2D(2*n+1,2*n+1);
	double **v2=double_alloc2D(2*n+1,2*n+1);
	double *w2=double_alloc1D(2*n+1);
	double *x2=double_alloc1D(2*n+1);

     /* find eigenvectors */

	int k=1;
	while ( k <= n ) {
		if( fabs(wi[k]) < TINY ) {
			for (i=1; i<=n; i++) {
				for (j=1; j<=n; j++) 
					u1[i][j]=a[i-1][j-1];
				u1[i][i]=u1[i][i]-wr[k];
				}
			
			if ( svdecomp(u1,n,n,w1,v1,msgOut) ) {
				lRet=FALSE;
				break;
				}	

			hmges(n,w1,v1,x1);
			for (i=1; i<=n; i++) {
				v[i-1][k-1]=x1[i];
				}
			k++;
			} 
		else {
			for (i=1; i<=n; i++) {
				for (j=1; j<=n; j++) 
					u2[i][j]=a[i-1][j-1];
				u2[i][i]=u2[i][i]-wr[k];
				for (j=n+1; j<=2*n; j++) 
					u2[i][j]=0.0;
				u2[i][n+i]=wi[k];
				}
			for (i=n+1; i<=2*n; i++) {
				for (j=1; j<=n; j++) 
					u2[i][j]=-u2[i-n][j+n];
				for (j=n+1; j<=2*n; j++) 
					u2[i][j]=u2[i-n][j-n];
				}
			if ( svdecomp(u2,2*n,2*n,w2,v2,msgOut) ) {
				lRet=FALSE;
				break;
				}	
			
			hmges(2*n,w2,v2,x2);
			for (i=1; i<=n; i++) {
				v[i-1][k-1]=x2[i];
				v[i-1][k+1-1]=x2[n+i];
				}	
			k=k+2;
			}
		}

	free_2D( (char **)u1 );
	free_2D( (char **)v1 );
	free(w1);
	free(x1);
	free_2D( (char **)u2 );
	free_2D( (char **)v2 );
	free(w2);
	free(x2);
	return lRet;
	}	

// BALANCE A MATRIX TO PREVENT ROUND OFF ERRORS ON EIGENVALUE DETERMINATION.
void balanc(double** a, int n){
	const double RADIX=2.0;
	const double SQRDX=RADIX*RADIX;
	
	int j,i;
	double s,r,g,f,c;
	bool last=false;
	
	while( ! last ) {
		last=true;
		for( i=1; i<=n; i++ ) {
			r=c=0.0;
			for( j=1; j<=n; j++)
				if (j != i) {
					c += fabs(a[j][i]);
					r += fabs(a[i][j]);
					}
			if( c!=0.0 && r!=0.0 ) {
				g=r/RADIX;
				f=1.0;
				s=c+r;
				while (c<g) {
					f *= RADIX;
					c *= SQRDX;
					}
				g=r*RADIX;
				while (c>g) {
					f /= RADIX;
					c /= SQRDX;
					}
				if ((c+r)/f < 0.95*s) {
					last=false;
					g=1.0/f;
					for( j=1; j<=n; j++ ) {
						a[i][j] *= g;
						a[j][i] *= f;
						}
					}
				}
			}
		}
	}

//-- Find all roots of an upper Hessenberg matrix a[][] is destroyed.
//-- wr[] and wi[] are filled with the real and imaginary parts of the eigenvalues.
BOOL 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=0.,q=0.,p=0.,anorm;
	
	//-- compute matrix norm for possible use in locating single small subdiagonal element
	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;					// gets changed only by an exceptional shift.
	while (nn >= 1) {		// begin search for next eigenvalue.
		its=0;
		do  {
			// begin iteration : look for a single small subdiagonal element
			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) {			// one root found
				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)) {	// two roots found ...
					p=0.5*(y-x);
					q=p*p+w;
					z=sqrt(fabs(q));
					x += t;
					if (q >= 0.0) {	// ...a real pair
						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 {			// ...a complex pair
						wr[nn-1]=wr[nn]=x+p;
						wi[nn-1]= -(wi[nn]=z);
						}
					nn -= 2;
					} 
				else {						// no roots found, continue iteration
					if (its == 30) 
						return FALSE;
					if (its == 10 || its == 20) {	// form exceptional shift 
						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;
					// form shift and then look for 2 consecutive small subdiagonal elements
					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);	// scale to prevent overflow / underflow
						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++) {
						// double QR step on rows l to nn and m to nn
						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++) {	// row modification
								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++) {		/// column modification 
								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);
		}
	return TRUE;
	}