AllEigenvalues.java

package ch.akuhn.matrix.eigenvalues;

import java.util.Arrays;

import org.netlib.lapack.LAPACK;
import org.netlib.util.intW;

import ch.akuhn.matrix.Matrix;
import ch.akuhn.matrix.Vector;


/**
 * Finds all eigenvalues of a matrix.
 *<P>
 * Computes for an <code>n</code>×<code>n</code> real nonsymmetric matrix
 * <code>A</code>, the eigenvalues (λ) and, optionally, the left and/or
 * right eigenvectors. The computed eigenvectors are normalized to have
 * Euclidean norm equal to 1 and largest component real.
 *<P>
 * 
 * @author Adrian Kuhn
 * 
 * @see http://www.netlib.org/lapack/double/dgeev.f
 * 
 */
public class AllEigenvalues extends Eigenvalues {

	private LAPACK lapack = LAPACK.getInstance();
	
	private boolean l = true;
	private boolean r = false;

	public AllEigenvalues(Matrix A) {
		super(A.columnCount());
		assert A.isSquare();
		this.A = A;
	}
	
	private Matrix A;

	@Override
	public AllEigenvalues run() {
		double[] wr = new double[n];
		double[] wi = new double[n];
		intW info = new intW(0);
		double[] a = A.asColumnMajorArray();
		double[] vl = new double[l ? n * n : 0];
		double[] vr = new double[r ? n * n : 0];
		double[] work = allocateWorkspace();
		lapack.dgeev(
				jobv(l),
				jobv(r),
				n,
				a, // overwritten on output!
				n,
				wr, // output: real eigenvalues 
				wi, // output: imaginary eigenvalues 
				vl, // output:: left eigenvectors
				n,
				vr, // output:: right eigenvectors
				n,
				work,
				work.length,
				info );
		if (info.val != 0) throw new Error("dgeev ERRNO=" + info.val);
		postprocess(wr, vl);
		return this;
	}

	/** <PRE>[wr,vl.enum_cons(n)]
	 *  .transpose
	 *  .sort_by(&:first)
	 *  .tranpose
	 *  .revert</PRE>
	 * 
	 */
	private void postprocess(double[] wr, double[] vl) {
		class Eigen implements Comparable<Eigen> {
			double value0;
			Vector vector0;
			@Override
			public int compareTo(Eigen eigen) {
				return Double.compare(value0, eigen.value0);
			}
		}
		Eigen[] eigen = new Eigen[n];
		for (int i = 0; i < n; i++) {
			eigen[i] = new Eigen();
			eigen[i].value0 = wr[i];
			eigen[i].vector0 = Vector.copy(vl, i*n, n); 
		}
		Arrays.sort(eigen);
		value = new double[nev];
		vector = new Vector[nev];
		for (int i = 0; i < nev; i++) {
			value[i] = eigen[n-nev+i].value0;
			vector[i] = eigen[n-nev+i].vector0;
		}
	}

	private String jobv(boolean canHasVectors) {
		return canHasVectors ? "V" : "N";
	}

	/** If LWORK = -1, then a workspace query is assumed; the routine
	 * only calculates the optimal size of the WORK array, returns
	 * this value as the first entry of the WORK array.
	 * 
	 */
	private double[] allocateWorkspace() {
		int lwork = ((l || r) ? 4 : 3) * n;
        double[] query = new double[1];
        intW info = new intW(0);
        lapack.dgeev(
				jobv(l),
				jobv(r),
				n,
				null,
				n,
				null,
				null,
				null,
				n,
				null,
				n,
				query,
				-1,
				info );
        if (info.val == 0) lwork = (int) query[0];
        return new double[lwork];
	}
	
}