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.la;
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    //BEGIN
116      // The code below has been automatically generated.
117      // DO NOT EDIT!
118            private static final FloatMatrix dummyFloat = new FloatMatrix(1);
119        
120            /** Compute the eigenvalues for a symmetric matrix. */
121            public static FloatMatrix symmetricEigenvalues(FloatMatrix A) {
122                    A.assertSquare();
123                    FloatMatrix eigenvalues = new FloatMatrix(A.rows);
124            int isuppz[] = new int[2*A.rows];
125                    SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyFloat, isuppz);
126                    return eigenvalues;
127            }
128    
129            /** 
130             * Computes the eigenvalues and eigenvectors for a symmetric matrix. 
131             *
132             * @return an array of FloatMatrix objects containing the eigenvectors 
133             * stored as the columns of the first matrix, and the eigenvalues as
134             * diagonal elements of the second matrix.
135             */
136            public static FloatMatrix[] symmetricEigenvectors(FloatMatrix A) {
137                    A.assertSquare();
138                    FloatMatrix eigenvalues = new FloatMatrix(A.rows);
139                    FloatMatrix eigenvectors = A.dup();
140            int isuppz[] = new int[2*A.rows];
141            SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz);
142                    return new FloatMatrix[] { eigenvectors, FloatMatrix.diag(eigenvalues) };
143            }
144            
145            /** Computes the eigenvalues of a general matrix. */
146            public static ComplexFloatMatrix eigenvalues(FloatMatrix A) {
147                A.assertSquare();
148                FloatMatrix WR = new FloatMatrix(A.rows);
149                FloatMatrix WI = WR.dup();
150                SimpleBlas.geev('N', 'N', A.dup(), WR, WI, dummyFloat, dummyFloat);
151                        
152                return new ComplexFloatMatrix(WR, WI);
153            }
154            
155            /** 
156             * Computes the eigenvalues and eigenvectors of a general matrix.
157             * 
158             * @return an array of ComplexFloatMatrix objects containing the eigenvectors
159             * stored as the columns of the first matrix, and the eigenvalues as the
160             * diagonal elements of the second matrix.
161             */
162            public static ComplexFloatMatrix[] eigenvectors(FloatMatrix A) {
163                A.assertSquare();
164                // setting up result arrays
165                FloatMatrix WR = new FloatMatrix(A.rows);
166                FloatMatrix WI = WR.dup();
167                FloatMatrix VR = new FloatMatrix(A.rows, A.rows);
168    
169                SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyFloat, VR);
170                
171                // transferring the result
172                ComplexFloatMatrix E = new ComplexFloatMatrix(WR, WI);
173                ComplexFloatMatrix V = new ComplexFloatMatrix(A.rows, A.rows);
174                System.err.printf("VR = %s\n", VR.toString());
175                for (int i = 0; i < A.rows; i++) {
176                    if (E.get(i).isReal()) {
177                        V.putColumn(i, new ComplexFloatMatrix(VR.getColumn(i)));
178                    } else {
179                        ComplexFloatMatrix v = new ComplexFloatMatrix(VR.getColumn(i), VR.getColumn(i+1));
180                        V.putColumn(i, v);
181                        V.putColumn(i+1, v.conji());
182                        i += 1;
183                    }
184                }
185                return new ComplexFloatMatrix[] { V, ComplexFloatMatrix.diag(E) };
186            }
187    //END
188    }