/*  timing module for LU factorization  */

import java.util.Random;
import jnt.linear_algebra.LU;
import jnt.linear_algebra.JBLAS;
import jnt.util.Stopwatch;
import jnt.io.NumericIO;


public class tLU
{
    static private final int DEBUG=0;

    public static void main(String args[])
    {
        double A[][];
        int N;
        
        if (args.length < 1)
        {
            System.out.println("Usage -f matrix-filename reps");
            System.out.println("      N reps (use random-matrices)");
            return;
        }
        
    
        if (args[0].compareTo("-f") == 0) // initialize usinge file
        {   
            A = NumericIO.readMatrixDouble(args[1]);
            N = A.length;
        }
        else // initialize using random matrices
        {
            N = Integer.valueOf(args[0]).intValue();
            A = new double[N][N];
        
            Random R = new Random();    
            for (int i=0; i<N; i++)
                for (int j=0; j<N; j++)
                    A[i][j] = 1.0 - 2.0 * R.nextDouble();   
							// range between -1 and 1
        }   
        
        int reps = Integer.valueOf(args[1]).intValue();
        System.out.println("LU benchmark: N= " + N + "   reps= " + reps);

        double b[] = new double[N];
        Stopwatch Q = new Stopwatch();

        double ops = (2.0 * N * N * N)/3.0   + 2.0 * N * N;
        


        // initailize the right hand side to the vector {1, 1, ... 1}
        //  
        for (int j=0; j<N; j++)
            b[j] = 1.0;

        if (DEBUG > 0)
        {
            System.out.println("A: ");
            NumericIO.write(A);
        
            System.out.println("b: ");
            NumericIO.write(b);
        }       
        
        // first, make a test to check for accuracy
        
        double tA[][] = JBLAS.new_copy(A);
        int piv[] = new int[N];
        double x[] = JBLAS.new_copy(b); 
        
        LU.factor(tA, piv);
        LU.solve(tA, piv, x);
    
        if (DEBUG > 0)
        {
        
            System.out.println("LU: ");
            NumericIO.write(tA);                
            
            System.out.println("solution x:"); 
            NumericIO.write(x);
            System.out.println("A*x: ");
            NumericIO.write( JBLAS.mult(A,x));
        
        }                   

        System.out.println("LU error |Ax-b|: "  +
             JBLAS.norm2(  JBLAS.minus(JBLAS.mult(A, x), b )));


        
        Q.reset();
        Q.start();
        for (int i=0; i<reps; i++)
        {
            LU tmp = new LU(A);
        }
        Q.stop();

        System.out.println("Time for " + reps+ " factors: "
					+ Q.read()+" seconds.");
        if (Q.read() > 0.0)
            System.out.println("        Mflops: " + 
                    (reps*1.0e-6)*ops / Q.read());
        else
            System.out.println("Time too small.");
            

        // now time triangular solve part
        //
        LU LUA = new LU(A);
        Q.reset();
        Q.start();
        for (int i=0; i<reps; i++)
        {
            x = LUA.solve(b);
        }
        Q.stop();
        
        System.out.println("Time for "+ reps + " solves: "
				+ Q.read()+" seconds.");
        if (Q.read() > 0.0) 
            System.out.println("       Mflops: " + 
                (reps*1.0e-6)* N*N / Q.read());
        else
            System.out.println("Time too small.");
            

    }
}