import java.io.*;
import javax.swing.*;

class SVD
{ 
double U[][];
double V[][];
double S[][];
double invS[][];

public static double[][] Transpose(double [][] left)
{
//transpose matrix (if A=a(i,j) Transpose(A)=a(j,i)
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[m][n];
for(i=0;i<n;i++)
  {
  for(j=0;j<m;j++)
    {
    b[j][i]=left[i][j];
    }
  }
return b;
}

public static double[][] T(double [][] left)
{return Transpose(left);}

public static double[][] I(int n)
{
//unit matrix
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=1.0;
return b;
}
//Converts vector to diagonal values of matrix
public static double[][] S(double[] left)
{
int n=left.length;
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=left[i];
return b;
}
//Converts double to diagonal values of matrix
public static double[][] S(double left,int n)
{
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=left;
return b;
}
//takes invert value of a diagonal matrix
public static double[][] invS(double[][] S)
{
int n=S.length;
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    {  if(S[i][i]==0) b[i][i]=0; 
	   else b[i][i]=1.0/S[i][i];
	}
return b;
}
//square root of a diagonal matrix
public static double[][] sqrtS(double[] left)
{
int n=left.length;
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=Math.sqrt(left[i]);
return b;
}
//multiply two matrices
public static double[][] multiply(double[][] left,double[][] right)
{
//multiplication of two matrices
int ii,jj,i,j,k;
int m1=left[0].length;
int n1=left.length;
int m2=right[0].length;
int n2=right.length;
double[][] b;
b=new double[n1][m2];
 if(m1 !=n2 )
 {
 System.out.println("inner matrix dimensions must agree"+"m1="+m1+"n2="+n2);
 return b;
 }
 for(i=0;i<n1;i++)
   {
   for(j=0;j<m2;j++)
     {
     for(k=0;k<m1;k++)
	    {b[i][j]+=left[i][k]*right[k][j];
	    }
     }
   }
return b;
//end of multiply of two matrices
}

//multiply a matrix with a double value
public static double[][] multiply(double[][] left,double right)
{
//multiplying a matrix with a constant
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[n][m];
for(i=0;i<n;i++)
  {
  for(j=0;j<m;j++)
    b[i][j]=left[i][j]*right;
  }
return b;
//end of multiplying a matrix with a constant double
}

public static double[] multiply(double[][] left,double[] right)
{
//multiplication of one matrix with one vector
int ii,jj,i,j,k;
int n1=left.length;
int m1=left[0].length;

int m2=right.length;
double[] b;
b=new double[m2];
 for(i=0;i<n1;i++)
   {
     b[i]=0;
     for(k=0;k<m1;k++)
	 {b[i]+=left[i][k]*right[k];
	 }
   }
return b;
//end of multiply of a matrix and a vector
}

//substract two matrices
public static 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 SVD(double A[][])
{initilize(A);}

public static double[] solve(double a[][],double b[])
{ SVD q=new SVD(a);
 double x[]=multiply(multiply(q.V,q.invS),multiply(SVD.T(q.U),b));
 return x;
}

public  double[] solve(double b[])
{
 double x[]=multiply(multiply(V,invS),multiply(T(U),b));
 return x;
}

public void initilize(double A[][])
{ 
  double ATA[][]=multiply(Transpose(A),A);
  double c[][];
  c=eigen_jacobi(ATA);
  int n=c[0].length;
  int n1=c.length;
  double dummy[]=new double[n1];
  for(int i=0;i<n;i++)
  {for(int k=i+1;k<n;k++)
    {if(c[0][i]<c[0][k])
      { for(int j=0;j<n1;j++)
        { dummy[j]=c[j][k];
          c[j][k]=c[j][i];
          c[j][i]=dummy[j];
        }
      }
    }
  }
  double x[]=new double[n];
  for(int i=0;i<n;i++)
  {x[i]=c[0][i];}
  System.out.println("x"+Matrix.toString(x));
  S=sqrtS(x);
  invS=invS(S);
  double xx[][]=new double[n][n];
  V=new double[n][n];
  U=new double[n][n];
  double L[]=new double[n];
  for(int i=0;i<n;i++)
  { for(int j=0;j<n;j++)
    {xx[i][j]=c[i+1][j];}
  }

  double bb[]=new double[n];  
  for(int j=0;j<n;j++)
  { for(int i=0;i<n;i++)
    {L[j]+=xx[i][j]*xx[i][j];}
    L[j]=Math.sqrt(L[j]);
  }
  for(int i=0;i<n;i++)
  { for(int j=0;j<n;j++)
    {V[i][j]=xx[i][j]/L[i];}
  }
  U=multiply(A,multiply(V,invS)); 
}

  public static void print(String s)
  {JOptionPane.showMessageDialog(null,s,
     "Singular value decomposition : ",JOptionPane.PLAIN_MESSAGE);}
  
  public static double[][] eigen_jacobi(double A[][])
  {	  
   double eps=1e-15;
   double err=1;
   int m=100;
   int n=A.length;
   double P[][]=new double[n][n];   
   double R[][]=new double[n][n];   
   int p=n,q=n;
   double max=-9.99e100;
   double fi;
   P=I(n);
   for(int i=0;i<n;i++)
   {for(int j=i+1;j<n;j++)
       {if(Math.abs(A[i][j]) > max) {max=Math.abs(A[i][j]);p=i;q=j;}}
   }  
   while(Math.abs(A[p][q]) > eps )
   {
    fi=0.5*Math.atan(2.0*A[p][q]/(A[p][p]-A[q][q]));
    R=I(n);
    R[p][p]=Math.cos(fi);
    R[q][q]=Math.cos(fi);
    R[p][q]=-Math.sin(fi);
    R[q][p]=Math.sin(fi);
    P=multiply(P,R);
    A=multiply(T(R),multiply(A,R));
    max=-9.99e100;
    for(int i=0;i<n;i++)
   {for(int j=i+1;j<n;j++){if(Math.abs(A[i][j])>max) {max=Math.abs(A[i][j]);p=i;q=j;}}}  
    }
   double lambda[]=new double[n];
   for(int i=0;i<n;i++) {lambda[i]=A[i][i];}
   double PP[][]=new double[n+1][n];
   PP[0]=lambda;
   for(int i=0;i<n;i++)
   {for(int j=0;j<n;j++)
     {PP[i+1][j]=P[i][j];}
   }
   return PP;
  }
  
  public double[][] A()
  {return multiply(U,multiply(S,Transpose(V)));}
   
  public static void main(String args[]) throws IOException
  {
  double a[][]={{3,2,-1},{2,-2,4},{-1,0.5,-1}};
  double b[]={1,-2,0};
  //Dynamic call to SVD
  SVD q=new SVD(a);
  double x[]=q.solve(b);
  String s="a=\n"+Matrix.toString(a)+"\n";
  s+="b=\n"+Matrix.toString(b)+"\n";
  s+="S:\n"+Matrix.toString(q.S)+"\n";
  s+="inv(S):\n"+Matrix.toString(q.invS)+"\n";
  s+="V:\n"+Matrix.toString(q.V)+"\n";
  s+="U:\n"+Matrix.toString(q.U)+"\n";
  s+="Reconstituted a matrix : A=U*S*VT\n"+Matrix.toString(q.A())+"\n";
  //Static call to SVD
  //double x2[]=solve(a,b);
  s+="SINGULAR VALUE DECOMPOSITION : \n"+Matrix.toString(x);
  print(s);
  }
}
