ARPACK arpack = ARPACK.getInstance();
/*
* Setting up parameters for DSAUPD call.
*
*/
intW ido = new intW(0); // must start zero
String bmat = "I"; // standard problem
doubleW tol = new doubleW(0); // uses machine precision
intW info = new intW(0); // request random starting vector
double[] resid = new double[n]; // allocate starting vector
/*
* NVC is the largest number of basis vectors that will
* be used in the Implicitly Restarted Arnoldi Process.
* Work per major iteration is proportional to N*NCV*NCV.
*
*/
int ncv = Math.min(nev * 4, n); // rule of thumb use twice nev
double[] v = new double[n * ncv];
double[] workd = new double[3*n];
double[] workl = new double[ncv*(ncv+8)];
// Parameters
int[] iparam = new int[11];
{
int ishfts = 1; // use the exact shift strategy
int maxitr = 300; // max number of Arnoldi iterations
int mode = 1; // use "mode 1"
iparam[1-1] = ishfts;
iparam[3-1] = maxitr;
iparam[7-1] = mode;
}
int[] ipntr = new int[11];
// debug setting
org.netlib.arpack.arpack_debug.ndigit.val = -3;
org.netlib.arpack.arpack_debug.logfil.val = 6;
org.netlib.arpack.arpack_debug.msgets.val = 0;
org.netlib.arpack.arpack_debug.msaitr.val = 0;
org.netlib.arpack.arpack_debug.msapps.val = 0;
org.netlib.arpack.arpack_debug.msaupd.val = 0;
org.netlib.arpack.arpack_debug.msaup2.val = 0;
org.netlib.arpack.arpack_debug.mseigt.val = 0;
org.netlib.arpack.arpack_debug.mseupd.val = 0;
/*
* Main loop (reverse communication loop)
*
* Repeatedly calls the routine DSAUPD and takes
* actions indicated by parameter IDO until either
* convergence is indicated or maxitr has been
* exceeded.
*
*/
assert n > 0;
assert nev > 0;
assert ncv > nev && ncv <= n;
while (true) {
arpack.dsaupd(
ido, // reverse communication parameter
bmat, // "I" = standard problem
n, // problem size
which.name(), // which values are requested?
nev, // who many values?
tol, // 0 = use machine precision
resid,
ncv,
v,
n,
iparam,
ipntr,
workd,
workl,
workl.length,
info );
if (!(ido.val == 1 || ido.val == -1)) break;
/* Perform matrix vector multiplication
*
* y <--- OP*x
*
* The user should supply his own matrix-
* vector multiplication routine here that
* takes workd(ipntr(1)) as the input, and
* return the result to workd(ipntr(2)).
*
*/
int x0 = ipntr[1-1]-1; // Fortran is off-by-one compared to Java!
int y0 = ipntr[2-1]-1;
Vector x = Vector.copy(workd, x0, n);
Vector y = this.callback(x);
assert y.size() == n;
y.storeOn(workd, y0);
}
/*
* Either we have convergence or there is an error.
*
*/
if (info.val != 0) throw new Error("dsaupd ERRNO = "+info.val
+", see http://www.caam.rice.edu/software/ARPACK/UG/node136.html");
/*
* Post-Process using DSEUPD.
*
* Computed eigenvalues may be extracted.
*
* The routine DSEUPD now called to do this
* post processing (other modes may require
* more complicated post processing than
* "mode 1").
*
*/
boolean rvec = true; // request vectors
boolean[] select = new boolean[ncv];
double[] d = new double[ncv * 2];
double sigma = 0; // not used in "mode 1"
intW ierr = new intW(0);
intW nevW = new intW(nev);
arpack.dseupd(
rvec,
"All",
select,
d,