src/jnt/scimark2/LU.java

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_[];
}