001package ch.akuhn.matrix.eigenvalues;
002
003import java.util.Arrays;
004
005import org.netlib.lapack.LAPACK;
006import org.netlib.util.intW;
007
008import ch.akuhn.matrix.Matrix;
009import ch.akuhn.matrix.Vector;
010
011/**
012 * Finds all eigenvalues of a matrix.
013 * <P>
014 * Computes for an <code>n</code>&times;<code>n</code> real nonsymmetric matrix
015 * <code>A</code>, the eigenvalues (&lambda;) and, optionally, the left and/or
016 * right eigenvectors. The computed eigenvectors are normalized to have
017 * Euclidean norm equal to 1 and largest component real.
018 * <P>
019 * 
020 * @author Adrian Kuhn
021 * 
022 * @see "http://www.netlib.org/lapack/double/dgeev.f"
023 * 
024 */
025public class AllEigenvalues extends Eigenvalues {
026
027        private LAPACK lapack = LAPACK.getInstance();
028
029        private boolean l = true;
030        private boolean r = false;
031
032        /**
033         * Construct with the given matrix
034         * 
035         * @param A
036         */
037        public AllEigenvalues(Matrix A) {
038                super(A.columnCount());
039                assert A.isSquare();
040                this.A = A;
041        }
042
043        private Matrix A;
044
045        @Override
046        public AllEigenvalues run() {
047                final double[] wr = new double[n];
048                final double[] wi = new double[n];
049                final intW info = new intW(0);
050                final double[] a = A.asColumnMajorArray();
051                final double[] vl = new double[l ? n * n : 0];
052                final double[] vr = new double[r ? n * n : 0];
053                final double[] work = allocateWorkspace();
054                lapack.dgeev(
055                                jobv(l),
056                                jobv(r),
057                                n,
058                                a, // overwritten on output!
059                                n,
060                                wr, // output: real eigenvalues
061                                wi, // output: imaginary eigenvalues
062                                vl, // output:: left eigenvectors
063                                n,
064                                vr, // output:: right eigenvectors
065                                n,
066                                work,
067                                work.length,
068                                info);
069                if (info.val != 0)
070                        throw new Error("dgeev ERRNO=" + info.val);
071                postprocess(wr, vl);
072                return this;
073        }
074
075        /**
076         * <PRE>
077         * [wr,vl.enum_cons(n)]
078         *  .transpose
079         *  .sort_by(&:first)
080         *  .tranpose
081         *  .revert
082         * </PRE>
083         * 
084         */
085        private void postprocess(double[] wr, double[] vl) {
086                class Eigen implements Comparable<Eigen> {
087                        double value0;
088                        Vector vector0;
089
090                        @Override
091                        public int compareTo(Eigen eigen) {
092                                return Double.compare(value0, eigen.value0);
093                        }
094                }
095                final Eigen[] eigen = new Eigen[n];
096                for (int i = 0; i < n; i++) {
097                        eigen[i] = new Eigen();
098                        eigen[i].value0 = wr[i];
099                        eigen[i].vector0 = Vector.copy(vl, i * n, n);
100                }
101                Arrays.sort(eigen);
102                value = new double[nev];
103                vector = new Vector[nev];
104                for (int i = 0; i < nev; i++) {
105                        value[i] = eigen[n - nev + i].value0;
106                        vector[i] = eigen[n - nev + i].vector0;
107                }
108        }
109
110        private String jobv(boolean canHasVectors) {
111                return canHasVectors ? "V" : "N";
112        }
113
114        /**
115         * If LWORK = -1, then a workspace query is assumed; the routine only
116         * calculates the optimal size of the WORK array, returns this value as the
117         * first entry of the WORK array.
118         * 
119         */
120        private double[] allocateWorkspace() {
121                int lwork = ((l || r) ? 4 : 3) * n;
122                final double[] query = new double[1];
123                final intW info = new intW(0);
124                lapack.dgeev(
125                                jobv(l),
126                                jobv(r),
127                                n,
128                                null,
129                                n,
130                                null,
131                                null,
132                                null,
133                                n,
134                                null,
135                                n,
136                                query,
137                                -1,
138                                info);
139                if (info.val == 0)
140                        lwork = (int) query[0];
141                return new double[lwork];
142        }
143
144}