Dsyrk.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 Dsyrk<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
=======
DSYRK performs one of the symmetric rank k operations
C := alpha*A*A' + beta*C,
or
C := alpha*A'*A + beta*C,
where alpha and beta are scalars, C is an n by n symmetric matrix
and A is an n by k matrix in the first case and a k by n matrix
in the second case.
Parameters
==========
UPLO - CHARACTER*1.
On entry, UPLO specifies whether the upper or lower
triangular part of the array C is to be referenced as
follows:
UPLO = 'U' or 'u' Only the upper triangular part of C
is to be referenced.
UPLO = 'L' or 'l' Only the lower triangular part of C
is to be referenced.
Unchanged on exit.
TRANS - CHARACTER*1.
On entry, TRANS specifies the operation to be performed as
follows:
TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
Unchanged on exit.
N - INTEGER.
On entry, N specifies the order of the matrix C. N must be
at least zero.
Unchanged on exit.
K - INTEGER.
On entry with TRANS = 'N' or 'n', K specifies the number
of columns of the matrix A, and on entry with
TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
of rows of the matrix A. K must be at least zero.
Unchanged on exit.
ALPHA - DOUBLE PRECISION.
On entry, ALPHA specifies the scalar alpha.
Unchanged on exit.
A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
k when TRANS = 'N' or 'n', and is n otherwise.
Before entry with TRANS = 'N' or 'n', the leading n by k
part of the array A must contain the matrix A, otherwise
the leading k by n part of the array A must contain the
matrix A.
Unchanged on exit.
LDA - INTEGER.
On entry, LDA specifies the first dimension of A as declared
in the calling (sub) program. When TRANS = 'N' or 'n'
then LDA must be at least max( 1, n ), otherwise LDA must
be at least max( 1, k ).
Unchanged on exit.
BETA - DOUBLE PRECISION.
On entry, BETA specifies the scalar beta.
Unchanged on exit.
C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
Before entry with UPLO = 'U' or 'u', the leading n by n
upper triangular part of the array C must contain the upper
triangular part of the symmetric matrix and the strictly
lower triangular part of C is not referenced. On exit, the
upper triangular part of the array C 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 C must contain the lower
triangular part of the symmetric matrix and the strictly
upper triangular part of C is not referenced. On exit, the
lower triangular part of the array C is overwritten by the
lower triangular part of the updated matrix.
LDC - INTEGER.
On entry, LDC specifies the first dimension of C as declared
in the calling (sub) program. LDC must be at least
max( 1, n ).
Unchanged on exit.
Level 3 Blas routine.
-- Written on 8-February-1989.
Jack Dongarra, Argonne National Laboratory.
Iain Duff, AERE Harwell.
Jeremy Du Croz, Numerical Algorithms Group Ltd.
Sven Hammarling, Numerical Algorithms Group Ltd.
Test the input parameters.
Parameter adjustments */
/**
* @param uplo uplo
* @param trans trans
* @param n n
* @param k k
* @param alpha alpha
* @param a a
* @param lda lda
* @param beta beta
* @param c c
* @param ldc ldc
* @return result
*/
public int execute(String uplo, String trans, int n, int k, RS alpha, RS[] a, int lda, RS beta, RS[] c, int ldc) {
final RS unit = a[0].createUnit();
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int c_dim1 = ldc;
int c_offset = 1 + c_dim1 * 1;
int pointer_c = -c_offset;
int nrowa;
if (BLAS.lsame(trans, "N")) { //$NON-NLS-1$
nrowa = n;
} else {
nrowa = k;
}
boolean upper = BLAS.lsame(uplo, "U"); //$NON-NLS-1$
int info = 0;
if (!upper && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$
info = 1;
} else if (!BLAS.lsame(trans, "N") && !BLAS.lsame(trans, "T") && !BLAS.lsame(trans, "C")) { //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
info = 2;
} else if (n < 0) {
info = 3;
} else if (k < 0) {
info = 4;
} else if (lda < Math.max(1, nrowa)) {
info = 7;
} else if (ldc < Math.max(1, n)) {
info = 10;
}
if (info != 0) {
BLAS.xerbla("DSYRK ", info); //$NON-NLS-1$
return 0;
}
// Quick return if possible.
if (n == 0 || (alpha.isZero() || k == 0) && beta.equals(unit.createUnit())) {
return 0;
}
// And when alpha.eq.zero.
if (alpha.isZero()) {
if (upper) {
if (beta.isZero()) {
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= j; ++i) {
c[j * c_dim1 + i + pointer_c] = unit.createZero();
}
}
} else {
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= j; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = beta.multiply(c[p]);
}
}
}
} else {
if (beta.isZero()) {
for (int j = 1; j <= n; ++j) {
for (int i = j; i <= n; ++i) {
c[j * c_dim1 + i + pointer_c] = unit.createZero();
}
}
} else {
for (int j = 1; j <= n; ++j) {
for (int i = j; i <= n; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = beta.multiply(c[p]);
}
}
}
}
return 0;
}
if (BLAS.lsame(trans, "N")) { //$NON-NLS-1$
// Form C := alphaAA' + betaC.
if (upper) {
for (int j = 1; j <= n; ++j) {
if (beta.isZero()) {
for (int i = 1; i <= j; ++i) {
c[j * c_dim1 + i + pointer_c] = unit.createZero();
}
} else if (!beta.equals(unit.create(1))) {
for (int i = 1; i <= j; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = beta.multiply(c[p]);
}
}
for (int l = 1; l <= k; ++l) {
if (!a[l * a_dim1 + j + pointer_a].isZero()) {
RS temp = alpha.multiply(a[l * a_dim1 + j + pointer_a]);
for (int i = 1; i <= j; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = c[p].add(temp.multiply(a[l * a_dim1 + i + pointer_a]));
}
}
}
}
} else {
for (int j = 1; j <= n; ++j) {
if (beta.isZero()) {
for (int i = j; i <= n; ++i) {
c[j * c_dim1 + i + pointer_c] = unit.createZero();
}
} else if (!beta.equals(unit.createUnit())) {
for (int i = j; i <= n; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = beta.multiply(c[p]);
}
}
for (int l = 1; l <= k; ++l) {
if (!a[l * a_dim1 + j + pointer_a].isZero()) {
RS temp = alpha.multiply(a[l * a_dim1 + j + pointer_a]);
for (int i = j; i <= n; ++i) {
int p = j * c_dim1 + i + pointer_c;
c[p] = c[p].add(temp.multiply(a[l * a_dim1 + i + pointer_a]));
}
}
}
}
}
} else {
// Form C := alphaA'A + betaC.
if (upper) {
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= j; ++i) {
RS temp = unit.createZero();
for (int l = 1; l <= k; ++l) {
int p = i * a_dim1 + l + pointer_a;
temp = temp.add(a[p].multiply(a[j * a_dim1 + l + pointer_a]));
}
int p = j * c_dim1 + i + pointer_c;
if (beta.isZero()) {
c[p] = alpha.multiply(temp);
} else {
c[p] = (alpha.multiply(temp)).add(beta.multiply(c[p]));
}
}
}
} else {
for (int j = 1; j <= n; ++j) {
for (int i = j; i <= n; ++i) {
RS temp = unit.createZero();
for (int l = 1; l <= k; ++l) {
int p = i * a_dim1 + l + pointer_a;
temp = temp.add(a[p].multiply(a[j * a_dim1 + l + pointer_a]));
}
int p = j * c_dim1 + i + pointer_c;
if (beta.isZero()) {
c[p] = alpha.multiply(temp);
} else {
c[p] = (alpha.multiply(temp)).add(beta.multiply(c[p]));
}
}
}
}
}
return 0;
}
}