content: package Jampack;
  
  /**
   * Zludpp implements the LU decomposition with partial
   * pivoting.  Specifically, given a matrix A, there
   * is a permunation matrix P, a unit lower triangular matrix
   * L whose subdiagonal elements are less than one in magnitude
   * and a upper triangular matrix U such that
   * <pre>
   *     A = PLU
   * </pre>
   * Zludpp represents P as a pivot array (see <a href="Jampack.Pivot.html">
   * Pivot.java </a>), L as a Zltmat, and U as a Zutmat.
   *
   * @author G. W. Stewart
   * @version Pre-alpha
   */
  
  public class Zludpp {
  
  
      /**
       * The number of rows in L
       */
      public int nrl;
      /**
       * The number of columns in L
       */
      public int ncl;
      /**
       * The pivot array (see <a href="Jampack.Pivot.html"> Pivot.java </a>)
       */
      public int pvt[];
      /**
       * The lower triangular matrix L
       */
      public Zltmat L;
      /**
       * The upper triangular matrix U
       */
      public Zutmat U;
      /**
       * The number of rows in U
       */
      int nru;
      /**
       * The number of columns in U
       */
      int ncu;
  
      /**
       * Computes the partially pivoted LU decompostion.
       *
       * @param A A Zmat
       * @return The Zludpp of A
       * @throws JampackException Passed from below.
       */
  
      public Zludpp(Zmat A)
              throws JampackException {
          int i, j, k, nr, nc;
          double absi, mx, t;
          Zmat T;
          Z Tk[];
  
          A.getProperties();
  
        /* Set up L and U */
  
          nr = A.nr;
          nrl = nr;
          nc = A.nc;
          ncl = Math.min(A.nr, A.nc);
          nru = ncl;
          ncu = nc;
  
  
          L = new Zltmat(nrl, ncl);
          U = new Zutmat(nru, ncu);
          pvt = new int[ncl];
  
  
        /* Set up the matrix T in which the elimination will be
           performed and copy A to it.*/
  
          if (nrl >= ncu)
              T = L;
          else
              T = U;
  
          for (i = 0; i < nr; i++) {
              for (j = 0; j < nc; j++) {
                  T.re[i][j] = A.re[i][j];
                  T.im[i][j] = A.im[i][j];
              }
          }
  
        /* Outer elimination loop. */
  
          Tk = new Z[nrl];   // This should be repaced by a Z1.
  
          for (k = 0; k < Math.min(nr, nc); k++) {
  
           /* Find the pivot row. */
  
              mx = 0.;
              pvt[k] = k;
              for (i = k; i < nr; i++) {
                  Tk[i] = T.get0(i, k);
                  if ((absi = Z.abs(Tk[i])) > mx) {
                      pvt[k] = i;
                      mx = absi;
                  }
              }
              if (mx == 0.0) continue;
  
           /* Perform the exchange. */
  
              Tk[k].Exch(Tk[pvt[k]]);
              for (j = 0; j < nc; j++) {
                  t = T.re[k][j];
                  T.re[k][j] = T.re[pvt[k]][j];
                  T.re[pvt[k]][j] = t;
                  t = T.im[k][j];
                  T.im[k][j] = T.im[pvt[k]][j];
                  T.im[pvt[k]][j] = t;
              }
  
           /* Compute multipliers and eliminate. */
  
              for (i = k + 1; i < nr; i++) {
                  T.put0(i, k, Tk[i].Div(Tk[i], Tk[k]));
                  for (j = k + 1; j < nc; j++) {
                      T.re[i][j] = T.re[i][j]
                              - T.re[i][k] * T.re[k][j] + T.im[i][k] * T.im[k][j];
                      T.im[i][j] = T.im[i][j]
                              - T.im[i][k] * T.re[k][j] - T.re[i][k] * T.im[k][j];
                  }
              }
          }
  
        /* Finalize L and U */
  
          if (nr >= nc) // Copy U from T.
              for (i = 0; i < nc; i++) {
                  for (j = 0; j < nc; j++)
                      if (i > j) {
                          U.re[i][j] = 0.0;
                          U.im[i][j] = 0.0;
                      } else {
                          U.re[i][j] = T.re[i][j];
                          U.im[i][j] = T.im[i][j];
                          L.re[i][j] = 0.0;
                          L.im[i][j] = 0.0;
                      }
                  L.re[i][i] = 1.0;
                  L.im[i][i] = 0.0;
              }
          else // Copy L from T.
              for (i = 0; i < nr; i++) {
                  for (j = 0; j < nr; j++)
                      if (i > j) {
                          L.re[i][j] = T.re[i][j];
                          L.im[i][j] = T.im[i][j];
                          U.re[i][j] = 0.0;
                          U.im[i][j] = 0.0;
                      } else {
                          L.re[i][j] = 0.0;
                          L.im[i][j] = 0.0;
                      }
                  L.re[i][i] = 1.0;
                  L.im[i][i] = 0.0;
              }
      }
  }
  