protected static double ludcmp(double[][] A,int[] index){
  final double TINY=1.0e-20;
  double[] vv=new double[4];
  double d=1;
  double temp, sum;
  for (int i=0; i < 4; i+=1) {
    double big=0;
    for (int j=0; j < 4; j+=1) {
      if ((temp=Math.abs(A[i][j])) > big) {
        big=temp;
      }
    }
    if (big == 0) {
      return 0;
    }
 else {
      vv[i]=1 / big;
    }
  }
  for (int j=0; j < 4; j+=1) {
    for (int i=0; i < j; i+=1) {
      sum=A[i][j];
      for (int k=0; k < i; k+=1) {
        sum-=A[i][k] * A[k][j];
      }
      A[i][j]=sum;
    }
    double big=0;
    double dum;
    int imax=-1;
    for (int i=j; i < 4; i+=1) {
      sum=A[i][j];
      for (int k=0; k < j; k++) {
        sum-=A[i][k] * A[k][j];
      }
      A[i][j]=sum;
      if ((dum=vv[i] * Math.abs(sum)) >= big) {
        big=dum;
        imax=i;
      }
    }
    if (j != imax) {
      for (int k=0; k < 4; k+=1) {
        dum=A[imax][k];
        A[imax][k]=A[j][k];
        A[j][k]=dum;
      }
      d=-d;
      vv[imax]=vv[j];
    }
    index[j]=imax;
    if (A[j][j] == 0.0)     A[j][j]=TINY;
    if (j != 3) {
      dum=1.0 / A[j][j];
      for (int i=j + 1; i < 4; i+=1) {
        A[i][j]*=dum;
      }
    }
  }
  return d;
}

protected static double ludcmp ( double [ ] [ ] A , int [ ] index ) { final double TINY = 1.0e-20 ; double [ ] vv = new double [ 4 ] ; double d = 1 ; double temp , sum ; for ( int i = 0 ; i < 4 ; i += 1 ) { double big = 0 ; for ( int j = 0 ; j < 4 ; j += 1 ) { if ( ( temp = Math . abs ( A [ i ] [ j ] ) ) > big ) { big = temp ; } } if ( big == 0 ) { return 0 ; } else { vv [ i ] = 1 / big ; } } for ( int j = 0 ; j < 4 ; j += 1 ) { for ( int i = 0 ; i < j ; i += 1 ) { sum = A [ i ] [ j ] ; for ( int k = 0 ; k < i ; k += 1 ) { sum -= A [ i ] [ k ] * A [ k ] [ j ] ; } A [ i ] [ j ] = sum ; } double big = 0 ; double dum ; int imax = - 1 ; for ( int i = j ; i < 4 ; i += 1 ) { sum = A [ i ] [ j ] ; for ( int k = 0 ; k < j ; k ++ ) { sum -= A [ i ] [ k ] * A [ k ] [ j ] ; } A [ i ] [ j ] = sum ; if ( ( dum = vv [ i ] * Math . abs ( sum ) ) >= big ) { big = dum ; imax = i ; } } if ( j != imax ) { for ( int k = 0 ; k < 4 ; k += 1 ) { dum = A [ imax ] [ k ] ; A [ imax ] [ k ] = A [ j ] [ k ] ; A [ j ] [ k ] = dum ; } d = - d ; vv [ imax ] = vv [ j ] ; } index [ j ] = imax ; if ( A [ j ] [ j ] == 0.0 ) A [ j ] [ j ] = TINY ; if ( j != 3 ) { dum = 1.0 / A [ j ] [ j ] ; for ( int i = j + 1 ; i < 4 ; i += 1 ) { A [ i ] [ j ] *= dum ; } } } return d ; }
utility method to perform an lu factorization of a matrix . algorithm derived from " numerical recipes in c " , press et al . , 1988 .
