001package ch.akuhn.matrix.eigenvalues; 002 003import java.util.Arrays; 004 005import org.netlib.arpack.ARPACK; 006import org.netlib.util.doubleW; 007import org.netlib.util.intW; 008 009import ch.akuhn.matrix.Matrix; 010import ch.akuhn.matrix.Vector; 011 012/** 013 * Finds a few eigenvalues of a matrix. 014 * <P> 015 * This class use ARPACK to find a few eigenvalues (λ) and corresponding 016 * eigenvectors (<b>x</b>) for the standard eigenvalue problem: 017 * 018 * <PRE> 019 * <CODE>A<b>x</b> = λ<b>x</b></CODE> 020 * </PRE> 021 * 022 * where <code>A</code> is an <code>n</code> × <code>n</code> real 023 * symmetric matrix. 024 * <P> 025 * The only thing that must be supplied in order to use this class on your 026 * problem is to change the array dimensions appropriately, to specify 027 * <em>which</em> eigenvalues you want to compute and to supply a matrix-vector 028 * product 029 * 030 * <PRE> 031 * <b>w</b> ← A<b>v</b> 032 * </PRE> 033 * 034 * in the {@link #callback(Vector)} method. 035 * <P> 036 * Please refer to the ARPACK guide for further information. 037 * <P> 038 * <b>Example:</b> 039 * 040 * <PRE> 041 * Matrix A = <i>…square matrix…</i>; 042 * Eigenvalues eigen = Eigenvalues.of(A).largest(4); 043 * eigen.run(); 044 * double[] l = eigen.values; 045 * Vector[] x = eigen.vectors; 046 * </PRE> 047 * 048 * @author Adrian Kuhn (Java) based on <CODE>ddsimp.f</CODE> by Richard Lehoucq, 049 * Danny Sorensen, Chao Yang (Fortran) 050 * 051 * @see "http://www.caam.rice.edu/software/ARPACK/UG" 052 * 053 */ 054@SuppressWarnings("javadoc") 055public abstract class FewEigenvalues extends Eigenvalues { 056 057 private enum Which { 058 LA, SA, LM, SM, BE 059 }; 060 061 private Which which; 062 063 public static FewEigenvalues of(final Matrix matrix) { 064 assert matrix.isSquare(); 065 return new FewEigenvalues(matrix.columnCount()) { 066 @Override 067 protected Vector callback(Vector vector) { 068 return matrix.mult(vector); 069 } 070 }; 071 } 072 073 public FewEigenvalues(int n) { 074 super(n); 075 this.greatest(20); 076 } 077 078 private FewEigenvalues which(Which which, int nev) { 079 this.which = which; 080 this.nev = nev < n ? nev : n - 1; 081 return this; 082 } 083 084 /** Compute the largest algebraic eigenvalues. */ 085 @Override 086 public FewEigenvalues largest(int nev0) { 087 return which(Which.LA, nev0); 088 } 089 090 /** Compute the smallest algebraic eigenvalues. */ 091 public FewEigenvalues smallest(int nev0) { 092 return which(Which.SA, nev0); 093 } 094 095 /** Compute the largest eigenvalues in magnitude. */ 096 public FewEigenvalues greatest(int nev0) { 097 return which(Which.LM, nev0); 098 } 099 100 /** Compute the smallest eigenvalues in magnitude. */ 101 public FewEigenvalues lowest(int nev0) { 102 return which(Which.SM, nev0); 103 } 104 105 /** 106 * Compute eigenvalues from both end of the spectrum. When the 107 * <CODE>nev</CODE> is odd, compute one more from the high end than from the 108 * low end. 109 */ 110 public FewEigenvalues fromBothEnds(int nev0) { 111 return which(Which.BE, nev0); 112 } 113 114 /** 115 * Runs the eigenvalue decomposition, using an implicitly restarted Arnoldi 116 * process (IRAP). Please refer to the ARPACK guide for more information. 117 * 118 */ 119 @Override 120 public Eigenvalues run() { 121 final ARPACK arpack = ARPACK.getInstance(); 122 /* 123 * Setting up parameters for DSAUPD call. 124 */ 125 final intW ido = new intW(0); // must start zero 126 final String bmat = "I"; // standard problem 127 final doubleW tol = new doubleW(0); // uses machine precision 128 final intW info = new intW(0); // request random starting vector 129 final double[] resid = new double[n]; // allocate starting vector 130 /* 131 * NVC is the largest number of basis vectors that will be used in the 132 * Implicitly Restarted Arnoldi Process. Work per major iteration is 133 * proportional to N*NCV*NCV. 134 */ 135 final int ncv = Math.min(nev * 4, n); // rule of thumb use twice nev 136 final double[] v = new double[n * ncv]; 137 final double[] workd = new double[3 * n]; 138 final double[] workl = new double[ncv * (ncv + 8)]; 139 // Parameters 140 final int[] iparam = new int[11]; 141 { 142 final int ishfts = 1; // use the exact shift strategy 143 final int maxitr = 300; // max number of Arnoldi iterations 144 final int mode = 1; // use "mode 1" 145 iparam[1 - 1] = ishfts; 146 iparam[3 - 1] = maxitr; 147 iparam[7 - 1] = mode; 148 } 149 final int[] ipntr = new int[11]; 150 // debug setting 151 org.netlib.arpack.arpack_debug.ndigit.val = -3; 152 org.netlib.arpack.arpack_debug.logfil.val = 6; 153 org.netlib.arpack.arpack_debug.msgets.val = 0; 154 org.netlib.arpack.arpack_debug.msaitr.val = 0; 155 org.netlib.arpack.arpack_debug.msapps.val = 0; 156 org.netlib.arpack.arpack_debug.msaupd.val = 0; 157 org.netlib.arpack.arpack_debug.msaup2.val = 0; 158 org.netlib.arpack.arpack_debug.mseigt.val = 0; 159 org.netlib.arpack.arpack_debug.mseupd.val = 0; 160 /* 161 * Main loop (reverse communication loop) 162 * 163 * Repeatedly calls the routine DSAUPD and takes actions indicated by 164 * parameter IDO until either convergence is indicated or maxitr has 165 * been exceeded. 166 */ 167 assert n > 0; 168 assert nev > 0; 169 assert ncv > nev && ncv <= n; 170 171 while (true) { 172 arpack.dsaupd( 173 ido, // reverse communication parameter 174 bmat, // "I" = standard problem 175 n, // problem size 176 which.name(), // which values are requested? 177 nev, // who many values? 178 tol, // 0 = use machine precision 179 resid, 180 ncv, 181 v, 182 n, 183 iparam, 184 ipntr, 185 workd, 186 workl, 187 workl.length, 188 info); 189 if (!(ido.val == 1 || ido.val == -1)) 190 break; 191 /* 192 * Perform matrix vector multiplication 193 * 194 * y <--- OP*x 195 * 196 * The user should supply his own matrix- vector multiplication 197 * routine here that takes workd(ipntr(1)) as the input, and return 198 * the result to workd(ipntr(2)). 199 */ 200 final int x0 = ipntr[1 - 1] - 1; // Fortran is off-by-one compared 201 // to Java! 202 final int y0 = ipntr[2 - 1] - 1; 203 final Vector x = Vector.copy(workd, x0, n); 204 final Vector y = this.callback(x); 205 assert y.size() == n; 206 y.storeOn(workd, y0); 207 } 208 /* 209 * Either we have convergence or there is an error. 210 */ 211 if (info.val != 0) 212 throw new Error("dsaupd ERRNO = " + info.val 213 + ", see http://www.caam.rice.edu/software/ARPACK/UG/node136.html"); 214 /* 215 * Post-Process using DSEUPD. 216 * 217 * Computed eigenvalues may be extracted. 218 * 219 * The routine DSEUPD now called to do this post processing (other modes 220 * may require more complicated post processing than "mode 1"). 221 */ 222 final boolean rvec = true; // request vectors 223 final boolean[] select = new boolean[ncv]; 224 final double[] d = new double[ncv * 2]; 225 final double sigma = 0; // not used in "mode 1" 226 final intW ierr = new intW(0); 227 final intW nevW = new intW(nev); 228 arpack.dseupd( 229 rvec, 230 "All", 231 select, 232 d, 233 v, 234 n, 235 sigma, 236 bmat, 237 n, 238 which.name(), 239 nevW, 240 tol.val, 241 resid, 242 ncv, 243 v, 244 n, 245 iparam, 246 ipntr, 247 workd, 248 workl, 249 workl.length, 250 ierr); 251 if (ierr.val < 0) 252 throw new Error("dseupd ERRNO = " + info.val 253 + ", see http://www.caam.rice.edu/software/ARPACK/UG/node136.html"); 254 /* 255 * Eigenvalues are returned in the first column of the two dimensional 256 * array D and the corresponding eigenvectors are returned in the first 257 * NCONV (=IPARAM(5)) columns of the two dimensional array V if 258 * requested. Otherwise, an orthogonal basis for the invariant subspace 259 * corresponding to the eigenvalues in D is returned in V. 260 */ 261 final int nconv = iparam[5 - 1]; 262 value = Arrays.copyOf(d, nconv); 263 vector = new Vector[nconv]; 264 for (int i = 0; i < value.length; i++) { 265 vector[i] = Vector.copy(v, i * n, n); 266 } 267 return this; 268 } 269 270 protected abstract Vector callback(Vector vector); 271 272}