Dlarfg.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.gpack.blaswrap.BLAS;
import org.mklab.sdpj.gpack.f2clibs.LibF77;


/**
 * @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 Dlarfg<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   
  September 30, 1994   


  Purpose   
  =======   

  DLARFG generates a real elementary reflector H of order n, such   
  that   

     H * ( alpha ) = ( beta ),   H' * H = I.   
         (   x   )   (   0  )   

  where alpha and beta are scalars, and x is an (n-1)-element real   
  vector. H is represented in the form   

     H = I - tau * ( 1 ) * ( 1 v' ) ,   
                   ( v )   

  where tau is a real scalar and v is a real (n-1)-element   
  vector.   

  If the elements of x are all zero, then tau = 0 and H is taken to be   
  the unit matrix.   

  Otherwise  1 <= tau <= 2.   

  Arguments   
  =========   

  N       (input) INTEGER   
       The order of the elementary reflector.   

  ALPHA   (input/output) DOUBLE PRECISION   
       On entry, the value alpha.   
       On exit, it is overwritten with the value beta.   

  X       (input/output) DOUBLE PRECISION array, dimension   
                      (1+(N-2)*abs(INCX))   
       On entry, the vector x.   
       On exit, it is overwritten with the vector v.   

  INCX    (input) INTEGER   
       The increment between elements of X. INCX > 0.   

  TAU     (output) DOUBLE PRECISION   
       The value tau.   

  =====================================================================   
  */
  /**
   * @param n n
   * @param alpha alpha
   * @param x x
   * @param incx incx
   * @param tau tau
   * @return result
   */
  RS[] execute(int n, RS alpha, RS[] x, int incx, RS tau) {
    final RS unit = alpha.createUnit();

    RS tauTemp = tau.createZero();
    RS alphaTemp = alpha;

    if (n <= 1) {
      RS[] ans = unit.createArray(2);
      ans[0] = alphaTemp;
      ans[1] = tauTemp;
      return ans;
    }

    int pointer_x = -1;
    RS[] x1 = unit.createArray(x.length - (1 + pointer_x));
    System.arraycopy(x, 1 + pointer_x, x1, 0, x1.length);
    RS xnorm = BLAS.dnrm2(n - 1, x1, incx);

    if (xnorm.isZero()) {
      /* H  =  I */
      tauTemp = unit.createZero();
    } else {
      /* general case */
      RS d1 = Clapack.dlapy2(alphaTemp, xnorm);
      RS beta = LibF77.d_sign(d1, alphaTemp).unaryMinus();
      RS rs = d1.createUnit();
      RS safmin = Clapack.<RS,RM,CS,CM> dlamch("S",rs).divide(Clapack.<RS,RM,CS,CM> dlamch("E",rs)); //$NON-NLS-1$ //$NON-NLS-2$
      if (beta.abs().isLessThan(safmin)) {
        /* XNORM, BETA may be inaccurate; scale X and recompute them */
        RS rsafmn = unit.createUnit().divide(safmin);

        int knt = 0;
        do {
          ++knt;
          System.arraycopy(x, 1 + pointer_x, x1, 0, x1.length);
          BLAS.dscal(n - 1, rsafmn, x1, incx);
          System.arraycopy(x1, 0, x, 1 + pointer_x, x1.length);
          beta = beta.multiply(rsafmn);
          alphaTemp = alphaTemp.multiply(rsafmn);
        } while (beta.abs().isLessThan(safmin));

        /* New BETA is at most 1, at least SAFMIN */
        System.arraycopy(x, 1 + pointer_x, x1, 0, x1.length);
        xnorm = BLAS.dnrm2(n - 1, x1, incx);
        RS d__1 = Clapack.dlapy2(alphaTemp, xnorm);
        beta = LibF77.d_sign(d__1, alphaTemp).unaryMinus();
        tauTemp = (beta.subtract(alphaTemp)).divide(beta);

        RS d__2 = unit.create(1).divide(alphaTemp.subtract(beta));
        System.arraycopy(x, 1 + pointer_x, x1, 0, x1.length);
        BLAS.dscal(n - 1, d__2, x1, incx);
        System.arraycopy(x1, 0, x, 1 + pointer_x, x1.length);

        /* If ALPHA is subnormal, it may lose relative accuracy */
        alphaTemp = beta;
        for (int j = 1; j <= knt; ++j) {
          alphaTemp = alphaTemp.multiply(safmin);
          /* L20: */
        }
      } else {
        tauTemp = (beta.subtract(alphaTemp)).divide(beta);
        RS d2 = unit.create(1).divide(alphaTemp.subtract(beta));
        System.arraycopy(x, 1 + pointer_x, x1, 0, x1.length);
        BLAS.dscal(n - 1, d2, x1, incx);
        System.arraycopy(x1, 0, x, 1 + pointer_x, x1.length);
        alphaTemp = beta;
      }
    }

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