import org.ojalgo.OjAlgoUtils;
import org.ojalgo.matrix.decomposition.Eigenvalue;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.R064Store;
import org.ojalgo.netio.BasicLogger;

/**
 * An example that outlines how to solve Generalised Symmetric Definite Eigenproblems
 *
 * @see https://www.ojalgo.org/2019/08/generalised-eigenvalue-problems/
 */
public class GeneralisedEigenvalueProblem {

    public static void main(final String[] args) {

        BasicLogger.debug();
        BasicLogger.debug(GeneralisedEigenvalueProblem.class);
        BasicLogger.debug(OjAlgoUtils.getTitle());
        BasicLogger.debug(OjAlgoUtils.getDate());
        BasicLogger.debug();

        int dim = 5;

        /*
         * Create 2 random Symmetric Positive Definite matrices
         */
        R064Store mtrxA = R064Store.FACTORY.makeSPD(dim);
        R064Store mtrxB = R064Store.FACTORY.makeSPD(dim);

        /*
         * There are several generalisations. 3 are supported by ojAlgo, specified by the enum:
         * Eigenvalue.Generalisation This factory method returns the most common alternative.
         */
        Eigenvalue.Generalised<Double> generalisedEvD = Eigenvalue.R064.makeGeneralised(mtrxA);
        // Generalisation: [A][V] = [B][V][D]

        // Use 2-args alternative

        generalisedEvD.computeValuesOnly(mtrxA, mtrxB);
        generalisedEvD.decompose(mtrxA, mtrxB); // Use one of these methods, not both

        // Alternatively, prepare and then use the usual 1-arg methods

        generalisedEvD.prepare(mtrxB);

        generalisedEvD.computeValuesOnly(mtrxA); // either
        generalisedEvD.decompose(mtrxA); // or
        // Can be called repeatedly with each "preparation"
        // Useful if the B matrix doesn't change but A does

        MatrixStore<Double> mtrxV = generalisedEvD.getV();
        // Eigenvectors in the columns

        MatrixStore<Double> mtrxD = generalisedEvD.getD();
        // Eigenvalues on the diagonal

        /*
         * Reconstruct the left and right hand sides of the eigenproblem
         */
        MatrixStore<Double> left = mtrxA.multiply(mtrxV);
        MatrixStore<Double> right = mtrxB.multiply(mtrxV).multiply(mtrxD);

        BasicLogger.debug();
        BasicLogger.debug("Generalised Eigenproblem [A][V] = [B][V][D]");
        BasicLogger.debugMatrix("[A]", mtrxA);
        BasicLogger.debugMatrix("[B]", mtrxB);
        BasicLogger.debugMatrix("Eigenvectors - [V]", mtrxV);
        BasicLogger.debugMatrix("Eigenvalues - [D]", mtrxD);
        BasicLogger.debugMatrix("Left - [A][V]", left);
        BasicLogger.debugMatrix("Right - [B][V][D]", right);

    }

}