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    
038    package 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     */
046    public 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    }