Dlartg.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;
import org.mklab.sdpj.tool.Tools;


/**
 * -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of
 * California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and
 * Rice University September 30, 1994 Purpose ======= DLARTG generate a plane
 * rotation so that [ CS SN ].[ F ] = [ R ] where CS**2 + SN**2 = 1. [ -SN CS ]
 * [ G ] [ 0 ] This is a slower, more accurate version of the BLAS1 routine
 * DROTG, with the following other differences: F and G are unchanged on return.
 * If G=0, then CS=1 and SN=0. If F=0 and (G .ne. 0), then CS=0 and SN=1 without
 * doing any floating point operations (saves work in DBDSQR when there are
 * zeros on the diagonal). If F exceeds G in magnitude, CS will be positive.
 * Arguments ========= F (input) DOUBLE PRECISION The first component of vector
 * to be rotated. G (input) DOUBLE PRECISION The second component of vector to
 * be rotated. CS (output) DOUBLE PRECISION The cosine of the rotation. SN
 * (output) DOUBLE PRECISION The sine of the rotation. R (output) DOUBLE
 * PRECISION The nonzero component of the rotated vector.
 * =====================================================================
 * 
   * @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 Dlartg<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>> {

  /** */
  private RS safmn2;
  /** */
  private RS safmx2;
  /** */
  private RS safmin;
  /** */
  private RS eps;

  /**
   * 新しく生成された{@link Dlartg}オブジェクトを初期化します。
   * @param rs RS
   */
  public Dlartg(RS rs) {
    RS unit = rs.createUnit();
    this.safmin = Clapack.<RS,RM,CS,CM> dlamch("S",unit); //$NON-NLS-1$
    this.eps = Clapack.<RS,RM,CS,CM> dlamch("E",unit); //$NON-NLS-1$
    RS d1 = Clapack.<RS,RM,CS,CM> dlamch("B",unit); //$NON-NLS-1$
    int i1 = Tools.NumericalScalarToInt((this.safmin.divide(this.eps).log().divide(d1.log()).divide(2)));
    // i__1 = (int)(log(safmin / eps) / log(dlamch_("B")) / 2.);
    this.safmn2 = d1.power(i1);
    this.safmx2 = this.safmn2.inverse();
  }

  /**
   * @param f f
   * @param g g
   * @return result
   */
  public RS[] execute(RS f, RS g) {
    final RS unit = f.createUnit();

    RS csTemp;
    RS snTemp;
    RS rTemp;

    if (g.isZero()) {
      csTemp = unit.create(1);
      snTemp = unit.create(0);
      rTemp = f;
    } else if (f.isZero()) {
      csTemp = unit.create(0);
      snTemp = unit.create(1);
      rTemp = g;
    } else {
      RS f1 = f;
      RS g1 = g;
      /* Computing MAX */
      RS d1 = f1.abs();
      RS d2 = g1.abs();
      RS scale = d1.max(d2);

      if (scale.isGreaterThanOrEquals(this.safmx2)) {
        int count = 0;
        //L10:
        do {
          ++count;
          f1 = f1.multiply(this.safmn2);
          g1 = g1.multiply(this.safmn2);
          /* Computing MAX */
          d1 = f1.abs();
          d2 = g1.abs();
          scale = d1.max(d2);
        } while (scale.isGreaterThan(this.safmx2));

        d1 = f1;
        d2 = g1;
        rTemp = ((d1.multiply(d1).add(d2.multiply(d2)))).sqrt();
        csTemp = f1.divide(rTemp);
        snTemp = g1.divide(rTemp);
        for (int i = 1; i <= count; ++i) {
          rTemp = rTemp.multiply(this.safmx2);
          /* L20: */
        }
      } else if (scale.isLessThanOrEquals(this.safmn2)) {
        int count = 0;
        //L30:
        do {
          ++count;
          f1 = f1.multiply(this.safmx2);
          g1 = g1.multiply(this.safmx2);
          d1 = f1.abs();
          d2 = g1.abs();
          scale = d1.max(d2);
        } while (scale.isLessThanOrEquals(this.safmn2));

        d1 = f1;
        d2 = g1;
        rTemp = ((d1.multiply(d1)).add(d2.multiply(d2))).sqrt();
        csTemp = f1.divide(rTemp);
        snTemp = g1.divide(rTemp);
        for (int i = 1; i <= count; ++i) {
          rTemp = rTemp.multiply(this.safmn2);
          /* L40: */
        }
      } else {
        d1 = f1;
        d2 = g1;
        rTemp = ((d1.multiply(d1)).add(d2.multiply(d2))).sqrt();
        csTemp = f1.divide(rTemp);
        snTemp = g1.divide(rTemp);
      }
      if ((f.abs()).isGreaterThan(g.abs()) && csTemp.isLessThan(0)) {
        csTemp = csTemp.unaryMinus();
        snTemp = snTemp.unaryMinus();
        rTemp = rTemp.unaryMinus();
      }
    }

    RS[] ans = unit.createArray(3);
    ans[0] = csTemp;
    ans[1] = snTemp;
    ans[2] = rTemp;
    return ans;
  }
}