content: package Jampack;
  
  /**
   * Zsvd implements the singular value decomposion of a Zmat.
   * Specifically if X is an mxn matrix with m&gt;=n there are unitary
   * matrices U and V such that
   * <pre>
   *     U^H*X*V = | S |
   *               | 0 |
   * </pre>
   * where S = diag(s1,...,sm) with
   * <pre>
   *     s1 >= s2 >= ... >= sn >=0.
   * </pre>
   * If m&lt;n the decomposition has the form
   * <pre>
   *     U^H*X*V = | S  0 |,
   * </pre>
   * <p/>
   * where S is diagonal of order m.  The diagonals of S are the
   * singular values of A.  The columns of U are the left singular
   * vectors of A and the columns of V are the right singular vectors.
   *
   * @author G. W. Stewart
   * @version Pre-alpha
   */
  
  public class Zsvd {
  
      /**
       * Limits the number of iterations in the SVD algorithm
       */
      public static int MAXITER = 30;
  
      /**
       * The matrix of left singular vectors
       */
      public Zmat U;
  
      /**
       * The matrix of right singular vectore
       */
      public Zmat V;
  
      /**
       * The diagonal matrix of singular values
       */
  
      public Zdiagmat S;
  
      /**
       * Computes the SVD of a Zmat XX.  Throws a JampackException
       * if the maximum number of iterations is exceeded.
       *
       * @param XX A Zmat
       * @return The Zsvd of XX
       * @throws JampackException Thrown if maximimum number of iterations is
       *                          exceeded.<br>
       *                          Passed from below.
       */
  
      public Zsvd(Zmat XX)
              throws JampackException {
  
          int i, il, iu, iter, j, k, kk, m, mc;
  
          double as, at, au, axkk, axkk1, dmax, dmin, ds, ea,
                  es, shift, ss, t, tre;
  
          Z xkk, xkk1, xk1k1, ukj, vik1;
  
          Rot P = new Rot();
  
  /*    Initialization */
  
          Z scale = new Z();
          Z zr = new Z();
  
          Zmat X = new Zmat(XX);
  
          Z1 h;
          Z1 temp = new Z1(Math.max(X.nr, X.nc));
  
          mc = Math.min(X.nr, X.nc);
          double d[] = new double[mc];
          double e[] = new double[mc];
  
          S = new Zdiagmat(mc);
          U = Eye.o(X.nr);
          V = Eye.o(X.nc);
  
          m = Math.min(X.rx, X.cx);
        
  
  /*
        Reduction to Bidiagonal form.
  */
  
          for (k = X.bx; k <= m; k++) {
  
              h = House.genc(X, k, X.rx, k);
              House.ua(h, X, k, X.rx, k + 1, X.cx, temp);
              House.au(U, h, U.bx, U.rx, k, U.cx, temp);
  
              if (k != X.cx) {
                  h = House.genr(X, k, k + 1, X.cx);
                  House.au(X, h, k + 1, X.rx, k + 1, X.cx, temp);
                  House.au(V, h, V.bx, V.rx, k + 1, V.cx, temp);
              }
          }
  
  /*
        Scale the bidiagonal matrix so that its elements are
        real.
  */
  
          for (k = X.bx; k <= m; k++) {
              kk = k - X.bx;
              xkk = X.get(k, k);
              axkk = Z.abs(xkk);
              X.put(k, k, new Z(axkk));
              d[kk] = axkk;
              scale.Div(scale.Conj(xkk), axkk);
              if (k < X.cx) {
                  xkk1 = X.get(k, k + 1);
                  X.put(k, k + 1, xkk1.Times(scale, xkk1));
              }
              scale.Conj(scale);
              for (i = U.bx; i <= U.rx; i++) {
                  U.put(i, k, zr.Times(U.get(i, k), scale));
              }
  
              if (k < X.cx) {
  
                  xkk1 = X.get(k, k + 1);
                  axkk1 = Z.abs(xkk1);
                  X.put(k, k + 1, new Z(axkk1));
                  e[kk] = axkk1;
                  scale.Div(scale.Conj(xkk1), axkk1);
                  if (k < X.rx) {
                      xk1k1 = X.get(k + 1, k + 1);
                      X.put(k + 1, k + 1, xk1k1.Times(scale, xk1k1));
                  }
                  for (i = V.bx; i <= V.rx; i++) {
                      V.put(i, k + 1, zr.Times(V.get(i, k + 1), scale));
                  }
              }
          }
  
          m = m - X.bx;  // Zero based loops from here on.
  /*
        If X has more columns than rows, rotate out the extra
        superdiagonal element.
  */
          if (X.nr < X.nc) {
              t = e[m];
              for (k = m; k >= 0; k--) {
                  Rot.genr(d[k], t, P);
                  d[k] = P.zr;
                  if (k != 0) {
                      t = P.sr * e[k - 1];
                      e[k - 1] = P.c * e[k - 1];
                  }
                  Rot.ap(V, P, V.bx, V.rx, k + V.bx, X.rx + 1);
                  Rot.ap(X, P, X.bx, X.rx, k + X.bx, X.rx + 1);
              }
          }
  /*
        Caculate the singular values of the bidiagonal matrix.
  */
          iu = m;
          iter = 0;
          while (true) {
  /*
           These two loops determine the rows (il to iu) to
           iterate on.
  */
              while (iu > 0) {
                  if (Math.abs(e[iu - 1]) >
                          1.0e-16 * (Math.abs(d[iu]) + Math.abs(d[iu - 1])))
                      break;
                  e[iu - 1] = 0.;
                  iter = 0;
                  iu = iu - 1;
              }
              iter = iter + 1;
              if (iter > MAXITER) {
                  throw new JampackException
                          ("Maximum number of iterations exceeded.");
              }
              if (iu == 0) break;
  
              il = iu - 1;
              while (il > 0) {
                  if (Math.abs(e[il - 1]) <=
                          1.0e-16 * (Math.abs(d[il]) + Math.abs(d[il - 1])))
                      break;
                  il = il - 1;
              }
              if (il != 0) {
                  e[il - 1] = 0.;
              }
  /*
           Compute the shift (formulas adapted from LAPACK).
  */
              dmax = Math.max(Math.abs(d[iu]), Math.abs(d[iu - 1]));
              dmin = Math.min(Math.abs(d[iu]), Math.abs(d[iu - 1]));
              ea = Math.abs(e[iu - 1]);
              if (dmin == 0.) {
                  shift = 0.;
              } else if (ea < dmax) {
                  as = 1. + dmin / dmax;
                  at = (dmax - dmin) / dmax;
                  au = ea / dmax;
                  au = au * au;
                  shift = dmin * (2. / (Math.sqrt(as * as + au) + Math.sqrt(at * at + au)));
              } else {
                  au = dmax / ea;
                  if (au == 0.) {
                      shift = (dmin * dmax) / ea;
                  } else {
                      as = 1. + dmin / dmax;
                      at = (dmax - dmin) / dmax;
                      t = 1. / (Math.sqrt(1. + (as * au) * (as * au)) +
                              Math.sqrt(1. + (at * au) * (at * au)));
                      shift = (t * dmin) * au;
                  }
              }
  /*
          Perform the implicitly shifted QR step.
  */
              t = Math.max(Math.max(Math.abs(d[il]), Math.abs(e[il])), shift);
              ds = d[il] / t;
              es = e[il] / t;
              ss = shift / t;
              Rot.genr((ds - ss) * (ds + ss), ds * es, P);
              for (i = il; i < iu; i++) {
                  t = P.c * d[i] - P.sr * e[i];
                  e[i] = P.sr * d[i] + P.c * e[i];
                  d[i] = t;
                  t = -P.sr * d[i + 1];
                  d[i + 1] = P.c * d[i + 1];
                  Rot.ap(V, P, V.bx, V.rx, V.bx + i, V.bx + i + 1);
                  Rot.genc(d[i], t, P);
                  d[i] = P.zr;
                  t = P.c * e[i] + P.sr * d[i + 1];
                  d[i + 1] = P.c * d[i + 1] - P.sr * e[i];
                  e[i] = t;
                  Rot.aph(U, P, U.bx, U.rx, U.bx + i, U.bx + i + 1);
                  if (i != iu - 1) {
                      t = P.sr * e[i + 1];
                      e[i + 1] = P.c * e[i + 1];
                      Rot.genr(e[i], t, P);
                      e[i] = P.zr;
                  }
              }
          }
  
  /*
        Sort the singular values, setting negative values of d
        to positive.
  */
          for (k = m; k >= 0; k--) {
              if (d[k] < 0) {
                  d[k] = -d[k];
                  for (i = 0; i < X.nc; i++) {
                      V.re[i][k] = -V.re[i][k];
                      V.im[i][k] = -V.im[i][k];
                  }
              }
              for (j = k; j < m; j++) {
                  if (d[j] < d[j + 1]) {
                      t = d[j];
                      d[j] = d[j + 1];
                      d[j + 1] = t;
                      for (i = 0; i < X.nr; i++) {
                          t = U.re[i][j];
                          U.re[i][j] = U.re[i][j + 1];
                          U.re[i][j + 1] = t;
                          t = U.im[i][j];
                          U.im[i][j] = U.im[i][j + 1];
                          U.im[i][j + 1] = t;
                      }
                      for (i = 0; i < X.nc; i++) {
                          t = V.re[i][j];
                          V.re[i][j] = V.re[i][j + 1];
                          V.re[i][j + 1] = t;
                          t = V.im[i][j];
                          V.im[i][j] = V.im[i][j + 1];
                          V.im[i][j + 1] = t;
                      }
                  }
              }
          }
  /*
        Return the decompostion;
  */
          S.re = d;
          return;
      }
  }
  