// simeq_lup.java   LU with pivot                                     
// solve for X,  A X = Y, using L U = A decomposition 

public class simeq_lup
{
  simeq_lup(int n, double A[], double Y[], double X[])
  {
    double B[] = new double[n*(n+1)];     //  WORKING MATRIX 
    int ROW[] = new int[n];               // ROW INTERCHANGE INDICIES 
    double by[] = new double[n];          // interchanged Y 
    int HOLD , I_PIVOT;  // PIVOT INDICIES 
    double pivot;        // PIVOT ELEMENT VALUE 
    int i, j, k, m, im, it;

    for(m=0; m<n; m++) ROW[m] = m;
    for(m=0; m<n-1; m++)
    {
      pivot = Math.abs(A[ROW[m]*n+m]);
      im = m;
      for(i=m+1; i<n; i++)
      {
        if(Math.abs(A[ROW[i]*n+m])>pivot)
        {
          im = i;
          pivot = Math.abs(A[ROW[i]*n+m]);
        }
      }
      it = ROW[im];
      ROW[im] = ROW[m];
      ROW[m] = it;
    }
  
    // copy A into B with largest diagonal 
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        B[i*n+j] = A[ROW[i]*n+j];
      }
      by[i] = Y[ROW[i]];
    }

    for(m=0; m<n; m++)
    {
      // compute column of L in B 
      for(i=m; i<n; i++)
      {
        B[i*n+m] = B[i*n+m];
        for(k=0; k<m; k++)
        {
          B[i*n+m] = B[i*n+m] - B[i*n+k]*B[k*n+m];
        }
      }
      // compute row of U in B 
      for(j=m+1; j<n; j++)
      {
        B[m*n+j] = B[m*n+j];
        for(k=0; k<m; k++)
        {
          B[m*n+j] = B[m*n+j] - B[m*n+k]*B[k*n+j];
        }
        B[m*n+j] = B[m*n+j]/B[m*n+m];
      }
    }

    // solve for X given by 
    for(i=0; i<n; i++)
    {
      for(k=0; k<i; k++)
      {
        by[i] = by[i] - B[i*n+k]*by[k];
      }
      by[i] = by[i]/B[i*n+i];
    }

    for(j=n-1; j>=0; j--)
    {
      X[j] = by[j];
      for(k=j+1; k<n; k++)
      {
        X[j] = X[j] - B[j*n+k]*X[k];
      }
    }
  } // end simeq_lup
} // end simeq_lup.java