Dnrm2.java

package org.mklab.sdpj.gpack.blaswrap;

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> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class Dnrm2<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>> {

  /**
   * DNRM2 returns the euclidean norm of a vector via the function name, so that
   * DNRM2 := sqrt( x'*x )
   * 
   * @param size 成分の数
   * @param x x
   * @param incx x成分の指数の増分
   * @return ユークリッドノルム
   */
  RS execute(int size, RS[] x, int incx) {
    final RS unit = x[0].createUnit();

    int pointer_x = -1;

    if (size < 1 || incx < 1) {
      return unit.createZero();
    }

    if (size == 1) {
      return x[1 + pointer_x].abs();
    }

    RS scale = unit.createZero();
    RS ssq = unit.createUnit();

    int i1 = (size - 1) * incx + 1;
    for (int ix = 1; incx < 0 ? ix >= i1 : ix <= i1; ix += incx) {
      if (!x[ix + pointer_x].isZero()) {
        RS d1 = x[ix + pointer_x];
        RS absxi = d1.abs();
        if (scale.isLessThan(absxi)) {
          // Computing 2nd power
          RS d2 = scale.divide(absxi);
          ssq = (ssq.multiply(d2.multiply(d2))).add(1);
          scale = absxi;
        } else {
          // Computing 2nd power
          RS d2 = absxi.divide(scale);
          ssq = ssq.add(d2.multiply(d2));
        }
      }
    }

    return scale.multiply(ssq.sqrt());
  }
}