
import java.io.*;
import javax.swing.*;
class residual
{
  public static void pivot(double a[][],int o[],double s[],int k)
  {// Gauss(Dolittle) LU decomposition with partial pivoting
    int n=a.length;
    int dummy1=0;
    double buyuk;
    double dummy=0;
    //pivotlama
	int p=k;
	buyuk=Math.abs(a[o[k]][k]);
	for(int ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[o[ii]][k]/s[o[ii]]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	dummy1=o[p];
	o[p]=o[k];
	o[k]=dummy1;	   
  }
  
  public static int[] sirali_indeks(int n)
  {   // sıralı rakamları taşıyan bir referans vektörü oluşturur 
	  int o[]=new int[n];
	  for(int i=0;i<n;i++) {o[i]=i;}
	  return o;
  }	  
  
  public static double[][] gaussLU(double ai[][],int o[])
  { //kısmi pivotlu gauss eleme yöntemi
  int n=ai.length;	    
  double a[][]=new double[n][n];
  double s[]=new double[n];
  double carpan=0;
  double toplam=0;
  double absaij;

  //gauss eleme
  int i,j,k,p,ii,jj;
  //System.out.println("Orijinal Matris :\n"+Matrix.toString(ai));
  for(i=0;i<n;i++)
  {for(j=0;j<n;j++){a[i][j]=ai[i][j];}}
  for(i=0;i<n;i++)
  { o[i]=i;s[i]=Math.abs(a[i][0]);
    for(j=1;j<n;j++)
      { absaij=Math.abs(a[i][j]);
        if(absaij>s[i]) {s[i]=absaij;}
      }
  }
  for(k=0;k<(n-1);k++)
  { pivot(a,o,s,k);
	//pivot referansını değiştir		
	//System.out.println("pivotlu Matris : k="+k+"\n"+Matrix.toString(a));
	//gauss elemeyi çözme  
    for(i=k+1;i<n;i++)
    {  carpan=a[o[i]][k]/a[o[k]][k];
       a[o[i]][k]=carpan;
       for(j=k+1;j<n;j++)
       {   a[o[i]][j]-=carpan*a[o[k]][j]; }
	}
  }
  return a;
  }
  
  public static double[] yerine_koyma(double ai[][],double bi[])
  {   int n=ai.length;
	  return yerine_koyma(ai,bi,sirali_indeks(n));
  }
  
  public static double[][] yerine_koyma(double ai[][],double bi[][])
  {   int n=ai.length;
	  return yerine_koyma(ai,bi,sirali_indeks(n));
  }  
  
  public static double[] yerine_koyma(double ai[][],double bi[],int o[])
  { //gauss LU çözümüyle LU indirgenmesi yapılmış Matris kullanarak lineer 
    //denklem sistemi kökleri
    int n=ai.length;
    double toplam;
    double x[]=new double[n];
     double a[][]=new double[n][n];
    double b[]=new double[n];
    int i,j;
  for(i=0;i<n;i++)
      {for(j=0;j<n;j++){a[i][j]=ai[i][j];};b[i]=bi[i];}   
  //System.out.println("a=\n "+Matrix.toString(a));
  //System.out.println("b=\n "+Matrix.toStringT(b));             
  for(i=1;i<n;i++)
      {toplam=b[o[i]];for(j=0;j<i;j++){toplam-=a[o[i]][j]*b[o[j]];};b[o[i]]=toplam;}
   //System.out.println("d=\n "+Matrix.toStringT(b));      
  x[n-1]=b[o[n-1]]/a[o[n-1]][n-1];
  for(i=(n-2);i>=0;i--)
      {toplam=0;      
       for(j=i+1;j<n;j++)      
        {toplam+=a[o[i]][j]*x[j];};
        x[i]=(b[o[i]]-toplam)/a[o[i]][i];
      }
  //System.out.println("d=\n "+Matrix.toStringT(x));        
  return x;              
  }  
  
  public static double[][] yerine_koyma(double ai[][],double bi[][],int o[])
  {
    int n=ai.length;
    int m=bi[0].length;
    double toplam;
    double x[][]=new double[n][m];
    double a[][]=new double[n][n];
    double b[][]=new double[n][m];
    int i,j,k;
  for(i=0;i<n;i++)
      {for(j=0;j<n;j++){a[i][j]=ai[i][j];};for(k=0;k<m;k++) {b[i][k]=bi[i][k];}}   
  //System.out.println("a=\n "+Matrix.toString(a));
  //System.out.println("b=\n "+Matrix.toStringT(b));             
  for(k=0;k<m;k++)
  {
  for(i=1;i<n;i++)
      {toplam=b[o[i]][k];for(j=0;j<i;j++){toplam-=a[o[i]][j]*b[o[j]][k];};b[o[i]][k]=toplam;}
   //System.out.println("d=\n "+Matrix.toStringT(b));      
  x[n-1][k]=b[o[n-1]][k]/a[o[n-1]][n-1];
  for(i=(n-2);i>=0;i--)
      {toplam=0;      
       for(j=i+1;j<n;j++)      
        {toplam+=a[o[i]][j]*x[j][k];};
        x[i][k]=(b[o[i]][k]-toplam)/a[o[i]][i];
      }
  }//kin sonu   
  //System.out.println("d=\n "+Matrix.toStringT(x));        
  return x;          
  }
  
  public static double [] AXB(double a[][],double b[])
  {
  //denklem sistemi çözümü b vektörü için	  
  int n=a.length;	  
  int q[]=new int[n];
  double c[][]=gaussLU(a,q);
  double s[]=yerine_koyma(c,b,q);
  return s;
  }
  
  public static double [][] AXB(double a[][],double b[][])
  {//denklem sistemi çözümü b matrisi için	  
  int n=a.length;	  
  int q[]=new int[n];
  double c[][]=gaussLU(a,q);
  double s[][]=yerine_koyma(c,b,q);
  return s;
  }

  
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
}

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 substraction method
}
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 error(double q[])
 { double err=0;
   for(int i=0;i<q.length;i++)
   {err+=q[i]*q[i];}
   return Math.sqrt(err);
 }
 
  public static double[] residual_correction(double a[][],double b[],double eps)
  {
  double eps1;
  System.out.println("a\n"+Matrix.toString(a));
  int n=a.length;	  
  int q[]=new int[n];
  double c[][]=gaussLU(a,q);
  double x[]=yerine_koyma(c,b,q);
  System.out.println("x\n"+Matrix.toString(x));
  double e[]=new double[x.length];
  double r[]=new double[x.length];
  do
  {	  
  r= substract(b,multiply(a,x));
  e=yerine_koyma(c,r,q);
  System.out.println("r\n"+Matrix.toString(r));
  System.out.println("e\n"+Matrix.toString(e));
  x=add(x,e);
  System.out.println("x\n"+Matrix.toString(x));
  eps1=error(e);
  } while(eps1>eps);
  return x;  
  }
    
  public static void main(String args[]) throws IOException
  {
  double eps=1e-10;
  double a[][]=Mathd.Hilbert(4);
  
  double b[]={1,0,0,0};
  double x[]=residual_correction(a,b,eps);
  Text.print(x,"Solution of Gauss(Doolittle) decomposition and elimination with residual correction method");
  }
}


