public class Inverse_power_method
{
//To approximate an eigenvalue and an associated eigenvector 
//of the n × n matrix A given a nonzero vector x:	

//  Numerical Analysis -- Richard L. Burden -- J.Douglas Faires
//  Page 585	9th edition 
double A[][];
double x[];
int N;
 
     public Inverse_power_method(double AA[][],double xx[],int NN)
     {
	 A=AA;
	 x=xx;
	 N=NN;    
	 }
     
	 
	 public static double[][] multiply(double left,double[][] right)
     {
//multiplying a matrix with a constant
     int i,j;
     int n=right.length;
     int m=right[0].length;
     double b[][];
     b=new double[n][m];
     for(i=0;i<n;i++)
     {
        for(j=0;j<m;j++)
        b[i][j]=right[i][j]*left;
     }
      return b;
//end of multiplying a matrix with a constant double
      } 
	 
	 public double multiply(double[] left,double right[])
	 {
	 int n=left.length;
	 double s=0;
	 for(int i=0;i<n;i++)
	 {s+=left[i]*right[i];}	 
	 return s;	 
	 }

     public  double[] multiply(double[][] left,double[] right)
     {
         //multiplication of one matrix with one vector
          int ii,jj,i,j,k;
          int m1=left[0].length;
          int n1=left.length;
          int m2=right.length;
          double[] b;
          b=new double[m2];
          if(n1 != m2)
          {
             System.out.println("inner matrix dimensions must agree");
             for(ii=0;ii<n1;ii++)
             {
               b[ii]=0;
             }
             return b;
          }
             for(i=0;i<m1;i++)
             {
                b[i]=0;
                for(k=0;k<n1;k++)
	            b[i]+=left[i][k]*right[k];
             }
             return b;
             //end of multiply of a matrix and a vector
      }

      public  double norm(double v[],int p)
      {
        //p vector norm
        if(p==0) p=1;
        double total=0;
        double x,y;
        for(int i=0;i<v.length;i++)
        {
         x=Math.abs(v[i]);
         total+=Math.pow(x,p);
        }
        return Math.pow(total,(1.0/p));
      }

      public  double[] pivotlugauss(double a[][],double b[])
      { 
        int n=b.length;
        double x[]=new double[n];
        double carpan=0;
        double toplam=0;
        double buyuk;
        double dummy=0;
     
        int i,j,k,p,ii,jj;
       
        for(k=0;k<(n-1);k++)
          { 
            p=k;
            buyuk=Math.abs(a[k][k]);
              for(ii=k+1;ii<n;ii++)
              { dummy=Math.abs(a[ii][k]);
               if(dummy > buyuk) {buyuk=dummy;p=ii;}
              }
              if(p!=k)
              { 
	             for(jj=k;jj<n;jj++)
                 { dummy=a[p][jj];
                   a[p][jj]=a[k][jj];
                   a[k][jj]=dummy;
                 }
              dummy=b[p];
              b[p]=b[k];
              b[k]=dummy;
              }
              
              for(i=k+1;i<n;i++)
              { carpan=a[i][k]/a[k][k];
                a[i][k]=0;
               for(j=k+1;j<n;j++)
               { a[i][j]-=carpan*a[k][j]; }
                 b[i] =b[i] -carpan*b[k];
              }
           }
            
            x[n-1]=b[n-1]/a[n-1][n-1];
            for(i=n-2;i>=0;i--)
            {
             toplam=0;
             for(j=i+1;j<n;j++)
             { toplam+=a[i][j]*x[j];}
             x[i]=(b[i]-toplam)/a[i][i];
            }
            // System.out.println(s);
             return x;
       }

       public double[][] substract(double[][] left,double[][] right)
       {
           //addition of two matrices
           int n1=left.length;
           int m1=left[0].length;
           int n2=right.length;
           int m2=right[0].length;
           int nMax,mMax;
           int i,j;
           if(m1>=m2) mMax=m1;
           else       mMax=m2;
           if(n1>=n2) nMax=n1;
           else       nMax=n2;
           double b[][];
           b=new double[nMax][mMax];

           for(i=0;i<n1;i++)
           {
              for(j=0;j<m1;j++)
              {
               b[i][j]=b[i][j]+left[i][j];
              }
           }
  

           for(i=0;i<n2;i++)
           {
               for(j=0;j<m2;j++)
               {
               b[i][j]=b[i][j]-right[i][j];
               }
           }
           return b;
//end of matrix substraction method
       }

       public static double norminf(double v[])
       {
     //infinite vector norm
       double x;
       double max=0;
       for(int i=0;i<v.length;i++)
        {
        x=Math.abs(v[i]);
        if(x>max) max=x;
        }
        return max;
       }
       
       public double[][] eye()
       {
	   int n=A.length;
	   double[][] I=new double[n][n];
	   for(int i=0;i<n;i++)
	   {   
	    I[i][i]=1; 	      
       }    
	   return I;    
	   }
       
       
       public void eigenvector()
       {
	   double I[][]=eye();    
	   int k=1;    
	   double q=0;
	   int n=A.length;
	   double Ax[]=new double[n];    
	   double qI[][]=new double[n][n];
	   double AqI[][]=new double[n][n];
	   Ax=multiply(A,x);
	   double xax=multiply(x,Ax);    
	   double xx=multiply(x,x);
	   q=xax/xx;
	   for(int i=0;i<n;i++)
	   {x[i]=x[i]/x[0];}
	   qI=multiply(q,I);
	   AqI=substract(A,qI);
	   double mu=0;
	   double[]y=new double[n];
	   double yp=0;
	   
	   while(k<N)
	   {	   
	       y=pivotlugauss(AqI,x);
	       	   
		   mu=-yp;
	       yp=norminf(y);  	   
		   for(int i=0;i<n;i++)
		   {x[i]=y[i]/-yp;}
		   k++;
      }
	     mu=(1/mu)+q;     
	        System.out.println("Maximum eigenvalue of a matrix = "+mu);
	        System.out.println("Eigenvectors of maximum eigenvalue");
	        for(int i=0;i<n;i++)
	        {System.out.println(x[i]);}  
	       
	   }

       public static void main(String args[])
       {
	   double A[][]={{-4,14,0},{-5,13,0},{-1,0,2}};
	   double x[]={1,1,1};
	   int N=100;  
	   Inverse_power_method ipm=new Inverse_power_method(A,x,N);
	      ipm.eigenvector();
	       
	       
	       
	   }
       
       
       

}