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