package com.company;

import Jama.*;
import java.util.Random;

public class Main {

    static double[][] LU;

    public static void main(String[] args) {
        int n = 4;

        // Macierz A
        Matrix A = Matrix.random(n, n);

        Random rand = new Random();
        int isU = rand.nextInt((1 - 0) + 1) + 0;


        // L lub U - uzupelnianie zerami
        if (isU == 0) {
            System.out.println("U");
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < i; j++) {
                    A.set(i, j, 0);
                }
            }
        } else {
            System.out.println("L");
            for (int j = 0; j < n; j++) {
                for (int i = 0; i < j; i++) {
                    A.set(i, j, 0);
                }
            }
        }

        // Macierz B - wektor
        Matrix B = Matrix.random(n, 1);

        System.out.println("Macierz A:");
        A.print(5,3);
        System.out.println("\nMacierz B:");
        System.out.println(" ");
        B.print(5,3);

        // Macierz odwrotna dla testu
        System.out.println("Rozwiazanie z biblioteki:");
        Matrix  G = A.solve(B);
        G.print(5,3);
        System.out.println(" ");

        // Solution
        System.out.println("Rozwiazanie:");
        Matrix S = FindInverse(A, B, n);

        S.print(5, 3);


        double Af = 0;
        for(int i = 0; i < S.getRowDimension(); i++) {
            for(int j = 0; j < S.getColumnDimension(); j++) {
                Af += S.get(i, j) * S.get(i, j);
            }
        }


        Af = Math.sqrt(Af);
        //System.out.println("\nWartość Af: " + Af);

        // --------------------------------

        Matrix identityMatrix = Matrix.identity(n, n);

        Matrix[] X = new Matrix[n];

        for(int i = 0; i < n; i++) {
            X[i] = FindInverse(A, identityMatrix.getMatrix(0, n - 1, i, i), n);
        }

        Matrix X1 = FindInverse(A, identityMatrix.getMatrix(0, n - 1, 0, 0), n);
        Matrix X2 = FindInverse(A, identityMatrix.getMatrix(0, n - 1, 1, 1), n);
        Matrix X3 = FindInverse(A, identityMatrix.getMatrix(0, n - 1, 2, 2), n);
        Matrix X4 = FindInverse(A, identityMatrix.getMatrix(0, n - 1, 3, 3), n);

        Matrix IKS = new Matrix(n,n);

        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                IKS.set(i, j, X[j].get(i, 0));
            }
        }

        IKS.print(5, 3);

        // po pomnozeniu otrzymanej odwrotnej X przez A
        Matrix newIdentity = A.times(IKS);

        newIdentity.print(10, 10);


        double Af2 = 0;
        for(int i = 0; i < newIdentity.getRowDimension(); i++) {
            for(int j = 0; j < newIdentity.getColumnDimension(); j++) {
                Af2 += newIdentity.get(i, j) * newIdentity.get(i, j);
            }
        }

        Af2 = Math.sqrt(Af2);
        System.out.println("\nWartość Af: " + Af2);

    }

    public static Matrix FindInverse(Matrix A, Matrix B, int n) {

        // pivot generowanie
        int piv[] = getPivot(A, n);

        int columnNumber = 1;
        Matrix matrixX = B.getMatrix(piv,0,columnNumber-1);
        double[][] X = matrixX.getArray();


        // Solve L*Y = B(piv,:)
        for (int k = 0; k < n; k++) {
            for (int i = k+1; i < n; i++) {
                for (int j = 0; j < columnNumber; j++) {
                    X[i][j] -= X[k][j]*LU[i][k];
                }
            }
        }
        // Solve U*X = Y;
        for (int k = n-1; k >= 0; k--) {
            for (int j = 0; j < columnNumber; j++) {
                X[k][j] /= LU[k][k];
            }
            for (int i = 0; i < k; i++) {
                for (int j = 0; j < columnNumber; j++) {
                    X[i][j] -= X[k][j]*LU[i][k];
                }
            }
        }
        return matrixX;
    }


    public static int[] getPivot(Matrix A, int n) {
        LU = A.getArrayCopy();
        int piv[] = new int[n];

        for (int i = 0; i < n; i++) {
            piv[i] = i;
        }
        double[] LUrowi;
        double[] LUcolj = new double[n];

        // Outer loop.
        for (int j = 0; j < n; j++) {

            // Make a copy of the j-th column to localize references.
            for (int i = 0; i < n; i++) {
                LUcolj[i] = LU[i][j];
            }

            // Apply previous transformations.
            for (int i = 0; i < n; i++) {
                LUrowi = LU[i];

                // Most of the time is spent in the following dot product.
                int kmax = Math.min(i,j);
                double s = 0.0;
                for (int k = 0; k < kmax; k++) {
                    s += LUrowi[k]*LUcolj[k];
                }

                LUrowi[j] = LUcolj[i] -= s;
            }

            // Find pivot and exchange if necessary.
            int p = j;
            for (int i = j+1; i < n; i++) {
                if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
                    p = i;
                }
            }
            if (p != j) {
                for (int k = 0; k < n; k++) {
                    double t = LU[p][k]; LU[p][k] = LU[j][k]; LU[j][k] = t;
                }
                int k = piv[p];
                piv[p] = piv[j];
                piv[j] = k;
            }

            if (j < n & LU[j][j] != 0.0) {
                for (int i = j+1; i < n; i++) {
                    LU[i][j] /= LU[j][j];
                }
            }
        }

        return  piv;
    }
}