/**************************************************************
 *		@author:		Manfred Vogel
 *		@date:			09.06.2005
 *		@filename:	    Wielandt
 **************************************************************
 */

/**
 * Hinweis : Mit der Verwendung der Indices i-1 etc. in den Matrizen
 *           können sowohl Algorithmus wie auch Initialiserung der
 *           Matrix A direkt übernommen werden !
 */
public class Wielandt {

	static int n = 10;
	static int wielang = 150;
	static double[][] A = new double[n][n];
	static double[] eigenWerte  = new double[n];
 
	public static void main(String[] args) {
		double[][] B;
		double[] x,y,a;
		double[] eigenVec1;
		double mu=0.0;
		double c,norm00;
		int mi;
		
		for (int i = 1; i <= n; i++)
		   for (int k = 1; k <= n; k++) 
			  A[i-1][k-1] = 1.0 + (i+k)*Math.cos(i*k);
 		
		printMatrix(A);

        for (int m = 1; m <= n-1; m++)  {
           x   = new double[n+1-m];
           x[0]=1.0;
		   y   = new double[x.length];
           for (int j = 1; j <= wielang; j++){
		      y = matrixMult( A , x);
              mu= skalar( x , y);
			  x = vectorNormalize( y );
           }
           
           eigenVec1 = new double[x.length];
           for(int i= 0; i<x.length; i++) eigenVec1[i]=x[i];
           mi=0;
		   norm00 = infinityNorm( x );
           for (int i = 1; i <= n-m+1; i++)
               if ( Math.abs(eigenVec1[i-1]) == norm00) mi=i-1;

           a = new double[x.length];
           for(int i=1;i<=n+1-m;i++) a[i-1] = A[mi][i-1]/(mu*eigenVec1[mi]);

           B = new double[x.length][x.length];
           for(int i=0; i<x.length; i++)
              for(int k=0; k<a.length; k++)
                 B[i][k]=A[i][k]-mu*x[i]*a[k];

           for(int i=1;i<=n-m+1;i++) {
              c = B[i-1][0];
              B[i-1][0] = B[i-1][mi];
              B[i-1][mi]= c;
              c = B[0][i-1];
              B[0][i-1] = B[mi][i-1];
              B[mi][i-1] = c;
           }
		   
           A = new double[A.length-1][A[0].length-1];
           for(int i=0; i<A.length; i++)
               for(int k=0; k<A[i].length; k++)
                  A[i][k]=B[i+1][k+1];
		   
           eigenWerte[m-1] = mu; 
        }
        eigenWerte[n-1] = A[0][0]; 
		
        printVector("Vektor der Eigenwerte :" , eigenWerte);
    }
	
//...... Methoden ...................................................
	
	static void printMatrix(double[][] matrix) {
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix.length; j++) {	
				InOut.print(matrix[i][j], 1, 3);
				System.out.print("  ");
			}
			System.out.print("\n");
		}
		System.out.print("\n");
		System.out.print("\n");
	}

	static void printVector(double[] vector) {
		for (int i = 0; i < vector.length; i++) {
			InOut.print(vector[i], 1, 4);
			System.out.print("\n");
		}
		System.out.print("\n");
		System.out.print("\n");
	}

	static void printVector(String was ,double[] vector) {
		System.out.println(was);
		for (int i = 0; i < vector.length; i++) {
			System.out.print("Vektorelement "+i+" :  ");
			InOut.print(vector[i], 1, 4);
			System.out.print("\n");
		}
		System.out.print("\n");
		System.out.print("\n");
	}

	private static double[] matrixMult(double[][] A, double[] x) {
		if (A[0].length != x.length ) {
		   System.out.println("dimension mismatch in maxtrixMult");
		   printMatrix(A);
		   printVector(x);
		   System.exit(0);
		}
		double[] y = new double[x.length];
		for( int i = 0; i<x.length; i++) 
		   for(int k= 0; k<x.length; k++) 
			  y[i] += A[i][k]*x[k];
		
		return y;
	}

	private static double[][] matrixSubtraction(double[][] A, double[][] B) {
		if (A.length != B.length ) {
		   System.out.println("dimension mismatch in maxtrixSubtraction");
		   printMatrix(A);
		   printMatrix(B);
		   System.exit(0);
		}
		double[][] C = new double[A.length][A[0].length];
		for( int i = 0; i<A.length; i++) 
			for(int k= 0; k<A[i].length; k++) C[i][k] = A[i][k]-B[i][k];
		
		return C;
	}

	private static double skalar(double[] x, double[] y) {
		double s=0.0;
			for(int i= 0; i<x.length; i++) s += x[i]*y[i];
		return s;
	}


	private static double[] vectorNormalize(double[] x) {
		double s=0.0;
			for(int i= 0; i<x.length; i++) s += x[i]*x[i];
			s = Math.sqrt(s);
		    for(int i= 0; i<x.length; i++) x[i] /= s;
		return x;
	}

	private static double infinityNorm(double[] x) {
		double max=0.0;
		    for(int i= 0; i<x.length; i++) max=Math.max( max, Math.abs(x[i]));
		return max;
	}


}