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 }