Dsyr2.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 Dsyr2<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   
  =======   
  DSYR2  performs the symmetric rank 2 operation   
     A := alpha*x*y' + alpha*y*x' + A,   
  where alpha is a scalar, x and y are n element vectors and A is an n   
  by n symmetric matrix.   
  Parameters   
  ==========   
  UPLO   - CHARACTER*1.   
           On entry, UPLO specifies whether the upper or lower   
           triangular part of the array A is to be referenced as   
           follows:   
              UPLO = 'U' or 'u'   Only the upper triangular part of A   
                                  is to be referenced.   
              UPLO = 'L' or 'l'   Only the lower triangular part of A   
                                  is to be referenced.   
           Unchanged on exit.   
  N      - INTEGER.   
           On entry, N specifies the order 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 + ( n - 1 )*abs( INCX ) ).   
           Before entry, the incremented array X must contain the n   
           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 with  UPLO = 'U' or 'u', the leading n by n   
           upper triangular part of the array A must contain the upper   
           triangular part of the symmetric matrix and the strictly   
           lower triangular part of A is not referenced. On exit, the   
           upper triangular part of the array A is overwritten by the   
           upper triangular part of the updated matrix.   
           Before entry with UPLO = 'L' or 'l', the leading n by n   
           lower triangular part of the array A must contain the lower   
           triangular part of the symmetric matrix and the strictly   
           upper triangular part of A is not referenced. On exit, the   
           lower triangular part of the array A is overwritten by the   
           lower triangular part of 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, n ).   
           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 uplo uplo
   * @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
   */
  int execute(String uplo, int n, RS alpha, RS[] x, int incx, RS[] y, int incy, RS[] a, int lda) {
    int a_dim1 = lda;
    int a_offset = 1 + a_dim1 * 1;
    int pointer_a = -a_offset;

    int info = 0;
    if (!BLAS.lsame(uplo, "U") && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$ //$NON-NLS-2$
      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, n)) {
      info = 9;
    }
    if (info != 0) {
      BLAS.xerbla("DSYR2 ", info); //$NON-NLS-1$
      return 0;
    }
    // Quick return if possible.
    if (n == 0 || alpha.isZero()) {
      return 0;
    }
    // Set up the start points in X and Y if the increments are not both unity.
    int kx = 0;
    int ky = 0;
    int jx = 0;
    int jy = 0;
    if (incx != 1 || incy != 1) {
      if (incx > 0) {
        kx = 1;
      } else {
        kx = 1 - (n - 1) * incx;
      }
      if (incy > 0) {
        ky = 1;
      } else {
        ky = 1 - (n - 1) * incy;
      }
      jx = kx;
      jy = ky;
    }
    /* Start the operations. In this version the elements of A are   
     * accessed sequentially with one pass through the triangular part of A.
     */
    if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
      // Form  A  when A is stored in the upper triangle.
      if (incx == 1 && incy == 1) {
        for (int j = 1; j <= n; ++j) {
          if (!x[j - 1].isZero() || !y[j - 1].isZero()) {
            RS temp1 = alpha.multiply(y[j - 1]);
            RS temp2 = alpha.multiply(x[j - 1]);
            for (int i = 1; i <= j; ++i) {
              int p = j * a_dim1 + i + pointer_a;
              a[p] = a[p].add(x[i - 1].multiply(temp1)).add(y[i - 1].multiply(temp2));
            }
          }
        }
      } else {
        for (int j = 1; j <= n; ++j) {
          if (!x[jx - 1].isZero() || !y[jy - 1].isZero()) {
            RS temp1 = alpha.multiply(y[jy - 1]);
            RS temp2 = alpha.multiply(x[jx - 1]);
            int ix = kx;
            int iy = ky;
            for (int i = 1; i <= j; ++i) {
              int p = j * a_dim1 + i + pointer_a;
              a[p] = a[p].add(x[ix - 1].multiply(temp1)).add(y[iy - 1].multiply(temp2));
              ix += incx;
              iy += incy;
            }
          }
          jx += incx;
          jy += incy;
        }
      }
    } else {
      // Form  A  when A is stored in the lower triangle.
      if (incx == 1 && incy == 1) {
        for (int j = 1; j <= n; ++j) {
          if (!x[j - 1].isZero() || !y[j - 1].isZero()) {
            RS temp1 = alpha.multiply(y[j - 1]);
            RS temp2 = alpha.multiply(x[j - 1]);
            for (int i = j; i <= n; ++i) {
              int p = j * a_dim1 + i + pointer_a;
              a[p] = a[p].add(x[i - 1].multiply(temp1)).add(y[i - 1].multiply(temp2));
            }
          }
        }
      } else {
        for (int j = 1; j <= n; ++j) {
          if (!x[jx - 1].isZero() || !y[jy - 1].isZero()) {
            RS temp1 = alpha.multiply(y[jy - 1]);
            RS temp2 = alpha.multiply(x[jx - 1]);
            int ix = jx;
            int iy = jy;
            for (int i = j; i <= n; ++i) {
              int p = j * a_dim1 + i + pointer_a;
              a[p] = a[p].add(x[ix - 1].multiply(temp1)).add(y[iy - 1].multiply(temp2));
              ix += incx;
              iy += incy;
            }
          }
          jx += incx;
          jy += incy;
        }
      }
    }
    return 0;
  }
}