IEEEck.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;


//TODO Nan の確認。-inf +inf での足し算とかでエラー出ないかとか。
/**
 * @author koga
 * @version $Revision$, 2009/04/25
   * @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 IEEEck<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   
  June 30, 1998   


  Purpose   
  =======   

  IEEECK is called from the ILAENV to verify that Infinity and   
  possibly NaN arithmetic is safe (i.e. will not trap).   

  Arguments   
  =========   

  ISPEC   (input) INTEGER   
       Specifies whether to test just for inifinity arithmetic   
       or whether to test for infinity and NaN arithmetic.   
       = 0: Verify infinity arithmetic only.   
       = 1: Verify infinity and NaN arithmetic.   

  ZERO    (input) REAL   
       Must contain the value 0.0   
       This is passed to prevent the compiler from optimizing   
       away this code.   

  ONE     (input) REAL   
       Must contain the value 1.0   
       This is passed to prevent the compiler from optimizing   
       away this code.   

  RETURN VALUE:  INTEGER   
       = 0:  Arithmetic failed to produce the correct answers   
       = 1:  Arithmetic produced the correct answers */

  /**
   * @param ispec ispec
   * @param zero zero
   * @param one one
   * @return result
   */
  int execute(int ispec, RS zero, RS one) {
    int ret_val = 1;

    RS posinf = one.divide(zero);
    if (posinf.isLessThanOrEquals(one)) {
      ret_val = 0;
      return ret_val;
    }

    RS neginf = (one.unaryMinus()).divide(zero);
    if (neginf.isGreaterThanOrEquals(zero)) {
      ret_val = 0;
      return ret_val;
    }

    RS negzro = one.divide((neginf.add(one)));
    if (!negzro.equals(zero)) {
      ret_val = 0;
      return ret_val;
    }

    neginf = one.divide(negzro);
    if (neginf.isGreaterThanOrEquals(zero)) {
      ret_val = 0;
      return ret_val;
    }

    RS newzro = negzro.add(zero);
    if (!newzro.equals(zero)) {
      ret_val = 0;
      return ret_val;
    }

    posinf = one.divide(newzro);
    if (posinf.isLessThanOrEquals(one)) {
      ret_val = 0;
      return ret_val;
    }

    neginf = neginf.multiply(posinf);
    if (neginf.isGreaterThanOrEquals(zero)) {
      ret_val = 0;
      return ret_val;
    }

    posinf = posinf.multiply(posinf);
    if (posinf.isLessThanOrEquals(one)) {
      ret_val = 0;
      return ret_val;
    }

    /*     Return if we were only asked to check infinity arithmetic */

    if (ispec == 0) {
      return ret_val;
    }

    RS nan1 = posinf.add(neginf);
    RS nan2 = posinf.divide(neginf);
    RS nan3 = posinf.divide(posinf);
    RS nan4 = posinf.multiply(zero);
    RS nan5 = neginf.multiply(negzro);
    RS nan6 = nan5.multiply(0);

    if (nan1 == nan1) {
      ret_val = 0;
      return ret_val;
    }

    if (nan2 == nan2) {
      ret_val = 0;
      return ret_val;
    }

    if (nan3 == nan3) {
      ret_val = 0;
      return ret_val;
    }

    if (nan4 == nan4) {
      ret_val = 0;
      return ret_val;
    }

    if (nan5 == nan5) {
      ret_val = 0;
      return ret_val;
    }

    if (nan6 == nan6) {
      ret_val = 0;
      return ret_val;
    }

    return ret_val;
  }
}