001// --- BEGIN LICENSE BLOCK ---
002/* 
003 * Copyright (c) 2009-2011, Mikio L. Braun
004 *               2011, Nicolas Oury
005 * All rights reserved.
006 * 
007 * Redistribution and use in source and binary forms, with or without
008 * modification, are permitted provided that the following conditions are
009 * met:
010 * 
011 *     * Redistributions of source code must retain the above copyright
012 *       notice, this list of conditions and the following disclaimer.
013 * 
014 *     * Redistributions in binary form must reproduce the above
015 *       copyright notice, this list of conditions and the following
016 *       disclaimer in the documentation and/or other materials provided
017 *       with the distribution.
018 * 
019 *     * Neither the name of the Technische Universit?t Berlin nor the
020 *       names of its contributors may be used to endorse or promote
021 *       products derived from this software without specific prior
022 *       written permission.
023 * 
024 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
025 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
026 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
027 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
028 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
029 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
030 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
031 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
032 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
033 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
034 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
035 */
036// --- END LICENSE BLOCK ---
037
038package org.jblas;
039
040/**
041 * <p>Eigenvalue and Eigenvector related functions.</p>
042 * <p/>
043 * <p>Methods exist for working with symmetric matrices or general eigenvalues.
044 * The symmetric versions are usually much faster on symmetric matrices.</p>
045 */
046public class Eigen {
047    private static final DoubleMatrix dummyDouble = new DoubleMatrix(1);
048
049    /**
050     * Compute the eigenvalues for a symmetric matrix.
051     */
052    public static DoubleMatrix symmetricEigenvalues(DoubleMatrix A) {
053        A.assertSquare();
054        DoubleMatrix eigenvalues = new DoubleMatrix(A.rows);
055        int isuppz[] = new int[2 * A.rows];
056        SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyDouble, isuppz);
057        return eigenvalues;
058    }
059
060    /**
061     * Computes the eigenvalues and eigenvectors for a symmetric matrix.
062     *
063     * @return an array of DoubleMatrix objects containing the eigenvectors
064     *         stored as the columns of the first matrix, and the eigenvalues as
065     *         diagonal elements of the second matrix.
066     */
067    public static DoubleMatrix[] symmetricEigenvectors(DoubleMatrix A) {
068        A.assertSquare();
069        DoubleMatrix eigenvalues = new DoubleMatrix(A.rows);
070        DoubleMatrix eigenvectors = A.dup();
071        int isuppz[] = new int[2 * A.rows];
072        SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz);
073        return new DoubleMatrix[]{eigenvectors, DoubleMatrix.diag(eigenvalues)};
074    }
075
076    /**
077     * Computes the eigenvalues of a general matrix.
078     */
079    public static ComplexDoubleMatrix eigenvalues(DoubleMatrix A) {
080        A.assertSquare();
081        DoubleMatrix WR = new DoubleMatrix(A.rows);
082        DoubleMatrix WI = WR.dup();
083        SimpleBlas.geev('N', 'N', A.dup(), WR, WI, dummyDouble, dummyDouble);
084
085        return new ComplexDoubleMatrix(WR, WI);
086    }
087
088    /**
089     * Computes the eigenvalues and eigenvectors of a general matrix.
090     *
091     * @return an array of ComplexDoubleMatrix objects containing the eigenvectors
092     *         stored as the columns of the first matrix, and the eigenvalues as the
093     *         diagonal elements of the second matrix.
094     */
095    public static ComplexDoubleMatrix[] eigenvectors(DoubleMatrix A) {
096        A.assertSquare();
097        // setting up result arrays
098        DoubleMatrix WR = new DoubleMatrix(A.rows);
099        DoubleMatrix WI = WR.dup();
100        DoubleMatrix VR = new DoubleMatrix(A.rows, A.rows);
101
102        SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyDouble, VR);
103
104        // transferring the result
105        ComplexDoubleMatrix E = new ComplexDoubleMatrix(WR, WI);
106        ComplexDoubleMatrix V = new ComplexDoubleMatrix(A.rows, A.rows);
107        //System.err.printf("VR = %s\n", VR.toString());
108        for (int i = 0; i < A.rows; i++) {
109            if (E.get(i).isReal()) {
110                V.putColumn(i, new ComplexDoubleMatrix(VR.getColumn(i)));
111            } else {
112                ComplexDoubleMatrix v = new ComplexDoubleMatrix(VR.getColumn(i), VR.getColumn(i + 1));
113                V.putColumn(i, v);
114                V.putColumn(i + 1, v.conji());
115                i += 1;
116            }
117        }
118        return new ComplexDoubleMatrix[]{V, ComplexDoubleMatrix.diag(E)};
119    }
120
121    /**
122     * Compute generalized eigenvalues of the problem A x = L B x.
123     *
124     * @param A symmetric Matrix A. Only the upper triangle will be considered.
125     * @param B symmetric Matrix B. Only the upper triangle will be considered.
126     * @return a vector of eigenvalues L.
127     */
128    public static DoubleMatrix symmetricGeneralizedEigenvalues(DoubleMatrix A, DoubleMatrix B) {
129        A.assertSquare();
130        B.assertSquare();
131        DoubleMatrix W = new DoubleMatrix(A.rows);
132        SimpleBlas.sygvd(1, 'N', 'U', A.dup(), B.dup(), W);
133        return W;
134    }
135
136    /**
137     * Solve a general problem A x = L B x.
138     *
139     * @param A symmetric matrix A
140     * @param B symmetric matrix B
141     * @return an array of matrices of length two. The first one is an array of the eigenvectors X
142     *         The second one is A vector containing the corresponding eigenvalues L.
143     */
144    public static DoubleMatrix[] symmetricGeneralizedEigenvectors(DoubleMatrix A, DoubleMatrix B) {
145        A.assertSquare();
146        B.assertSquare();
147        DoubleMatrix[] result = new DoubleMatrix[2];
148        DoubleMatrix dA = A.dup();
149        DoubleMatrix dB = B.dup();
150        DoubleMatrix W = new DoubleMatrix(dA.rows);
151        SimpleBlas.sygvd(1, 'V', 'U', dA, dB, W);
152        result[0] = dA;
153        result[1] = W;
154        return result;
155    }
156
157//BEGIN
158  // The code below has been automatically generated.
159  // DO NOT EDIT!
160    private static final FloatMatrix dummyFloat = new FloatMatrix(1);
161
162    /**
163     * Compute the eigenvalues for a symmetric matrix.
164     */
165    public static FloatMatrix symmetricEigenvalues(FloatMatrix A) {
166        A.assertSquare();
167        FloatMatrix eigenvalues = new FloatMatrix(A.rows);
168        int isuppz[] = new int[2 * A.rows];
169        SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyFloat, isuppz);
170        return eigenvalues;
171    }
172
173    /**
174     * Computes the eigenvalues and eigenvectors for a symmetric matrix.
175     *
176     * @return an array of FloatMatrix objects containing the eigenvectors
177     *         stored as the columns of the first matrix, and the eigenvalues as
178     *         diagonal elements of the second matrix.
179     */
180    public static FloatMatrix[] symmetricEigenvectors(FloatMatrix A) {
181        A.assertSquare();
182        FloatMatrix eigenvalues = new FloatMatrix(A.rows);
183        FloatMatrix eigenvectors = A.dup();
184        int isuppz[] = new int[2 * A.rows];
185        SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz);
186        return new FloatMatrix[]{eigenvectors, FloatMatrix.diag(eigenvalues)};
187    }
188
189    /**
190     * Computes the eigenvalues of a general matrix.
191     */
192    public static ComplexFloatMatrix eigenvalues(FloatMatrix A) {
193        A.assertSquare();
194        FloatMatrix WR = new FloatMatrix(A.rows);
195        FloatMatrix WI = WR.dup();
196        SimpleBlas.geev('N', 'N', A.dup(), WR, WI, dummyFloat, dummyFloat);
197
198        return new ComplexFloatMatrix(WR, WI);
199    }
200
201    /**
202     * Computes the eigenvalues and eigenvectors of a general matrix.
203     *
204     * @return an array of ComplexFloatMatrix objects containing the eigenvectors
205     *         stored as the columns of the first matrix, and the eigenvalues as the
206     *         diagonal elements of the second matrix.
207     */
208    public static ComplexFloatMatrix[] eigenvectors(FloatMatrix A) {
209        A.assertSquare();
210        // setting up result arrays
211        FloatMatrix WR = new FloatMatrix(A.rows);
212        FloatMatrix WI = WR.dup();
213        FloatMatrix VR = new FloatMatrix(A.rows, A.rows);
214
215        SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyFloat, VR);
216
217        // transferring the result
218        ComplexFloatMatrix E = new ComplexFloatMatrix(WR, WI);
219        ComplexFloatMatrix V = new ComplexFloatMatrix(A.rows, A.rows);
220        //System.err.printf("VR = %s\n", VR.toString());
221        for (int i = 0; i < A.rows; i++) {
222            if (E.get(i).isReal()) {
223                V.putColumn(i, new ComplexFloatMatrix(VR.getColumn(i)));
224            } else {
225                ComplexFloatMatrix v = new ComplexFloatMatrix(VR.getColumn(i), VR.getColumn(i + 1));
226                V.putColumn(i, v);
227                V.putColumn(i + 1, v.conji());
228                i += 1;
229            }
230        }
231        return new ComplexFloatMatrix[]{V, ComplexFloatMatrix.diag(E)};
232    }
233
234    /**
235     * Compute generalized eigenvalues of the problem A x = L B x.
236     *
237     * @param A symmetric Matrix A. Only the upper triangle will be considered.
238     * @param B symmetric Matrix B. Only the upper triangle will be considered.
239     * @return a vector of eigenvalues L.
240     */
241    public static FloatMatrix symmetricGeneralizedEigenvalues(FloatMatrix A, FloatMatrix B) {
242        A.assertSquare();
243        B.assertSquare();
244        FloatMatrix W = new FloatMatrix(A.rows);
245        SimpleBlas.sygvd(1, 'N', 'U', A.dup(), B.dup(), W);
246        return W;
247    }
248
249    /**
250     * Solve a general problem A x = L B x.
251     *
252     * @param A symmetric matrix A
253     * @param B symmetric matrix B
254     * @return an array of matrices of length two. The first one is an array of the eigenvectors X
255     *         The second one is A vector containing the corresponding eigenvalues L.
256     */
257    public static FloatMatrix[] symmetricGeneralizedEigenvectors(FloatMatrix A, FloatMatrix B) {
258        A.assertSquare();
259        B.assertSquare();
260        FloatMatrix[] result = new FloatMatrix[2];
261        FloatMatrix dA = A.dup();
262        FloatMatrix dB = B.dup();
263        FloatMatrix W = new FloatMatrix(dA.rows);
264        SimpleBlas.sygvd(1, 'V', 'U', dA, dB, W);
265        result[0] = dA;
266        result[1] = W;
267        return result;
268    }
269
270//END
271}