public class LUDecomp {

	double a[][];
	double b[];
	double x[];
	double ops,norma,resid,normx;
	int ipvt[];
	int n, lda,ldaa,info;
	
	public void initialize(int size){
		n= size;
	    ldaa = n; 
	    lda = ldaa + 1;
	    a = new double[ldaa][lda];
	    b = new double [ldaa];
	    x = new double [ldaa];
	    ipvt = new int [ldaa];
	    long nl = (long) size;   //avoid integer overflow
	    ops = (2.0*(nl*nl*nl))/3.0 + 2.0*(nl*nl);
	    norma = matgen(a,lda,size,b);   
	}
	
	public void run(){
		 info = dgefa(a,lda,n,ipvt);
	     dgesl(a,lda,n,ipvt,b,0);
	}
	
	public void validate(){
	   int i;
	   double eps,residn;
  	   final double ref = 6.0; 

	   for (i = 0; i < n; i++) {
	      x[i] = b[i];
	   }
	   norma = matgen(a,lda,n,b);
	   for (i = 0; i < n; i++) {
	      b[i] = -b[i];
       }

	   dmxpy(n,b,n,lda,x,a);
	   resid = 0.0;
	   normx = 0.0;
       for (i = 0; i < n; i++) {
	      resid = (resid > abs(b[i])) ? resid : abs(b[i]);
	      normx = (normx > abs(x[i])) ? normx : abs(x[i]);
	   }

	   eps =  epslon((double)1.0);
	   residn = resid/( n*norma*normx*eps );  
	   if (residn > ref) {
		  System.out.println("Validation failed");
          System.out.println("Computed Norm Res = " + residn);
	      System.out.println("Reference Norm Res = " + ref); 
	   }
	}
	
	// estimate unit roundoff in quantities of size x.
	private double epslon (double x) {
	    double a,b,c,eps;
	    a = 4.0e0/3.0e0;
	    eps = 0;
	    while (eps == 0) {
	       b = a - 1.0;
	       c = b + b + b;
	       eps = abs(c-1.0);
	    }
	    return(eps*abs(x));
	}
	
