package Jampack; /** Zsvd implements the singular value decomposion of a Zmat. Specifically if X is an mxn matrix with m>=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<n the decomposition has the form <pre> * U^H*X*V = | S 0 |, </pre> 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. @version Pre-alpha @author G. W. Stewart */ 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 @exception 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; } }