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

  //  Purpose   
  //  =======   
  //  DGER   performs the rank 1 operation   
  //     A := alpha*x*y' + A,   
  //  where alpha is a scalar, x is an m element vector, y is an n element   
  //  vector and A is an m by n matrix.   
  //  Parameters   
  //  ==========   
  //  M      - INTEGER.   
  //           On entry, M specifies the number of rows of the matrix A.   
  //           M must be at least zero.   
  //           Unchanged on exit.   
  //  N      - INTEGER.   
  //           On entry, N specifies the number of columns of the matrix A.   
  //           N must be at least zero.   
  //           Unchanged on exit.   
  //  ALPHA  - DOUBLE PRECISION.   
  //           On entry, ALPHA specifies the scalar alpha.   
  //           Unchanged on exit.   
  //  X      - DOUBLE PRECISION array of dimension at least   
  //           ( 1 + ( m - 1 )*abs( INCX ) ).   
  //           Before entry, the incremented array X must contain the m   
  //           element vector x.   
  //           Unchanged on exit.   
  //  INCX   - INTEGER.   
  //           On entry, INCX specifies the increment for the elements of   
  //           X. INCX must not be zero.   
  //           Unchanged on exit.   
  //  Y      - DOUBLE PRECISION array of dimension at least   
  //           ( 1 + ( n - 1 )*abs( INCY ) ).   
  //           Before entry, the incremented array Y must contain the n   
  //           element vector y.   
  //           Unchanged on exit.   
  //  INCY   - INTEGER.   
  //           On entry, INCY specifies the increment for the elements of   
  //           Y. INCY must not be zero.   
  //           Unchanged on exit.   
  //  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).   
  //           Before entry, the leading m by n part of the array A must   
  //           contain the matrix of coefficients. On exit, A is   
  //           overwritten by the updated matrix.   
  //  LDA    - INTEGER.   
  //           On entry, LDA specifies the first dimension of A as declared   
  //           in the calling (sub) program. LDA must be at least   
  //           max( 1, m ).   
  //           Unchanged on exit.   
  //  Level 2 Blas routine.   
  //  -- Written on 22-October-1986.   
  //     Jack Dongarra, Argonne National Lab.   
  //     Jeremy Du Croz, Nag Central Office.   
  //     Sven Hammarling, Nag Central Office.   
  //     Richard Hanson, Sandia National Labs.   
  //     Test the input parameters.   
  //     Parameter adjustments

  /**
   * @param m m
   * @param n n
   * @param alpha alpha
   * @param x x
   * @param incx incx
   * @param y y
   * @param incy incy
   * @param a a
   * @param lda lda
   * @return result
   */
  public int execute(int m, int n, RS alpha, RS[] x, int incx, RS[] y, int incy, RS[] a, int lda) {
    int pointer_x = -1;
    int pointer_y = -1;

    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    int info = 0;
    if (m < 0) {
      info = 1;
    } else if (n < 0) {
      info = 2;
    } else if (incx == 0) {
      info = 5;
    } else if (incy == 0) {
      info = 7;
    } else if (lda < Math.max(1, m)) {
      info = 9;
    }
    if (info != 0) {
      BLAS.xerbla("DGER  ", info); //$NON-NLS-1$
      return 0;
    }
    /* Quick return if possible. */
    if (m == 0 || n == 0 || alpha.isZero()) {
      return 0;
    }
    /*
     * Start the operations. In this version the elements of A are accessed
     * sequentially with one pass through A.
     */
    int jy;
    if (incy > 0) {
      jy = 1;
    } else {
      jy = 1 - (n - 1) * incy;
    }
    if (incx == 1) {
      for (int j = 1; j <= n; ++j) {
        if (!y[jy + pointer_y].isZero()) {
          RS temp = alpha.multiply(y[jy + pointer_y]);
          for (int i = 1; i <= m; ++i) {
            int p = j * a_dim1 + i + pointer_a;
            a[p] = a[p].add(x[i + pointer_x].multiply(temp));
          }
        }
        jy += incy;
      }
    } else {
      int kx;
      if (incx > 0) {
        kx = 1;
      } else {
        kx = 1 - (m - 1) * incx;
      }
      for (int j = 1; j <= n; ++j) {
        if (!y[jy + pointer_y].isZero()) {
          RS temp = alpha.multiply(y[jy + pointer_y]);
          int ix = kx;
          for (int i = 1; i <= m; ++i) {
            int p = j * a_dim1 + i + pointer_a;
            a[p] = a[p].add(x[ix + pointer_x].multiply(temp));
            ix += incx;
          }
        }
        jy += incy;
      }
    }
    return 0;
  }
}