void stpone(SMat A, double[][] wrkptr, double[] rnmp, double[] tolp, int n) { double t, rnm, anorm; double[] alf = wrkptr[6]; /* get initial vector; default is random */ rnm = startv(A, wrkptr, 0, n); if (rnm == 0.0 || ierr != 0) return; /* normalize starting vector */ t = 1.0 / rnm; svd_datx(n, t, wrkptr[0], 1, wrkptr[1], 1); svd_dscal(n, t, wrkptr[3], 1); /* take the first step */ svd_opb(A, wrkptr[3], wrkptr[0], OPBTemp); alf[0] = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1); svd_daxpy(n, -alf[0], wrkptr[1], 1, wrkptr[0], 1); t = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1); svd_daxpy(n, -t, wrkptr[1], 1, wrkptr[0], 1); alf[0] += t; svd_dcopy(n, wrkptr[0], 1, wrkptr[4], 1); rnm = Math.sqrt(svd_ddot(n, wrkptr[0], 1, wrkptr[4], 1)); anorm = rnm + fabs(alf[0]); rnmp[0] = rnm; tolp[0] = reps * anorm; return; }
enum storeVals {STORQ, RETRQ, STORP, RETRP};
mid = svd_idamax(step + 1, bnd, 0, 1); if (fabs(ritz[i-1] - ritz[i]) < eps34 * fabs(ritz[i])) if (bnd[i] > tol && bnd[i-1] > tol) { bnd[i-1] = Math.sqrt(bnd[i] * bnd[i] + bnd[i-1] * bnd[i-1]); if (fabs(ritz[i+1] - ritz[i]) < eps34 * fabs(ritz[i])) if (bnd[i] > tol && bnd[i+1] > tol) { bnd[i+1] = Math.sqrt(bnd[i] * bnd[i] + bnd[i+1] * bnd[i+1]); gap = gapl; if (i < step) gapl = ritz[i+1] - ritz[i]; gap = svd_dmin(gap, gapl); if (gap > bnd[i]) bnd[i] = bnd[i] * (bnd[i] / gap); if (bnd[i] <= 16.0 * eps * fabs(ritz[i])) { neig++; if (!enough[0]) enough[0] = endl < ritz[i] && ritz[i] < endr;
m = svd_imin(A.rows, A.cols); if (dimensions <= 0 || dimensions > m) dimensions = m; write_header(iterations, dimensions, end[0], end[1], true, kappa, A.rows, A.cols, A.vals); if (0 != check_parameters(A, dimensions, iterations, end[0], end[1], true)) { if (A.rows == 0 || A.cols == 0) { R = new SVDRec(); if (SVDVerbosity > 0) printf("TRANSPOSING THE MATRIX FOR SPEED\n"); transpose = true; A = svdTransposeS(A); machar(); // XXX has side effect, computes machine precision fake_memset_127(bnd); OPBTemp = svd_doubleArray(A.rows, false, "las2: OPBTemp"); steps = lanso(A, iterations, dimensions, end[0], end[1], ritz, bnd, wptr, ref_neig, n); int neig = ref_neig[0]; // XXX unwrap neig printf("NUMBER OF LANCZOS STEPS = %6d\n" + "RITZ VALUES STABILIZED = %6d\n", steps + 1, neig); printf("\nCOMPUTED RITZ VALUES (ERROR BNDS)\n");
s = svd_doubleArray(jsq, true, "ritvec: s"); xv2 = svd_doubleArray(n, false, "ritvec: xv2"); svd_dcopy(js, alf, 1, w1, -1); svd_dcopy(steps, bet, 1, 1, w2, 1, -1); // WAS svd_dcopy(steps, &bet[1], 1, &w2[1], -1); imtql2(js, js, w1, w2, s); if (0 != ierr) return 0; // TODO use exception here? for (i = 0; i < n; i++) w1[i] = 0.0; for (i = 0; i < js; i++) { store(n, storeVals.RETRQ, i, w2); svd_daxpy(n, s[tmp], w2, 1, w1, 1); tmp -= js; rotateArray(R.Vt.value, R.Vt.rows * R.Vt.cols, x * R.Vt.cols); R.d = svd_imin(R.d, nsig); for (x = 0; x < R.d; x++) { svd_opb(A, R.Vt.value[x], xv2, OPBTemp); tmp0 = svd_ddot(n, R.Vt.value[x], 1, xv2, 1); svd_daxpy(n, -tmp0, R.Vt.value[x], 1, xv2, 1); tmp0 = Math.sqrt(tmp0); xnorm = Math.sqrt(svd_ddot(n, xv2, 1, xv2, 1)); svd_opa(A, R.Vt.value[x], R.Ut.value[x]); tmp1 = 1.0 / tmp0; svd_dscal(A.rows, tmp1, R.Ut.value[x], 1);
stpone(A, wptr, ref_rnm, ref_tol, n); double tol = ref_tol[0]; // XXX unwrap double rnm = ref_rnm[0]; // XXX unwrap ll = 0; first = 1; last = svd_imin(dimensions + svd_imax(8, dimensions), iterations); ENOUGH = false; double[] ref2_rnm = new double[] { rnm }; // XXX wrap double[] ref2_tol = new double[] { tol }; // XXX wrap j = lanczos_step(A, first, last, wptr, alf, eta, oldeta, bet, ref_ll, ref_ENOUGH, ref2_rnm, ref2_tol, n); ll = ref_ll[0]; // XXX unwrap svd_dcopy(i-l+1, alf, l, 1, ritz, l, -1); // WAS svd_dcopy(i-l+1, &alf[l], 1, &ritz[l], -1); svd_dcopy(i-l, bet, l+1, 1, wrk, l+1, -1); // WAS svd_dcopy(i-l, &bet[l+1], 1, &wrk[l+1], -1); imtqlb(i-l+1, ritz, wrk, bnd, l); // TODO start at l svd_error("svdLAS2: imtqlb failed to converge (ierr = %ld)\n", ierr); svd_error(" l = %ld i = %ld\n", l, i); for (id3 = l; id3 <= i; id3++) svd_error(" %ld %lg %lg %lg\n", id3, ritz[id3], wrk[id3], bnd[id3]); bnd[id3] = rnm * fabs(bnd[id3]); l = i + 1;
public PrincipalComponents (ObjectVector[] vectorInput) { this.vectorInput = vectorInput; this.dimension = vectorInput[0].getVector().getDimension(); double[][] vectorArray = new double[vectorInput.length][dimension]; for (int i = 0; i < vectorInput.length; ++i) { if (vectorInput[i].getVector().getClass() != RealVector.class) { throw new IncompatibleVectorsException( "Principal components class only works with Real Vectors so far!"); } if (vectorInput[i].getVector().getDimension() != dimension) { throw new IncompatibleVectorsException("Dimensions must all be equal!"); } RealVector realVector = (RealVector) vectorInput[i].getVector(); float[] tempVec = realVector.getCoordinates().clone(); for (int j = 0; j < dimension; ++j) { vectorArray[i][j] = (double) tempVec[j]; } } this.matrix = new DMat(vectorArray.length, vectorArray[0].length); matrix.value = vectorArray; System.err.println("Created matrix ... performing svd ..."); Svdlib svd = new Svdlib(); System.err.println("Starting SVD using algorithm LAS2"); svdR = svd.svdLAS2A(Svdlib.svdConvertDtoS(matrix), matrix.cols); }
/***************************************************************** * Function finds the index of element having max. absolute value* based on * FORTRAN 77 routine from Linpack by J. Dongarra * *****************************************************************/ static int svd_idamax(int n, double[] dx, int ix0, int incx) { int ix,imax; double dmax; if (n < 1) return -1; if (n == 1) return 0; if (incx == 0) return -1; ix = (incx < 0) ? ix0 + ((-n+1) * incx) : ix0; imax = ix; dmax = fabs(dx[ix]); for (int i=1; i < n; i++) { ix += incx; double dtemp = fabs(dx[ix]); if (dtemp > dmax) { dmax = dtemp; imax = ix; } } return imax; }
double[] bndn = new double[n]; System.arraycopy(bnd, offset, bndn, 0, n); imtqlb(n, dn, en, bndn); System.arraycopy(dn, 0, d, offset, n); System.arraycopy(en, 0, e, offset, n);
m = svd_imin(A.rows, A.cols); if (dimensions <= 0 || dimensions > m) dimensions = m; write_header(iterations, dimensions, end[0], end[1], true, kappa, A.rows, A.cols, A.vals); if (0 != check_parameters(A, dimensions, iterations, end[0], end[1], true)) { if (A.rows == 0 || A.cols == 0) { R = new SVDRec(); if (SVDVerbosity > 0) printf("TRANSPOSING THE MATRIX FOR SPEED\n"); transpose = true; A = svdTransposeS(A); machar(); // XXX has side effect, computes machine precision fake_memset_127(bnd); OPBTemp = svd_doubleArray(A.rows, false, "las2: OPBTemp"); steps = lanso(A, iterations, dimensions, end[0], end[1], ritz, bnd, wptr, ref_neig, n); int neig = ref_neig[0]; // XXX unwrap neig printf("NUMBER OF LANCZOS STEPS = %6d\n" + "RITZ VALUES STABILIZED = %6d\n", steps + 1, neig); printf("\nCOMPUTED RITZ VALUES (ERROR BNDS)\n");
s = svd_doubleArray(jsq, true, "ritvec: s"); xv2 = svd_doubleArray(n, false, "ritvec: xv2"); svd_dcopy(js, alf, 1, w1, -1); svd_dcopy(steps, bet, 1, 1, w2, 1, -1); // WAS svd_dcopy(steps, &bet[1], 1, &w2[1], -1); imtql2(js, js, w1, w2, s); if (0 != ierr) return 0; // TODO use exception here? for (i = 0; i < n; i++) w1[i] = 0.0; for (i = 0; i < js; i++) { store(n, storeVals.RETRQ, i, w2); svd_daxpy(n, s[tmp], w2, 1, w1, 1); tmp -= js; rotateArray(R.Vt.value, R.Vt.rows * R.Vt.cols, x * R.Vt.cols); R.d = svd_imin(R.d, nsig); for (x = 0; x < R.d; x++) { svd_opb(A, R.Vt.value[x], xv2, OPBTemp); tmp0 = svd_ddot(n, R.Vt.value[x], 1, xv2, 1); svd_daxpy(n, -tmp0, R.Vt.value[x], 1, xv2, 1); tmp0 = Math.sqrt(tmp0); xnorm = Math.sqrt(svd_ddot(n, xv2, 1, xv2, 1)); svd_opa(A, R.Vt.value[x], R.Ut.value[x]); tmp1 = 1.0 / tmp0; svd_dscal(A.rows, tmp1, R.Ut.value[x], 1);
stpone(A, wptr, ref_rnm, ref_tol, n); double tol = ref_tol[0]; // XXX unwrap double rnm = ref_rnm[0]; // XXX unwrap ll = 0; first = 1; last = svd_imin(dimensions + svd_imax(8, dimensions), iterations); ENOUGH = false; double[] ref2_rnm = new double[] { rnm }; // XXX wrap double[] ref2_tol = new double[] { tol }; // XXX wrap j = lanczos_step(A, first, last, wptr, alf, eta, oldeta, bet, ref_ll, ref_ENOUGH, ref2_rnm, ref2_tol, n); ll = ref_ll[0]; // XXX unwrap svd_dcopy(i-l+1, alf, l, 1, ritz, l, -1); // WAS svd_dcopy(i-l+1, &alf[l], 1, &ritz[l], -1); svd_dcopy(i-l, bet, l+1, 1, wrk, l+1, -1); // WAS svd_dcopy(i-l, &bet[l+1], 1, &wrk[l+1], -1); imtqlb(i-l+1, ritz, wrk, bnd, l); // TODO start at l svd_error("svdLAS2: imtqlb failed to converge (ierr = %ld)\n", ierr); svd_error(" l = %ld i = %ld\n", l, i); for (id3 = l; id3 <= i; id3++) svd_error(" %ld %lg %lg %lg\n", id3, ritz[id3], wrk[id3], bnd[id3]); bnd[id3] = rnm * fabs(bnd[id3]); l = i + 1;
enum storeVals {STORQ, RETRQ, STORP, RETRP};
mid = svd_idamax(step + 1, bnd, 0, 1); if (fabs(ritz[i-1] - ritz[i]) < eps34 * fabs(ritz[i])) if (bnd[i] > tol && bnd[i-1] > tol) { bnd[i-1] = Math.sqrt(bnd[i] * bnd[i] + bnd[i-1] * bnd[i-1]); if (fabs(ritz[i+1] - ritz[i]) < eps34 * fabs(ritz[i])) if (bnd[i] > tol && bnd[i+1] > tol) { bnd[i+1] = Math.sqrt(bnd[i] * bnd[i] + bnd[i+1] * bnd[i+1]); gap = gapl; if (i < step) gapl = ritz[i+1] - ritz[i]; gap = svd_dmin(gap, gapl); if (gap > bnd[i]) bnd[i] = bnd[i] * (bnd[i] / gap); if (bnd[i] <= 16.0 * eps * fabs(ritz[i])) { neig++; if (!enough[0]) enough[0] = endl < ritz[i] && ritz[i] < endr;
/***************************************************************** * Function finds the index of element having max. absolute value* based on * FORTRAN 77 routine from Linpack by J. Dongarra * *****************************************************************/ static int svd_idamax(int n, double[] dx, int ix0, int incx) { int ix,imax; double dmax; if (n < 1) return -1; if (n == 1) return 0; if (incx == 0) return -1; ix = (incx < 0) ? ix0 + ((-n+1) * incx) : ix0; imax = ix; dmax = fabs(dx[ix]); for (int i=1; i < n; i++) { ix += incx; double dtemp = fabs(dx[ix]); if (dtemp > dmax) { dmax = dtemp; imax = ix; } } return imax; }
double[] bndn = new double[n]; System.arraycopy(bnd, offset, bndn, 0, n); imtqlb(n, dn, en, bndn); System.arraycopy(dn, 0, d, offset, n); System.arraycopy(en, 0, e, offset, n);
void stpone(SMat A, double[][] wrkptr, double[] rnmp, double[] tolp, int n) { double t, rnm, anorm; double[] alf = wrkptr[6]; /* get initial vector; default is random */ rnm = startv(A, wrkptr, 0, n); if (rnm == 0.0 || ierr != 0) return; /* normalize starting vector */ t = 1.0 / rnm; svd_datx(n, t, wrkptr[0], 1, wrkptr[1], 1); svd_dscal(n, t, wrkptr[3], 1); /* take the first step */ svd_opb(A, wrkptr[3], wrkptr[0], OPBTemp); alf[0] = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1); svd_daxpy(n, -alf[0], wrkptr[1], 1, wrkptr[0], 1); t = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1); svd_daxpy(n, -t, wrkptr[1], 1, wrkptr[0], 1); alf[0] += t; svd_dcopy(n, wrkptr[0], 1, wrkptr[4], 1); rnm = Math.sqrt(svd_ddot(n, wrkptr[0], 1, wrkptr[4], 1)); anorm = rnm + fabs(alf[0]); rnmp[0] = rnm; tolp[0] = reps * anorm; return; }
if (m == last) break; else { test = fabs(d[m]) + fabs(d[m+1]); if (test + fabs(e[m]) == test) convergence = true; r = svd_pythag(g, 1.0); g = d[m] - p + e[l] / (g + svd_fsign(r, g)); s = 1.0; c = 1.0; f = s * e[i]; b = c * e[i]; r = svd_pythag(f, g); e[i+1] = r; if (r == 0.0) underflow = true;