package jnt.scimark2;/**    LU matrix factorization. (Based on TNT implementation.)    Decomposes a matrix A  into a triangular lower triangular    factor (L) and an upper triangular factor (U) such that    A = L*U.  By convnetion, the main diagonal of L consists    of 1's so that L and U can be stored compactly in    a NxN matrix.*/public class LU {    /**        Returns a <em>copy</em> of the compact LU factorization.        (useful mainly for debugging.)        @return the compact LU factorization.  The U factor        is stored in the upper triangular portion, and the L        factor is stored in the lower triangular portion.        The main diagonal of L consists (by convention) of        ones, and is not explicitly stored.    */	public static final double num_flops(int N)	{		// rougly 2/3*N^3		double Nd = (double) N;		return (2.0 * Nd *Nd *Nd/ 3.0);	}    protected static double[] new_copy(double x[])    {        int N = x.length;        double T[] = new double[N];        for (int i=0; i<N; i++)            T[i] = x[i];        return T;    }    protected static double[][] new_copy(double A[][])    {        int M = A.length;        int N = A[0].length;        double T[][] = new double[M][N];        for (int i=0; i<M; i++)        {            double Ti[] = T[i];            double Ai[] = A[i];            for (int j=0; j<N; j++)                Ti[j] = Ai[j];        }        return T;    }    public static int[] new_copy(int x[])    {        int N = x.length;        int T[] = new int[N];        for (int i=0; i<N; i++)            T[i] = x[i];        return T;    }    protected static final void insert_copy(double B[][], double A[][])    {        int M = A.length;        int N = A[0].length;		int remainder = N & 3;		 // N mod 4;        for (int i=0; i<M; i++)        {            double Bi[] = B[i];            double Ai[] = A[i];			for (int j=0; j<remainder; j++)                Bi[j] = Ai[j];            for (int j=remainder; j<N; j+=4)			{				Bi[j] = Ai[j];				Bi[j+1] = Ai[j+1];				Bi[j+2] = Ai[j+2];				Bi[j+3] = Ai[j+3];			}		}            }    public double[][] getLU()    {        return new_copy(LU_);    }    /**        Returns a <em>copy</em> of the pivot vector.        @return the pivot vector used in obtaining the        LU factorzation.  Subsequent solutions must        permute the right-hand side by this vector.    */    public int[] getPivot()    {        return new_copy(pivot_);    }        /**        Initalize LU factorization from matrix.        @param A (in) the matrix to associate with this                factorization.    */    public LU( double A[][] )    {        int M = A.length;        int N = A[0].length;        //if ( LU_ == null || LU_.length != M || LU_[0].length != N)            LU_ = new double[M][N];        insert_copy(LU_, A);        //if (pivot_.length != M)            pivot_ = new int[M];        factor(LU_, pivot_);    }    /**        Solve a linear system, with pre-computed factorization.        @param b (in) the right-hand side.        @return solution vector.    */    public double[] solve(double b[])    {        double x[] = new_copy(b);        solve(LU_, pivot_, x);        return x;    }    /**    LU factorization (in place).    @param A (in/out) On input, the matrix to be factored.        On output, the compact LU factorization.    @param pivit (out) The pivot vector records the        reordering of the rows of A during factorization.            @return 0, if OK, nozero value, othewise.*/public static int factor(double A[][],  int pivot[]){     int N = A.length;    int M = A[0].length;    int minMN = Math.min(M,N);    for (int j=0; j<minMN; j++)    {        // find pivot in column j and  test for singularity.        int jp=j;                double t = Math.abs(A[j][j]);        for (int i=j+1; i<M; i++)        {            double ab = Math.abs(A[i][j]);            if ( ab > t)            {                jp = i;                t = ab;            }        }                pivot[j] = jp;        // jp now has the index of maximum element         // of column j, below the diagonal        if ( A[jp][j] == 0 )                             return 1;       // factorization failed because of zero pivot        if (jp != j)        {            // swap rows j and jp            double tA[] = A[j];            A[j] = A[jp];            A[jp] = tA;        }        if (j<M-1)                // compute elements j+1:M of jth column        {            // note A(j,j), was A(jp,p) previously which was            // guarranteed not to be zero (Label #1)            //            double recp =  1.0 / A[j][j];            for (int k=j+1; k<M; k++)                A[k][j] *= recp;        }        if (j < minMN-1)        {            // rank-1 update to trailing submatrix:   E = E - x*y;            //            // E is the region A(j+1:M, j+1:N)            // x is the column vector A(j+1:M,j)            // y is row vector A(j,j+1:N)            for (int ii=j+1; ii<M; ii++)            {                double Aii[] = A[ii];                double Aj[] = A[j];                double AiiJ = Aii[j];                for (int jj=j+1; jj<N; jj++)                  Aii[jj] -= AiiJ * Aj[jj];            }        }    }    return 0;}    /**        Solve a linear system, using a prefactored matrix            in LU form.        @param LU (in) the factored matrix in LU form.         @param pivot (in) the pivot vector which lists            the reordering used during the factorization            stage.        @param b    (in/out) On input, the right-hand side.                    On output, the solution vector.    */    public static void solve(double LU[][], int pvt[], double b[])    {        int M = LU.length;        int N = LU[0].length;        int ii=0;        for (int i=0; i<M; i++)        {            int ip = pvt[i];            double sum = b[ip];            b[ip] = b[i];            if (ii==0)                for (int j=ii; j<i; j++)                    sum -= LU[i][j] * b[j];            else                 if (sum == 0.0)                    ii = i;            b[i] = sum;        }        for (int i=N-1; i>=0; i--)        {            double sum = b[i];            for (int j=i+1; j<N; j++)                sum -= LU[i][j] * b[j];            b[i] = sum / LU[i][i];        }    }                   private double LU_[][];    private int pivot_[];}