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
*     A = PLU

   Zludpp represents P as a pivot array (see 
   Pivot.java ), L as a Zltmat, and U as a Zutmat.

   @version Pre-alpha
   @author G. W. Stewart
*/

public class Zludpp{


/** The number of rows in L */
   public int nrl;
/** The number of columns in L */
   public int ncl;
/** The number of rows in U */
   int nru;
/** The number of columns in U */
   int ncu;
/** The pivot array (see  Pivot.java ) */
   public int pvt[];
/** The lower triangular matrix L */
   public Zltmat L;
/** The upper triangular matrix U */
   public Zutmat U;

/**
   Computes the partially pivoted LU decompostion.

   @param     A A Zmat
   @return    The Zludpp of A
   @exception 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 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) // Copy U from T.
         for (i=0; 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 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;
         }
   }
}