001// --- BEGIN LICENSE BLOCK --- 002/* 003 * Copyright (c) 2009-13, Mikio L. Braun 004 * 2011, Nicola Oury 005 * 2013, Alexander Sehlström 006 * 007 * All rights reserved. 008 * 009 * Redistribution and use in source and binary forms, with or without 010 * modification, are permitted provided that the following conditions are 011 * met: 012 * 013 * * Redistributions of source code must retain the above copyright 014 * notice, this list of conditions and the following disclaimer. 015 * 016 * * Redistributions in binary form must reproduce the above 017 * copyright notice, this list of conditions and the following 018 * disclaimer in the documentation and/or other materials provided 019 * with the distribution. 020 * 021 * * Neither the name of the Technische Universität Berlin nor the 022 * names of its contributors may be used to endorse or promote 023 * products derived from this software without specific prior 024 * written permission. 025 * 026 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 027 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 028 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 029 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 030 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 031 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 032 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 033 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 034 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 035 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 036 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 037 */ 038// --- END LICENSE BLOCK --- 039 040package org.jblas; 041 042import org.jblas.exceptions.NoEigenResultException; 043import org.jblas.ranges.IntervalRange; 044import org.jblas.ranges.Range; 045 046/** 047 * <p>Eigenvalue and Eigenvector related functions.</p> 048 * <p/> 049 * <p>Methods exist for working with symmetric matrices or general eigenvalues. 050 * The symmetric versions are usually much faster on symmetric matrices.</p> 051 */ 052public class Eigen { 053 private static final DoubleMatrix dummyDouble = new DoubleMatrix(1); 054 055 /** 056 * Compute the eigenvalues for a symmetric matrix. 057 */ 058 public static DoubleMatrix symmetricEigenvalues(DoubleMatrix A) { 059 A.assertSquare(); 060 DoubleMatrix eigenvalues = new DoubleMatrix(A.rows); 061 int isuppz[] = new int[2 * A.rows]; 062 SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyDouble, isuppz); 063 return eigenvalues; 064 } 065 066 /** 067 * Computes the eigenvalues and eigenvectors for a symmetric matrix. 068 * 069 * @return an array of DoubleMatrix objects containing the eigenvectors 070 * stored as the columns of the first matrix, and the eigenvalues as 071 * diagonal elements of the second matrix. 072 */ 073 public static DoubleMatrix[] symmetricEigenvectors(DoubleMatrix A) { 074 A.assertSquare(); 075 DoubleMatrix eigenvalues = new DoubleMatrix(A.rows); 076 DoubleMatrix eigenvectors = A.dup(); 077 int isuppz[] = new int[2 * A.rows]; 078 SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz); 079 return new DoubleMatrix[]{eigenvectors, DoubleMatrix.diag(eigenvalues)}; 080 } 081 082 /** 083 * Computes the eigenvalues of a general matrix. 084 */ 085 public static ComplexDoubleMatrix eigenvalues(DoubleMatrix A) { 086 A.assertSquare(); 087 DoubleMatrix WR = new DoubleMatrix(A.rows); 088 DoubleMatrix WI = WR.dup(); 089 SimpleBlas.geev('N', 'N', A.dup(), WR, WI, dummyDouble, dummyDouble); 090 091 return new ComplexDoubleMatrix(WR, WI); 092 } 093 094 /** 095 * Computes the eigenvalues and eigenvectors of a general matrix. 096 * 097 * @return an array of ComplexDoubleMatrix objects containing the eigenvectors 098 * stored as the columns of the first matrix, and the eigenvalues as the 099 * diagonal elements of the second matrix. 100 */ 101 public static ComplexDoubleMatrix[] eigenvectors(DoubleMatrix A) { 102 A.assertSquare(); 103 // setting up result arrays 104 DoubleMatrix WR = new DoubleMatrix(A.rows); 105 DoubleMatrix WI = WR.dup(); 106 DoubleMatrix VR = new DoubleMatrix(A.rows, A.rows); 107 108 SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyDouble, VR); 109 110 // transferring the result 111 ComplexDoubleMatrix E = new ComplexDoubleMatrix(WR, WI); 112 ComplexDoubleMatrix V = new ComplexDoubleMatrix(A.rows, A.rows); 113 //System.err.printf("VR = %s\n", VR.toString()); 114 for (int i = 0; i < A.rows; i++) { 115 if (E.get(i).isReal()) { 116 V.putColumn(i, new ComplexDoubleMatrix(VR.getColumn(i))); 117 } else { 118 ComplexDoubleMatrix v = new ComplexDoubleMatrix(VR.getColumn(i), VR.getColumn(i + 1)); 119 V.putColumn(i, v); 120 V.putColumn(i + 1, v.conji()); 121 i += 1; 122 } 123 } 124 return new ComplexDoubleMatrix[]{V, ComplexDoubleMatrix.diag(E)}; 125 } 126 127 /** 128 * Compute generalized eigenvalues of the problem A x = L B x. 129 * 130 * @param A symmetric Matrix A. Only the upper triangle will be considered. 131 * @param B symmetric Matrix B. Only the upper triangle will be considered. 132 * @return a vector of eigenvalues L. 133 */ 134 public static DoubleMatrix symmetricGeneralizedEigenvalues(DoubleMatrix A, DoubleMatrix B) { 135 A.assertSquare(); 136 B.assertSquare(); 137 DoubleMatrix W = new DoubleMatrix(A.rows); 138 SimpleBlas.sygvd(1, 'N', 'U', A.dup(), B.dup(), W); 139 return W; 140 } 141 142 /** 143 * Solve a general problem A x = L B x. 144 * 145 * @param A symmetric matrix A 146 * @param B symmetric matrix B 147 * @return an array of matrices of length two. The first one is an array of the eigenvectors X 148 * The second one is A vector containing the corresponding eigenvalues L. 149 */ 150 public static DoubleMatrix[] symmetricGeneralizedEigenvectors(DoubleMatrix A, DoubleMatrix B) { 151 A.assertSquare(); 152 B.assertSquare(); 153 DoubleMatrix[] result = new DoubleMatrix[2]; 154 DoubleMatrix dA = A.dup(); 155 DoubleMatrix dB = B.dup(); 156 DoubleMatrix W = new DoubleMatrix(dA.rows); 157 SimpleBlas.sygvd(1, 'V', 'U', dA, dB, W); 158 result[0] = dA; 159 result[1] = W; 160 return result; 161 } 162 163 /** 164 * Computes selected eigenvalues of the real generalized symmetric-definite eigenproblem of the form A x = L B x 165 * or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric and B is also positive definite. 166 * The selection is based on the given range of values for the desired eigenvalues. 167 * <p/> 168 * The range is half open: (vl,vu]. 169 * 170 * @param A symmetric Matrix A. Only the upper triangle will be considered. 171 * @param B symmetric Matrix B. Only the upper triangle will be considered. 172 * @param vl lower bound of the smallest eigenvalue to return 173 * @param vu upper bound of the largest eigenvalue to return 174 * @throws IllegalArgumentException if <code>vl > vu</code> 175 * @throws NoEigenResultException if no eigevalues are found for the selected range: (vl,vu] 176 * @return a vector of eigenvalues L 177 */ 178 public static DoubleMatrix symmetricGeneralizedEigenvalues(DoubleMatrix A, DoubleMatrix B, double vl, double vu) { 179 A.assertSquare(); 180 B.assertSquare(); 181 A.assertSameSize(B); 182 if (vu <= vl) { 183 throw new IllegalArgumentException("Bound exception: make sure vu > vl"); 184 } 185 double abstol = (double) 1e-9; // What is a good tolerance? 186 int[] m = new int[1]; 187 DoubleMatrix W = new DoubleMatrix(A.rows); 188 DoubleMatrix Z = new DoubleMatrix(A.rows, A.rows); 189 SimpleBlas.sygvx(1, 'N', 'V', 'U', A.dup(), B.dup(), vl, vu, 0, 0, abstol, m, W, Z); 190 if (m[0] == 0) { 191 throw new NoEigenResultException("No eigenvalues found for selected range"); 192 } 193 return W.get(new IntervalRange(0, m[0]), 0); 194 } 195 196 /** 197 * Computes selected eigenvalues of the real generalized symmetric-definite eigenproblem of the form A x = L B x 198 * or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric and B is also positive definite. 199 * The selection is based on the given range of indices for the desired eigenvalues. 200 * 201 * @param A symmetric Matrix A. Only the upper triangle will be considered. 202 * @param B symmetric Matrix B. Only the upper triangle will be considered. 203 * @param il lower index (in ascending order) of the smallest eigenvalue to return (index is 0-based) 204 * @param iu upper index (in ascending order) of the largest eigenvalue to return (index is 0-based) 205 * @throws IllegalArgumentException if <code>il > iu</code> or <code>il < 0 </code> or <code>iu > A.rows - 1</code> 206 * @return a vector of eigenvalues L 207 */ 208 public static DoubleMatrix symmetricGeneralizedEigenvalues(DoubleMatrix A, DoubleMatrix B, int il, int iu) { 209 A.assertSquare(); 210 B.assertSquare(); 211 A.assertSameSize(B); 212 if (iu < il) { 213 throw new IllegalArgumentException("Index exception: make sure iu >= il"); 214 } 215 if (il < 0) { 216 throw new IllegalArgumentException("Index exception: make sure il >= 0"); 217 } 218 if (iu > A.rows - 1) { 219 throw new IllegalArgumentException("Index exception: make sure iu <= A.rows - 1"); 220 } 221 double abstol = (double) 1e-9; // What is a good tolerance? 222 int[] m = new int[1]; 223 DoubleMatrix W = new DoubleMatrix(A.rows); 224 DoubleMatrix Z = new DoubleMatrix(A.rows, A.columns); 225 SimpleBlas.sygvx(1, 'N', 'I', 'U', A.dup(), B.dup(), 0.0, 0.0, il + 1, iu + 1, abstol, m, W, Z); 226 return W.get(new IntervalRange(0, m[0]), 0); 227 } 228 229 /** 230 * Computes selected eigenvalues and their corresponding eigenvectors of the real generalized symmetric-definite 231 * eigenproblem of the form A x = L B x or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric 232 * and B is also positive definite. The selection is based on the given range of values for the desired eigenvalues. 233 * 234 * The range is half open: (vl,vu]. 235 * 236 * @param A symmetric Matrix A. Only the upper triangle will be considered. 237 * @param B symmetric Matrix B. Only the upper triangle will be considered. 238 * @param vl lower bound of the smallest eigenvalue to return 239 * @param vu upper bound of the largest eigenvalue to return 240 * @return an array of matrices of length two. The first one is an array of the eigenvectors x. 241 * The second one is a vector containing the corresponding eigenvalues L. 242 * @throws IllegalArgumentException if <code>vl > vu</code> 243 * @throws NoEigenResultException if no eigevalues are found for the selected range: (vl,vu] 244 */ 245 public static DoubleMatrix[] symmetricGeneralizedEigenvectors(DoubleMatrix A, DoubleMatrix B, double vl, double vu) { 246 A.assertSquare(); 247 B.assertSquare(); 248 A.assertSameSize(B); 249 if (vu <= vl) { 250 throw new IllegalArgumentException("Bound exception: make sure vu > vl"); 251 } 252 double abstol = (double) 1e-9; // What is a good tolerance? 253 int[] m = new int[1]; 254 DoubleMatrix W = new DoubleMatrix(A.rows); 255 DoubleMatrix Z = new DoubleMatrix(A.rows, A.columns); 256 SimpleBlas.sygvx(1, 'V', 'V', 'U', A.dup(), B.dup(), vl, vu, 0, 0, abstol, m, W, Z); 257 if (m[0] == 0) { 258 throw new NoEigenResultException("No eigenvalues found for selected range"); 259 } 260 DoubleMatrix[] result = new DoubleMatrix[2]; 261 Range r = new IntervalRange(0, m[0]); 262 result[0] = Z.get(new IntervalRange(0, A.rows), r); 263 result[1] = W.get(r, 0); 264 return result; 265 } 266 267 /** 268 * Computes selected eigenvalues and their corresponding 269 * eigenvectors of the real generalized symmetric-definite 270 * eigenproblem of the form A x = L B x or, equivalently, 271 * (A - L B)x = 0. Here A and B are assumed to be symmetric 272 * and B is also positive definite. The selection is based 273 * on the given range of values for the desired eigenvalues. 274 * 275 * @param A symmetric Matrix A. Only the upper triangle will be considered. 276 * @param B symmetric Matrix B. Only the upper triangle will be considered. 277 * @param il lower index (in ascending order) of the smallest eigenvalue to return (index is 0-based) 278 * @param iu upper index (in ascending order) of the largest eigenvalue to return (index is 0-based) 279 * @throws IllegalArgumentException if <code>il > iu</code> or <code>il < 0 </code> 280 * or <code>iu > A.rows - 1</code> 281 * @return an array of matrices of length two. The first one is an array of the eigenvectors x. 282 * The second one is a vector containing the corresponding eigenvalues L. 283 */ 284 public static DoubleMatrix[] symmetricGeneralizedEigenvectors(DoubleMatrix A, DoubleMatrix B, int il, int iu) { 285 A.assertSquare(); 286 B.assertSquare(); 287 A.assertSameSize(B); 288 if (iu < il) { 289 throw new IllegalArgumentException("Index exception: make sure iu >= il"); 290 } 291 if (il < 0) { 292 throw new IllegalArgumentException("Index exception: make sure il >= 0"); 293 } 294 if (iu > A.rows - 1) { 295 throw new IllegalArgumentException("Index exception: make sure iu <= A.rows - 1"); 296 } 297 double abstol = (double) 1e-9; // What is a good tolerance? 298 int[] m = new int[1]; 299 DoubleMatrix W = new DoubleMatrix(A.rows); 300 DoubleMatrix Z = new DoubleMatrix(A.rows, A.columns); 301 SimpleBlas.sygvx(1, 'V', 'I', 'U', A.dup(), B.dup(), 0, 0, il + 1, iu + 1, abstol, m, W, Z); 302 DoubleMatrix[] result = new DoubleMatrix[2]; 303 Range r = new IntervalRange(0, m[0]); 304 result[0] = Z.get(new IntervalRange(0, A.rows), r); 305 result[1] = W.get(r, 0); 306 return result; 307 } 308 309//BEGIN 310 // The code below has been automatically generated. 311 // DO NOT EDIT! 312 private static final FloatMatrix dummyFloat = new FloatMatrix(1); 313 314 /** 315 * Compute the eigenvalues for a symmetric matrix. 316 */ 317 public static FloatMatrix symmetricEigenvalues(FloatMatrix A) { 318 A.assertSquare(); 319 FloatMatrix eigenvalues = new FloatMatrix(A.rows); 320 int isuppz[] = new int[2 * A.rows]; 321 SimpleBlas.syevr('N', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, dummyFloat, isuppz); 322 return eigenvalues; 323 } 324 325 /** 326 * Computes the eigenvalues and eigenvectors for a symmetric matrix. 327 * 328 * @return an array of FloatMatrix objects containing the eigenvectors 329 * stored as the columns of the first matrix, and the eigenvalues as 330 * diagonal elements of the second matrix. 331 */ 332 public static FloatMatrix[] symmetricEigenvectors(FloatMatrix A) { 333 A.assertSquare(); 334 FloatMatrix eigenvalues = new FloatMatrix(A.rows); 335 FloatMatrix eigenvectors = A.dup(); 336 int isuppz[] = new int[2 * A.rows]; 337 SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz); 338 return new FloatMatrix[]{eigenvectors, FloatMatrix.diag(eigenvalues)}; 339 } 340 341 /** 342 * Computes the eigenvalues of a general matrix. 343 */ 344 public static ComplexFloatMatrix eigenvalues(FloatMatrix A) { 345 A.assertSquare(); 346 FloatMatrix WR = new FloatMatrix(A.rows); 347 FloatMatrix WI = WR.dup(); 348 SimpleBlas.geev('N', 'N', A.dup(), WR, WI, dummyFloat, dummyFloat); 349 350 return new ComplexFloatMatrix(WR, WI); 351 } 352 353 /** 354 * Computes the eigenvalues and eigenvectors of a general matrix. 355 * 356 * @return an array of ComplexFloatMatrix objects containing the eigenvectors 357 * stored as the columns of the first matrix, and the eigenvalues as the 358 * diagonal elements of the second matrix. 359 */ 360 public static ComplexFloatMatrix[] eigenvectors(FloatMatrix A) { 361 A.assertSquare(); 362 // setting up result arrays 363 FloatMatrix WR = new FloatMatrix(A.rows); 364 FloatMatrix WI = WR.dup(); 365 FloatMatrix VR = new FloatMatrix(A.rows, A.rows); 366 367 SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyFloat, VR); 368 369 // transferring the result 370 ComplexFloatMatrix E = new ComplexFloatMatrix(WR, WI); 371 ComplexFloatMatrix V = new ComplexFloatMatrix(A.rows, A.rows); 372 //System.err.printf("VR = %s\n", VR.toString()); 373 for (int i = 0; i < A.rows; i++) { 374 if (E.get(i).isReal()) { 375 V.putColumn(i, new ComplexFloatMatrix(VR.getColumn(i))); 376 } else { 377 ComplexFloatMatrix v = new ComplexFloatMatrix(VR.getColumn(i), VR.getColumn(i + 1)); 378 V.putColumn(i, v); 379 V.putColumn(i + 1, v.conji()); 380 i += 1; 381 } 382 } 383 return new ComplexFloatMatrix[]{V, ComplexFloatMatrix.diag(E)}; 384 } 385 386 /** 387 * Compute generalized eigenvalues of the problem A x = L B x. 388 * 389 * @param A symmetric Matrix A. Only the upper triangle will be considered. 390 * @param B symmetric Matrix B. Only the upper triangle will be considered. 391 * @return a vector of eigenvalues L. 392 */ 393 public static FloatMatrix symmetricGeneralizedEigenvalues(FloatMatrix A, FloatMatrix B) { 394 A.assertSquare(); 395 B.assertSquare(); 396 FloatMatrix W = new FloatMatrix(A.rows); 397 SimpleBlas.sygvd(1, 'N', 'U', A.dup(), B.dup(), W); 398 return W; 399 } 400 401 /** 402 * Solve a general problem A x = L B x. 403 * 404 * @param A symmetric matrix A 405 * @param B symmetric matrix B 406 * @return an array of matrices of length two. The first one is an array of the eigenvectors X 407 * The second one is A vector containing the corresponding eigenvalues L. 408 */ 409 public static FloatMatrix[] symmetricGeneralizedEigenvectors(FloatMatrix A, FloatMatrix B) { 410 A.assertSquare(); 411 B.assertSquare(); 412 FloatMatrix[] result = new FloatMatrix[2]; 413 FloatMatrix dA = A.dup(); 414 FloatMatrix dB = B.dup(); 415 FloatMatrix W = new FloatMatrix(dA.rows); 416 SimpleBlas.sygvd(1, 'V', 'U', dA, dB, W); 417 result[0] = dA; 418 result[1] = W; 419 return result; 420 } 421 422 /** 423 * Computes selected eigenvalues of the real generalized symmetric-definite eigenproblem of the form A x = L B x 424 * or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric and B is also positive definite. 425 * The selection is based on the given range of values for the desired eigenvalues. 426 * <p/> 427 * The range is half open: (vl,vu]. 428 * 429 * @param A symmetric Matrix A. Only the upper triangle will be considered. 430 * @param B symmetric Matrix B. Only the upper triangle will be considered. 431 * @param vl lower bound of the smallest eigenvalue to return 432 * @param vu upper bound of the largest eigenvalue to return 433 * @throws IllegalArgumentException if <code>vl > vu</code> 434 * @throws NoEigenResultException if no eigevalues are found for the selected range: (vl,vu] 435 * @return a vector of eigenvalues L 436 */ 437 public static FloatMatrix symmetricGeneralizedEigenvalues(FloatMatrix A, FloatMatrix B, float vl, float vu) { 438 A.assertSquare(); 439 B.assertSquare(); 440 A.assertSameSize(B); 441 if (vu <= vl) { 442 throw new IllegalArgumentException("Bound exception: make sure vu > vl"); 443 } 444 float abstol = (float) 1e-9; // What is a good tolerance? 445 int[] m = new int[1]; 446 FloatMatrix W = new FloatMatrix(A.rows); 447 FloatMatrix Z = new FloatMatrix(A.rows, A.rows); 448 SimpleBlas.sygvx(1, 'N', 'V', 'U', A.dup(), B.dup(), vl, vu, 0, 0, abstol, m, W, Z); 449 if (m[0] == 0) { 450 throw new NoEigenResultException("No eigenvalues found for selected range"); 451 } 452 return W.get(new IntervalRange(0, m[0]), 0); 453 } 454 455 /** 456 * Computes selected eigenvalues of the real generalized symmetric-definite eigenproblem of the form A x = L B x 457 * or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric and B is also positive definite. 458 * The selection is based on the given range of indices for the desired eigenvalues. 459 * 460 * @param A symmetric Matrix A. Only the upper triangle will be considered. 461 * @param B symmetric Matrix B. Only the upper triangle will be considered. 462 * @param il lower index (in ascending order) of the smallest eigenvalue to return (index is 0-based) 463 * @param iu upper index (in ascending order) of the largest eigenvalue to return (index is 0-based) 464 * @throws IllegalArgumentException if <code>il > iu</code> or <code>il < 0 </code> or <code>iu > A.rows - 1</code> 465 * @return a vector of eigenvalues L 466 */ 467 public static FloatMatrix symmetricGeneralizedEigenvalues(FloatMatrix A, FloatMatrix B, int il, int iu) { 468 A.assertSquare(); 469 B.assertSquare(); 470 A.assertSameSize(B); 471 if (iu < il) { 472 throw new IllegalArgumentException("Index exception: make sure iu >= il"); 473 } 474 if (il < 0) { 475 throw new IllegalArgumentException("Index exception: make sure il >= 0"); 476 } 477 if (iu > A.rows - 1) { 478 throw new IllegalArgumentException("Index exception: make sure iu <= A.rows - 1"); 479 } 480 float abstol = (float) 1e-9; // What is a good tolerance? 481 int[] m = new int[1]; 482 FloatMatrix W = new FloatMatrix(A.rows); 483 FloatMatrix Z = new FloatMatrix(A.rows, A.columns); 484 SimpleBlas.sygvx(1, 'N', 'I', 'U', A.dup(), B.dup(), 0.0f, 0.0f, il + 1, iu + 1, abstol, m, W, Z); 485 return W.get(new IntervalRange(0, m[0]), 0); 486 } 487 488 /** 489 * Computes selected eigenvalues and their corresponding eigenvectors of the real generalized symmetric-definite 490 * eigenproblem of the form A x = L B x or, equivalently, (A - L B)x = 0. Here A and B are assumed to be symmetric 491 * and B is also positive definite. The selection is based on the given range of values for the desired eigenvalues. 492 * 493 * The range is half open: (vl,vu]. 494 * 495 * @param A symmetric Matrix A. Only the upper triangle will be considered. 496 * @param B symmetric Matrix B. Only the upper triangle will be considered. 497 * @param vl lower bound of the smallest eigenvalue to return 498 * @param vu upper bound of the largest eigenvalue to return 499 * @return an array of matrices of length two. The first one is an array of the eigenvectors x. 500 * The second one is a vector containing the corresponding eigenvalues L. 501 * @throws IllegalArgumentException if <code>vl > vu</code> 502 * @throws NoEigenResultException if no eigevalues are found for the selected range: (vl,vu] 503 */ 504 public static FloatMatrix[] symmetricGeneralizedEigenvectors(FloatMatrix A, FloatMatrix B, float vl, float vu) { 505 A.assertSquare(); 506 B.assertSquare(); 507 A.assertSameSize(B); 508 if (vu <= vl) { 509 throw new IllegalArgumentException("Bound exception: make sure vu > vl"); 510 } 511 float abstol = (float) 1e-9; // What is a good tolerance? 512 int[] m = new int[1]; 513 FloatMatrix W = new FloatMatrix(A.rows); 514 FloatMatrix Z = new FloatMatrix(A.rows, A.columns); 515 SimpleBlas.sygvx(1, 'V', 'V', 'U', A.dup(), B.dup(), vl, vu, 0, 0, abstol, m, W, Z); 516 if (m[0] == 0) { 517 throw new NoEigenResultException("No eigenvalues found for selected range"); 518 } 519 FloatMatrix[] result = new FloatMatrix[2]; 520 Range r = new IntervalRange(0, m[0]); 521 result[0] = Z.get(new IntervalRange(0, A.rows), r); 522 result[1] = W.get(r, 0); 523 return result; 524 } 525 526 /** 527 * Computes selected eigenvalues and their corresponding 528 * eigenvectors of the real generalized symmetric-definite 529 * eigenproblem of the form A x = L B x or, equivalently, 530 * (A - L B)x = 0. Here A and B are assumed to be symmetric 531 * and B is also positive definite. The selection is based 532 * on the given range of values for the desired eigenvalues. 533 * 534 * @param A symmetric Matrix A. Only the upper triangle will be considered. 535 * @param B symmetric Matrix B. Only the upper triangle will be considered. 536 * @param il lower index (in ascending order) of the smallest eigenvalue to return (index is 0-based) 537 * @param iu upper index (in ascending order) of the largest eigenvalue to return (index is 0-based) 538 * @throws IllegalArgumentException if <code>il > iu</code> or <code>il < 0 </code> 539 * or <code>iu > A.rows - 1</code> 540 * @return an array of matrices of length two. The first one is an array of the eigenvectors x. 541 * The second one is a vector containing the corresponding eigenvalues L. 542 */ 543 public static FloatMatrix[] symmetricGeneralizedEigenvectors(FloatMatrix A, FloatMatrix B, int il, int iu) { 544 A.assertSquare(); 545 B.assertSquare(); 546 A.assertSameSize(B); 547 if (iu < il) { 548 throw new IllegalArgumentException("Index exception: make sure iu >= il"); 549 } 550 if (il < 0) { 551 throw new IllegalArgumentException("Index exception: make sure il >= 0"); 552 } 553 if (iu > A.rows - 1) { 554 throw new IllegalArgumentException("Index exception: make sure iu <= A.rows - 1"); 555 } 556 float abstol = (float) 1e-9; // What is a good tolerance? 557 int[] m = new int[1]; 558 FloatMatrix W = new FloatMatrix(A.rows); 559 FloatMatrix Z = new FloatMatrix(A.rows, A.columns); 560 SimpleBlas.sygvx(1, 'V', 'I', 'U', A.dup(), B.dup(), 0, 0, il + 1, iu + 1, abstol, m, W, Z); 561 FloatMatrix[] result = new FloatMatrix[2]; 562 Range r = new IntervalRange(0, m[0]); 563 result[0] = Z.get(new IntervalRange(0, A.rows), r); 564 result[1] = W.get(r, 0); 565 return result; 566 } 567 568//END 569}