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 }