/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package geomExtension;

/**
 *
 * @author jyounishi
 */
import java.text.*;
//import java.awt.geom.*;
    
public class Matrix {
    static double eps=1e-12;
    public Matrix(){}
    
    public static double[] GaussianElimination(double[][] A0, double[] b0) {
        
        int N=A0.length;
        double[][] A=new double[N][N];
        double[] b=new double[N];
        for(int i=0;i<N;i++){
            for(int j=0;j<N;j++) A[i][j]=A0[i][j];
            b[i]=b0[i];
        }
        
        for (int i=0; i<N; i++) {
            // find pivot row and swap
            int max=i;
            for (int ip=i+1;ip<N;ip++) {
                if (Math.abs(A[ip][i])>Math.abs(A[max][i])) {
                    max=ip;
                }
            }
            double[] temp=A[i]; 
            A[i]=A[max]; 
            A[max]=temp;
            double  btemp=b[i]; 
            b[i]=b[max]; 
            b[max]=btemp;
            
            // singular or nearly singular
            if (Math.abs(A[i][i]) <= eps) {
                System.out.println("** aussianElimination : Matrix is singular or nearly singular");
                throw new RuntimeException("Matrix is singular or nearly singular");
            }

            for (int ip=i+1;ip<N;ip++) {
                double alpha=A[ip][i]/A[i][i];
                b[ip] -= alpha * b[i];
                for (int j=i;j<N;j++) {
                    A[ip][j]-=alpha*A[i][j];
                }
            }
        }

        // back substitution
        double[] x = new double[N];
        for (int i=N-1;i>=0;i--) {
            double sum=0.0;
            for (int j=i+1;j<N;j++) {
                sum += A[i][j]*x[j];
            }
            x[i]=(b[i]-sum)/A[i][i];
        }
        return x;
    }
    
    public static double[][] InverseMatrix(double[][] A){
        int N=A.length;
        double[][] Ainv=new double[N][N];
        double[] b=new double[N];
        double[] x=new double[N];
        for(int i=0;i<N;i++) b[i]=0.0;
        for(int i=0;i<N;i++){
            b[i]=1.0;
            if(i>=1) b[i-1]=0.0;
            x=GaussianElimination(A,b);
            for(int j=0;j<N;j++) Ainv[j][i]=x[j];
        }
        return Ainv;
    }
    
    public static double[][] multMatrix(double[][] A, double[][] B){
        int N=A.length;
        double[][] C=new double[N][N];
        for(int i=0;i<N;i++){
            for(int j=0;j<N;j++){
                double sum=0.0d;
                for(int k=0;k<N;k++) sum+=A[i][k]*B[k][j];
                C[i][j]=sum;
            }
        }
        return C;
    }
/*
    public static void printMatrix(double[][] A){
        NumberFormat nf = NumberFormat.getNumberInstance();
        nf.setMaximumIntegerDigits(4);
        nf.setMinimumFractionDigits(4);
        int len=A.length;
        for(int i=0;i<len;i++){
            String str="";
            for(int j=0;j<len;j++){
                str+=" Mat["+i+","+j+"]="+nf.format(A[i][j])+",";
            }
            System.out.println(str);
        }
    }
*/
    public static void print(double[][] A, double[] b, double[] x){
        NumberFormat nf = NumberFormat.getNumberInstance();
        nf.setMaximumIntegerDigits(4);
        nf.setMinimumFractionDigits(4);
        System.out.println("** Print result of GaussianElimination");
        int len=A.length;
        for(int i=0;i<len;i++){
            String str="";
            for(int j=0;j<len;j++){
                str+=" A["+i+","+j+"]="+nf.format(A[i][j])+",";
            }
            System.out.println(str);
        }
        
        String str="";
        for(int j=0;j<len;j++){
            str+=" b["+j+"]="+nf.format(b[j])+",";
        }
        System.out.println(str);
        
        str=" Solution;";
        for(int j=0;j<len;j++){
            str+=" x["+j+"]="+nf.format(x[j])+",";
        }
        System.out.println(str);
        
        str=" Accuracy ";
        for(int i=0;i<len;i++){
            double Ax=0;
            for(int j=0;j<len;j++) Ax+=A[i][j]*x[j];
            double err=Math.abs(Ax-b[i]);
            str+="err["+i+"]="+err+", ";
        }
        System.out.println(str);
    }
}