	// multiply matrix m times vector x and add the result to vector y
    private void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][]) {
	    int j,i;

	    // cleanup odd vector
	    for (j = 0; j < n2; j++) {
	       for (i = 0; i < n1; i++) {
		      y[i] += x[j]*m[j][i];
	       }
	    }
	}
	
    //dgefa factors a double precision matrix by gaussian elimination.
	private int dgefa( double a[][], int lda, int n, int ipvt[]) {
	    double[] col_k, col_j;
	    double t;
	    int j,k,kp1,l,nm1;
	    int info;
	    
	    // gaussian elimination with partial pivoting
	    info = 0;
	    nm1 = n - 1;
	    if (nm1 >=  0) {
	       for (k = 0; k < nm1; k++) {
		      col_k = a[k];
		      kp1 = k + 1;
		
		      // find l = pivot index
		      l = idamax(n-k,col_k,k,1) + k;
		      ipvt[k] = l;
		
		      // zero pivot implies this column already triangularized
              if (col_k[l] != 0) {
		  
		         // interchange if necessary
		         if (l != k) {
		            t = col_k[l];
		            col_k[l] = col_k[k];
		            col_k[k] = t;
		         }
		  
		      // compute multipliers
		      t = -1.0/col_k[k];
		      dscal(n-(kp1),t,col_k,kp1,1);
		  
		      // row elimination with column indexing
		      for (j = kp1; j < n; j++) {
		         col_j = a[j];
		         t = col_j[l];
		         if (l != k) {
		            col_j[l] = col_j[k];
		            col_j[k] = t;
		         }
		         daxpy(n-(kp1),t,col_k,kp1,1,
			     col_j,kp1,1);
		      }
		   } else {
		      info = k;
		   }
	    }
	  }
	  ipvt[n-1] = n-1;
	  if (a[(n-1)][(n-1)] == 0) info = n-1;
	    
	  return info;
	}
	
	/* dgesl solves the double precision system
       a * x = b  or  trans(a) * x = b
       using the factors computed by dgeco or dgefa. */
	private void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job) {
	    double t;
	    int k,kb,l,nm1,kp1;

	    nm1 = n - 1;
	    if (job == 0) {

	       // job = 0 , solve  a * x = b.  first solve  l*y = b
	       if (nm1 >= 1) {
		      for (k = 0; k < nm1; k++) {
		         l = ipvt[k];
		         t = b[l];
		         if (l != k){
		            b[l] = b[k];
		            b[k] = t;
		         }
		         kp1 = k + 1;
		         daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
		      }
	       }

	       // now solve  u*x = y
	       for (kb = 0; kb < n; kb++) {
		      k = n - (kb + 1);
		      b[k] /= a[k][k];
		      t = -b[k];
		      daxpy(k,t,a[k],0,1,b,0,1);
	       }
	    } else {

	       // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
	       for (k = 0; k < n; k++) {
		      t = ddot(k,a[k],0,1,b,0,1);
		      b[k] = (b[k] - t)/a[k][k];
	       }

   	       // now solve trans(l)*x = y 
           if (nm1 >= 1) {
		      for (kb = 1; kb < nm1; kb++) {
		         k = n - (kb+1);
		         kp1 = k + 1;
		         b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
		         l = ipvt[k];
		         if (l != k) {
		            t = b[l];
		            b[l] = b[k];
		            b[k] = t;
		         }
		      }
	       }
	    }
	}
	  
    //forms the dot product of two vectors.
	private double ddot( int n, double dx[], int dx_off, int incx, double dy[],
	       int dy_off, int incy) {
    
		double dtemp;
		int i,ix,iy;
        dtemp = 0;

        if (n > 0) {
           if (incx != 1 || incy != 1) {

      	      // code for unequal increments or equal increments not equal to 1
              ix = 0;
        	  iy = 0;
	          if (incx < 0) ix = (-n+1)*incx;
	          if (incy < 0) iy = (-n+1)*incy;
	          for (i = 0;i < n; i++) {
	             dtemp += dx[ix +dx_off]*dy[iy +dy_off];
	             ix += incx;
	             iy += incy;
	          }
           } else {

	          // code for both increments equal to 1
	          for (i=0;i < n; i++)
	             dtemp += dx[i +dx_off]*dy[i +dy_off];
           }
        }
        return(dtemp);
    }

	//finds the index of element having max absolute value.
	private  int idamax( int n, double dx[], int dx_off, int incx){
    
		double dmax, dtemp;
		int i, ix, itemp=0;

		if (n < 1) {
           itemp = -1;
		} else if (n ==1) {
           itemp = 0;
		} else if (incx != 1) {
           // code for increment not equal to 1
           dmax = abs(dx[0 +dx_off]);
           ix = 1 + incx;
           for (i = 1; i < n; i++) {
	          dtemp = abs(dx[ix + dx_off]);
	          if (dtemp > dmax)  {
	             itemp = i;
	             dmax = dtemp;
	          }
	          ix += incx;
           }
        } else {
           // code for increment equal to 1
           itemp = 0;
           dmax = abs(dx[0 +dx_off]);
           for (i = 1; i < n; i++) {
	          dtemp = abs(dx[i + dx_off]);
	          if (dtemp > dmax) {
	             itemp = i;
	             dmax = dtemp;
	          }
           }
       }
       return (itemp);
	}

	private double abs (double d) {
	    return (d >= 0) ? d : -d;
	}
	
    //scales a vector by a constant
    private void dscal( int n, double da, double dx[], int dx_off, int incx){
    
    	int i,nincx;

        if (n > 0) {
           if (incx != 1) {

	          // code for increment not equal to 1
              nincx = n*incx;
	          for (i = 0; i < nincx; i += incx)
	             dx[i +dx_off] *= da;
           } else {

	          // code for increment equal to 1
	          for (i = 0; i < n; i++)
	             dx[i +dx_off] *= da;
           }
        }
    } 
    
    //constant times a vector plus a vector.
    private void daxpy( int n, double da, double dx[], int dx_off, int incx,
	      double dy[], int dy_off, int incy) {
    	
    	int i,ix,iy;

        if ((n > 0) && (da != 0)) {
           if (incx != 1 || incy != 1) {

	          // code for unequal increments or equal increments not equal to 1
              ix = 0;
	          iy = 0;
	          if (incx < 0) ix = (-n+1)*incx;
	          if (incy < 0) iy = (-n+1)*incy;
	          for (i = 0;i < n; i++) {
	             dy[iy +dy_off] += da*dx[ix +dx_off];
	             ix += incx;
	             iy += incy;
	          }
	          return;
           } else {

	          // code for both increments equal to 1
              for (i=0; i < n; i++)
	             dy[i +dy_off] += da*dx[i +dx_off];
           }
        }
    }
	
	private double matgen (double a[][], int lda, int n, double b[]){
	    double norma;
	    int init, i, j;
	    init = 1325;
	    norma = 0.0;
	    for (i = 0; i < n; i++) {
	       for (j = 0; j < n; j++) {
		      init = 3125*init % 65536;
		      a[j][i] = (init - 32768.0)/16384.0;
		      norma = (a[j][i] > norma) ? a[j][i] : norma;
	       }
	    }
	    for (i = 0; i < n; i++) {
	       b[i] = 0.0;
	    }
	    for (j = 0; j < n; j++) {
	       for (i = 0; i < n; i++) {
		      b[i] += a[j][i];
	       }
	    }
	    return norma;
	}
	
	public static void main(String[] args) {	
		LUDecomp lu = new LUDecomp();
		
		lu.initialize(500);
		lu.run();
		lu.validate();
	}
}