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 (&lambda;) and corresponding
016 * eigenvectors (<b>x</b>) for the standard eigenvalue problem:
017 * 
018 * <PRE>
019 * <CODE>A<b>x</b> = &lambda;<b>x</b></CODE>
020 * </PRE>
021 * 
022 * where <code>A</code> is an <code>n</code> &times; <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> &larr; 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>&hellip;square matrix&hellip;</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}