eigen.cpp

Computes the eigenvalues and eigenvectors of a real upper-Hessenberg Matrix using the QR method

					tst1 = t;
					tst2 = tst1 + 1.0/tst1;
					if (tst2 > tst1)
						continue; //go to top of for j loop
					
					for (ii = j; ii <= en; ii++){
						A[ii][na] /= t;
						A[ii][en] /= t;
					} //End for ii loop
					
				} // End for j
				
			}//End else q < 0.0
			
		}//End for en
		
		//End back substitution. Vectors of isolated roots.
		
		for (i = 0; i < N; i++){
			if ((i < low) || (i > igh)) {
				for (j = i; j < N; j++){
					B[i][j] = A[i][j];
				}
			}//End if i
		}//End for i
		
		// Multiply by transformation matrix to give vectors of original full matrix.
		
		//Step from (N - 1) to low in steps of -1.
		
		for (i = (N - 1); i >= low; i--){
			
			m = ((i < igh) ? i : igh);
			
			for (ii = low; ii <= igh; ii++){
				zz = 0.0;
				for (jj = low; jj <= m; jj++)
					zz += B[ii][jj]*A[jj][i];
				B[ii][i] = zz;
			}//End for ii
		}//End of for i loop
		return;
	} //End of function hqr2
	
	void Eigen::norVec(double **Z, double *wi, int N){// Normalize the eigenvectors
		
		double s;
		int i,j = 0;
		
		do {
			s = 0.0;
			if (wi[j]  != 0){
				for (i = 0; i < N; i++)
					s += Z[i][j]*Z[i][j] + Z[i][j + 1]*Z[i][j + 1];
				s = sqrt(s);
				for (i = 0; i < N; i++){
					Z[i][j] /= s;
					Z[i][j + 1] /= s;
				} // End for i
				j++;
			} // End if (wi[j]  != 0)
			else { //else wi[j]  == 0
				for (i = 0; i < N; i++)
					s += Z[i][j]*Z[i][j];
				s = sqrt(s);
				for (i = 0; i < N; i++)
					Z[i][j] /= s;
			} // End else wi[j]  == 0
			j++;
		} while (j < N);
		return;
	} // End norVec
	
	void Eigen::balance (double **A, int N, int *trace, double *scale, int &low, int &igh){
		int i,j;
		/* Balance the matrix to improve accuracy of eigenvalues. Introduces no rounding errors, since 
		it scales A by powers of the radix. */
		
		int radix = 2;               //Base of machine floating point representation.
		double c, f, g, r, s, b2 = radix*radix;
		int k = 0, M = N - 1, noconv;
		
		//Search through rows, isolating eigenvalues and pushing them down.
		
		noconv = M;
		
		while (noconv >= 0){
			r = 0;
			
			for (j = 0; j <= M; j++) {
				if (j == noconv) continue;
				if (A[noconv][j] != 0.0){
					r = 1;
					break;
				}
			} //End for j
			
			if (r == 0){
				scale[M] = noconv;
				
				if (noconv != M){
					for (i = 0; i <= M; i++){
						f = A[i][noconv];
						A[i][noconv] = A[i][M];
						A[i][M] = f;
					}//End for i
					for (i = 0; i < N; i++){
						f = A[noconv][i];
						A[noconv][i] = A[M][i];
						A[M][i] = f;
					}//End for i
				}//End if (noconv != M)
				
				if (M == 0)
					break;  //break out of while loop
				
				noconv = --M;
				
			}//End if (r == 0)
			
			else //else (r != 0)
				noconv--;
			
		}//End while (noconv >= 0)
		
		if (M > 0) {  //Search through columns, isolating eigenvalues and pushing them left.
			
			noconv = 0;
			
			while (noconv <= M){  
				c = 0;
				
				for (i = k; i <= M; i++){
					if (i == noconv) continue;
					if (A[i][noconv] != 0.0){
						c = 1;
						break;
					}
				}//End for i
				
				if (c == 0){
					scale[k] = noconv;
					
					if (noconv != k){
						for (i = 0; i <= M; i++){
							f = A[i][noconv];
							A[i][noconv] = A[i][k];
							A[i][k] = f;
						}//End for i
						for (i = k; i < N; i++){
							f = A[noconv][i];
							A[noconv][i] = A[k][i];
							A[k][i] = f;
						}//End for i
						
					}//End if (noconv != k)
					
					noconv = ++k;
				}//End if (c == 0)
				
				else  //else (c != 0)
					noconv++;
				
			}//End while noconv
			
			//Balance the submatrix in rows k through M.
			
			for (i = k; i <= M; i++)
				scale[i] = 1.0;
			
			//Iterative loop for norm reduction
			do {
				noconv = 0;
				for (i = k; i <= M; i++) {
					c = r = 0.0;
					for (j = k; j <= M; j++){
						if (j == i) continue;
						c += fabs(A[j][i]);
						r += fabs(A[i][j]);
					} // End for j
					if ((c == 0.0) || (r == 0.0)) continue;   //guard against zero c or r due to underflow
					g = r/radix;
					f = 1.0;
					s = c + r;
					while (c < g) {
						f *= radix;
						c *= b2;
					} // End while (c < g)
					g = r*radix;
					while (c >= g) {
						f /= radix;
						c /= b2;
					} // End while (c >= g)
					
					//Now balance
					if ((c + r)/f < 0.95*s) {
						g = 1.0/f;
						scale[i] *= f;
						noconv = 1;
						for (j = k; j < N; j++)
							A[i][j] *= g;
						for (j = 0; j <= M; j++)
							A[j][i] *= f;
					} //End if ((c + r)/f < 0.95*s)
				} // End for i
			} while (noconv);  // End of do-while loop.
			
		} //End if (M > 0)
		
		low = k;
		igh = M;
		
		//End of balanc
	}
	
	void Eigen::upperHessenberg(double **A, int N, int *trace, double *scale, int low, int igh){
		int i, j;
		/* Now reduce the real general Matrix to upper-Hessenberg form using stabilized elementary 
		similarity transformations. */
		
		for (i = (low + 1); i < igh; i++){
			int k = i;
			double c = 0.0;
			
			for (j = i; j <= igh; j++){
				if (fabs(A[j][i - 1]) > fabs(c)){
					c = A[j][i - 1];
					k = j;
				}//End if
			}//End for j
			
			trace[i] = k;
			
			if (k != i){
				for (j = (i - 1); j < N; j++){
					double r = A[k][j];
					A[k][j] = A[i][j];
					A[i][j] = r;
				}//End for j
				
				for (j = 0; j <= igh; j++){
					double r = A[j][k];
					A[j][k] = A[j][i];
					A[j][i] = r;
				}//End for j
			}//End if (k != i)
			
			if (c != 0.0){
				for (int m = (i + 1); m <= igh; m++){
					double r = A[m][i - 1];
					
					if (r != 0.0){
						r = A[m][i - 1] = r/c;
						for (j = i; j < N; j++)
							A[m][j] += -(r*A[i][j]);
						for (j = 0; j <= igh; j++)
							A[j][i] += r*A[j][m];
					}//End if (r != 0.0)
				}//End for m
			}//End if (c != 0)
		}  //End for i.
		
	}
	
	void Eigen::Compute(double **A, double **B, double *wr, double *wi, int N, bool sort) {
		int i, j;
		MemUtil u;

		double *scale = new double[0];    //Contains information about transformations.
		u.GetMemory(scale, N);
		int *trace = new int[0];    //Records row and column interchanges
		u.GetMemory(trace, N);
		int ierr = -1, igh, low;
		
		balance(A, N, trace, scale, low, igh);
		upperHessenberg(A, N, trace, scale, low, igh);
		
		/* Accumulate the stabilized elementary similarity transformations used in the reduction 
		of A to upper-Hessenberg form. Introduces no rounding errors since it only transfers the 
		multipliers used in the reduction process into the eigenvector matrix. */
		
		for (i = 0; i < N; i++){ //Initialize B to the identity Matrix.
			for (j = 0; j < N; j++)
				B[i][j] = 0.0;
			B[i][i] = 1.0;
		} //End for i
		
		for (i = (igh - 1); i >= (low + 1); i--){
			int k = trace[i];
			for (j = (i + 1); j <= igh; j++)
				B[j][i] = A[j][i - 1];
			
			if (i == k)
				continue;
			
			for (j = i; j <= igh; j++){
				B[i][j] = B[k][j];
				B[k][j] = 0.0;
			} //End for j
			
			B[k][i] = 1.0; 
			
		} // End for i
		
		hqr2(A, N, B, low, igh, wr, wi, ierr);
		
		if (ierr == -1){
			
			if (low != igh){
				for (i = low; i <= igh; i++){
					double s = scale[i];
					for (j = 0; j < N; j++)
						B[i][j] *= s;
				}//End for i
			}//End if (low != igh)
			
			for (i = (low - 1); i >= 0; i--){
				int k = round(scale[i]);
				if (k != i){
					for (j = 0; j < N; j++){
						u.swap(B[i][j], B[k][j]);
					}//End for j
				}//End if k != i
			}//End for i
			
			for (i = (igh + 1); i < N; i++){
				int k = round(scale[i]);
				if (k != i){
					for (j = 0; j < N; j++){
						u.swap(B[i][j], B[k][j]);
					}//End for j
				}//End if k != i
			}//End for i
			
			norVec(B, wi, N);  //Normalize the eigenvectors
			
		}//End if ierr = -1
		else printf("Error: %d\n", ierr);
		
		//Sorting of Eigen Values and Vectors
		if(sort){
			for(i=0;i<N;i++){
				int maxI = i;
				for(j=i+1;j<N;j++){
					if(wr[j]>wr[maxI]) maxI=j;
				}
				if(maxI!=i){
					//swap eigen values
					u.swap(wr[i], wr[maxI]);
					u.swap(wi[i], wi[maxI]);
					//swap eigen vectors
					for(j=0;j<N;j++){
						u.swap(B[j][i], B[j][maxI]);
					}
				}
			}
		}
		for(i=0;i<N;i++){
			double sum=0.0;
			for(j=0;j<N;j++){ 
				sum+=B[j][i];
			}
			if(sum<0){
				for(j=0;j<N;j++){ 
					B[j][i] = -1*B[j][i];
				}
			}
		}
		delete [] scale;
		delete [] trace;
		return;
	}  
	
	/*
	The four by four matrix H[i,j] = i^j + j^i  is symmetric, so the
	eigenvalues are real.  The output of this program should be:
	
	  Eigenvalue  0: 0.151995 + i 0
	  Eigenvector 0: (0.610828, -0.72048, 0.324065, -0.0527271)
	  Eigenvalue  1: 1.80498 + i 0
	  Eigenvector 1: (-0.7614, -0.42123, 0.481884, -0.103072)
	  Eigenvalue  2: 17.6352 + i 0
	  Eigenvector 2: (-0.216885, -0.547084, -0.764881, 0.26195)
	  Eigenvalue  3: 556.408 + i 0
	  Eigenvector 3: (0.0110019, 0.064609, 0.278795, 0.958112)
	  */