/***************************************************************** * 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; }
/***************************************************************** * 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; }
enum storeVals {STORQ, RETRQ, STORP, RETRP};
enum storeVals {STORQ, RETRQ, STORP, RETRP};
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 = 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;
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 = 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;
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; }
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; }
kappa = svd_dmax(fabs(kappa), eps34);
kappa = svd_dmax(fabs(kappa), eps34);
bnd[id3] = rnm * fabs(bnd[id3]); l = i + 1;
bnd[id3] = rnm * fabs(bnd[id3]); l = i + 1;