// SCO5E Broyden secant method with Sherman-Morrison formula to solve inverse matrix
// Solution of non-linear system of equations.

// Referans : Numerical Analysis, 8th edition Richard L. Burden, J. Douglas Faires 3
// page 617

import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Tanımı gereken fonksiyon ================



class fa extends fi_xi
{
public double[] func(double x[])
{	
//function to be solved.
double ff[]=new double[3];
ff[0]=3.0*x[0]-Math.cos(x[1]*x[2])-0.5;
ff[1]=x[0]*x[0]-81.0*(x[1]+0.1)*(x[1]+0.1)+Math.sin(x[2])+1.06;
ff[2]=Math.exp(-x[0]*x[1])+20.0*x[2]+(10.0*Math.PI-3.0)/3.0;
return ff;
} 
}

class fb extends fij_xi
{
fi_xi f;	

public fb(fi_xi fi)
{f=fi;}

public double[][] func(double x[])
{ 
  int n=x.length;
  double a[][]=new double[n][n];
  for(int i=0;i<n;i++)
  {
    for(int j=0;j<n;j++)
    {
        a[i][j]=derivative(x,i,j);
    }
  }
return a;
} 

public double derivative(double x[],int denklem_ref,int x_ref)
{
// df(x)/dx[x_ref]  x_ref=0...n 
double h0=0.0256;
int i,m;
int n=7;
double f1[];
f1=new double[x.length];
double f2[];
f2=new double[x.length];
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
f1=f.func(x1);
f2=f.func(x2);
T[i][0]=( f1[denklem_ref] - f2[denklem_ref])/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

}

public class NA30A
{
// SCO5E Sherman-Morrison formülü destekli broyden sekant metoduyla
//denklem sistemlerinin köklerini (çözümünü) bulmak

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

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

public static double  VT_X(double [] left,double [] right)
{
//multiplys a vector transpose with a vector	
int n=left.length;
double tot=0;
for(int i=0;i<n;i++)
  {
  tot+=left[i]*right[i];
  }
return tot;  
}

public static double[] add(double[] left,double[] right)
{
//addition of two vectors
int n1=left.length;
int n2=right.length;
int nMax;
int i;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[];
b=new double[nMax];
for(i=0;i<n1;i++)
  {
    b[i]=b[i]+left[i];
  }
for(i=0;i<n2;i++)
  {
    b[i]=b[i]+right[i];
  }
return b;
//end of vector addition method
}

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 static double[][] add(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 addition method
}

public static double[] substract(double[] left,double[] right)
{
//addition of two vectors
int n1=left.length;
int n2=right.length;
int nMax;
int i;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[];
b=new double[nMax];
for(i=0;i<n1;i++)
  {
    b[i]=b[i]+left[i];
  }
for(i=0;i<n2;i++)
  {
    b[i]=b[i]-right[i];
  }
return b;
//end of vector addition method
}

public static double[][] inv(double[][] a)
{
// INVERSION OF A MATRIX
// inversion by using gaussian elimination
// with full pivoting
int n=a.length;
int m=a[0].length;
double b[][];
b=new double[n][n];
int indxc[];
int indxr[];
double ipiv[];
indxc=new int[n];
indxr=new int[n];
ipiv=new double[n];
int i,j,k,l,ll,ii,jj;
int icol=0;
int irow=0;
double big,dum,pivinv,temp;
if(n!=m)
{
   System.out.println("Matrix must be square ");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
	  return b;
}
for(i=0;i<n;i++)
   for(j=0;j<n;j++)
       b[i][j]=a[i][j];
for(i=0;i<n;i++)
{
  big=0.0;
  for(j=0;j<n;j++)
  {
    if(ipiv[j] != 1)
      for(k=0;k<n;k++)
      {
	if(ipiv[k] == 0)
	{
	  if(Math.abs(b[j][k]) >= big)
	  {
	    big=Math.abs(b[j][k]);
	    irow=j;
	    icol=k;
	  }
	}
	else if(ipiv[k] > 1 )
	  {
	  System.out.println("error : inverse of the matrix : singular matrix-1");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
	  return b;
	  }
      }
  }
  ++ ipiv[icol];
  if(irow != icol)
    for(l=0;l<n;l++)
    {
    temp=b[irow][l];
    b[irow][l]=b[icol][l];
    b[icol][l]=temp;
    }
  indxr[i]=irow;
  indxc[i]=icol;
  if(b[icol][icol] == 0.0)
  {
    System.out.println("error : inverse of the matrix : singular matrix-2");
	  for(ii=0;ii<n;ii++)
	    for(jj=0;jj<n;jj++)
		b[ii][jj]=0.0;
    return b;
  }
  pivinv=1.0/b[icol][icol];
  b[icol][icol]=1.0;
  for(l=0;l<n;l++) b[icol][l] *=pivinv;
  for(ll=0;ll<n;ll++)
  if(ll != icol)
  {
    dum=b[ll][icol];
    b[ll][icol]=0.0;
    for(l=0;l<n;l++) b[ll][l]-= b[icol][l]*dum;
  }
}
for(l=n-1;l>=0;l--)
{
  if(indxr[l] != indxc[l])
     for(k=0;k<n;k++)
     {
       temp=b[k][indxc[l]];
       b[k][indxc[l]]=b[k][indxr[l]];
       b[k][indxr[l]]=temp;
     }
}
return b;
}
public static double norm(double v[])
{
// vector norm
double total=0;;
for(int i=0;i<v.length;i++)
  {total+=v[i]*v[i];}
return Math.sqrt(total);
}

public static double[][] multiply(double[] left,double[] right)
{
//multiplication of one vector with one matrix
int ii,jj,i,j,k;
int m2=left.length;
int n2=right.length;
int m1=left.length;
double[][] b;
b=new double[m2][n2];
 for(i=0;i<m2;i++)
   {
     for(k=0;k<n2;k++)
	   b[i][k]=left[i]*right[k];
   }
return b;
//end of multiply of a vector and a matrix
}
  
public static double[] broyden(fi_xi f,double x[])
{
//====================================================	
// Roots of non-linear system of equations
// Broyden secant method
fb J=new fb(f);
int i;
int nmax=400;
double tolerance=1.0e-15;
int n=x.length;
double b[]=new double[n];
double v[]=new double[n];
double s[]=new double[n];
double w[]=new double[n];
double y[]=new double[n];
double z[]=new double[n];
double p;
double ut[]=new double[n];
double x1[]=new double[n];
double A[][]=new double[n][n];
x1=x;	
A=J.func(x1);
v=f.func(x1);
A=inv(A);
s=multiply(-1.0,multiply(A,v));
x=add(x,s);
int k=1;
while(k<nmax)
{
w=v;
v=f.func(x);
y=substract(v,w);
z=multiply(-1.0,multiply(A,y));
p=-VT_X(s,z);
ut=multiply(s,A);  //*********
A=add(A,multiply((1.0/p),multiply(add(s,z),ut)));
s=multiply(-1.0,multiply(A,v));
x=add(x,s);
if(norm(s)<tolerance) break;
k++;	
}
return x;
}

public static double[] enterdata()
{
    String s=JOptionPane.showInputDialog("enter first gueses for roots with a sapce in between: ");        
    StringTokenizer token=new StringTokenizer(s);
    int n=token.countTokens()-1;
    int m=n+1;
    double a[]=new double[m];
    int j=0;           
    while(token.hasMoreTokens())
    {
    Double ax=new Double(token.nextToken());
    a[j++]=ax.doubleValue();
    }
    return a; 
}

     public static void main (String args[])  
     {
	 fa f =new fa();
     double [] x0=enterdata();
     double [] r= broyden(f,x0);
     String s=" roots  : \n"+Matrix.toStringT(r);
     String s1="Broyden secant method with Sherman-Morrison iterative inverse matrix formula : ";
     JOptionPane.showMessageDialog(null,s,s1,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}
