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 &gt; 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 &gt; iu</code> or <code>il &lt; 0 </code> or <code>iu &gt; 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 &gt; 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 &gt; iu</code> or <code>il &lt; 0 </code>
280   *         or <code>iu &gt; 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 &gt; 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 &gt; iu</code> or <code>il &lt; 0 </code> or <code>iu &gt; 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 &gt; 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 &gt; iu</code> or <code>il &lt; 0 </code>
539   *         or <code>iu &gt; 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}