Dlae2.java

package org.mklab.sdpj.gpack.lapackwrap;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;


/**
 * @author koga
 * @version $Revision$, 2009/04/24
   * @param <RS> type of real scalar
   * @param <RM> type of real matrix
   * @param <CS> type of complex scalar
   * @param <CM> type of complex Matrix
 */
public class Dlae2<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {

  /**
   * -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of
   * California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and
   * Rice University October 31, 1992
   * <p>
   * <p>
   * Purpose ======= DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
   * [ A B ] [ B C ]. On return, RT1 is the eigenvalue of larger absolute value,
   * and RT2 is the eigenvalue of smaller absolute value. Arguments ========= A
   * (input) DOUBLE PRECISION The (1,1) element of the 2-by-2 matrix. B (input)
   * DOUBLE PRECISION The (1,2) and (2,1) elements of the 2-by-2 matrix. C
   * (input) DOUBLE PRECISION The (2,2) element of the 2-by-2 matrix. RT1
   * (output) DOUBLE PRECISION The eigenvalue of larger absolute value. RT2
   * (output) DOUBLE PRECISION The eigenvalue of smaller absolute value. Further
   * Details =============== RT1 is accurate to a few ulps barring
   * over/underflow. RT2 may be inaccurate if there is massive cancellation in
   * the determinant A*C-B*B; higher precision or correctly rounded or correctly
   * truncated arithmetic would be needed to compute RT2 accurately in all
   * cases. Overflow is possible only if RT1 is within a factor of 5 of
   * overflow. Underflow is harmless if the input data is 0 or exceeds
   * underflow_threshold / macheps.
   * =====================================================================
   */
  /**
   * @param a a
   * @param b b
   * @param c c
   * @return {rt1, rt2}
   */
  RS[] execute(RS a, RS b, RS c) {
    final RS unit = a.createUnit();

    RS sm = a.add(c);
    RS df = a.subtract(c);
    RS adf = df.abs();
    RS tb = b.add(b);
    RS ab = tb.abs();

    RS acmn;
    RS acmx;
    if ((a.abs()).isGreaterThan(c.abs())) {
      acmx = a;
      acmn = c;
    } else {
      acmx = c;
      acmn = a;
    }

    RS rt;
    if (adf.isGreaterThan(ab)) {
      RS d1 = ab.divide(adf);
      rt = adf.multiply(((d1.multiply(d1)).add(1)).sqrt());
    } else if (adf.isLessThan(ab)) {
      RS d1 = adf.divide(ab);
      rt = ab.multiply(d1.multiply(d1).add(1).sqrt());
    } else {
      // Includes case AB=ADF=0 
      rt = ab.multiply(unit.create(2).sqrt());
    }

    RS rt1Temp;
    RS rt2Temp;

    if (sm.isLessThan(0)) {
      rt1Temp = (sm.subtract(rt)).divide(2);

      /* Order of execution important.   
       * To get fully accurate smaller eigenvalue,   
       * next line needs to be executed in higher precision.
       */
      rt2Temp = (acmx.divide(rt1Temp).multiply(acmn)).subtract(b.divide(rt1Temp).multiply(b));
    } else if (sm.isGreaterThan(0)) {
      rt1Temp = (sm.add(rt)).divide(2);

      /* Order of execution important.   
       * To get fully accurate smaller eigenvalue,   
       * next line needs to be executed in higher precision.
       */
      rt2Temp = (acmx.divide(rt1Temp).multiply(acmn)).subtract(b.divide(rt1Temp).multiply(b));
    } else {
      // Includes case RT1 = RT2 = 0 
      rt1Temp = rt.divide(2);
      rt2Temp = rt.divide(2).unaryMinus();
    }

    RS[] ans = unit.createArray(2);
    ans[0] = rt1Temp;
    ans[1] = rt2Temp;
    return ans;
  }
}