Dtrmv.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 Dtrmv<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
=======
DTRMV performs one of the matrix-vector operations
x := A*x, or x := A'*x,
where x is an n element vector and A is an n by n unit, or non-unit,
upper or lower triangular matrix.
Parameters
==========
UPLO - CHARACTER*1.
On entry, UPLO specifies whether the matrix is an upper or
lower triangular matrix as follows:
UPLO = 'U' or 'u' A is an upper triangular matrix.
UPLO = 'L' or 'l' A is a lower triangular matrix.
Unchanged on exit.
TRANS - CHARACTER*1.
On entry, TRANS specifies the operation to be performed as
follows:
TRANS = 'N' or 'n' x := A*x.
TRANS = 'T' or 't' x := A'*x.
TRANS = 'C' or 'c' x := A'*x.
Unchanged on exit.
DIAG - CHARACTER*1.
On entry, DIAG specifies whether or not A is unit
triangular as follows:
DIAG = 'U' or 'u' A is assumed to be unit triangular.
DIAG = 'N' or 'n' A is not assumed to be unit
triangular.
Unchanged on exit.
N - INTEGER.
On entry, N specifies the order of the matrix A.
N must be at least 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 matrix and the strictly lower triangular part of
A is not referenced.
Before entry with UPLO = 'L' or 'l', the leading n by n
lower triangular part of the array A must contain the lower
triangular matrix and the strictly upper triangular part of
A is not referenced.
Note that when DIAG = 'U' or 'u', the diagonal elements of
A are not referenced either, but are assumed to be unity.
Unchanged on exit.
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.
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. On exit, X is overwritten with the
tranformed vector x.
INCX - INTEGER.
On entry, INCX specifies the increment for the elements of
X. INCX must not be zero.
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 trans trans
* @param diag diag
* @param n n
* @param a a
* @param lda lda
* @param x x
* @param incx incx
* @return result
*/
public int execute(String uplo, String trans, String diag, int n, RS[] a, int lda, RS[] x, int incx) {
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 (!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 (!BLAS.lsame(diag, "U") && !BLAS.lsame(diag, "N")) { //$NON-NLS-1$ //$NON-NLS-2$
info = 3;
} else if (n < 0) {
info = 4;
} else if (lda < Math.max(1, n)) {
info = 6;
} else if (incx == 0) {
info = 8;
}
if (info != 0) {
BLAS.xerbla("DTRMV ", info); //$NON-NLS-1$
return 0;
}
// Quick return if possible.
if (n == 0) {
return 0;
}
boolean nounit = BLAS.lsame(diag, "N"); //$NON-NLS-1$
/*
* Set up the start point in X if the increment is not unity.
* This will be (N - 1 )INCX too small for descending loops.
*/
int kx = 0;
if (incx <= 0) {
kx = 1 - (n - 1) * incx;
} else if (incx != 1) {
kx = 1;
}
/*
* Start the operations. In this version the elements of A are accessed
* sequentially with one pass through A.
*/
if (BLAS.lsame(trans, "N")) { //$NON-NLS-1$
// Form x := Ax.
if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
if (incx == 1) {
for (int j = 1; j <= n; ++j) {
if (!x[j - 1].isZero()) {
RS temp = x[j - 1];
for (int i = 1; i <= j - 1; ++i) {
x[i - 1] = x[i - 1].add(temp.multiply(a[j * a_dim1 + i + pointer_a]));
}
if (nounit) {
x[j - 1] = temp.multiply(a[j * a_dim1 + j + pointer_a]);
}
}
}
} else {
int jx = kx;
for (int j = 1; j <= n; ++j) {
if (!x[jx - 1].isZero()) {
RS temp = x[jx - 1];
int ix = kx;
for (int i = 1; i <= j - 1; ++i) {
x[ix - 1] = x[ix - 1].add(temp.multiply(a[j * a_dim1 + i + pointer_a]));
ix += incx;
}
if (nounit) {
x[jx - 1] = x[jx - 1].multiply(a[j * a_dim1 + j + pointer_a]);
}
}
jx += incx;
}
}
} else {
if (incx == 1) {
for (int j = n; j >= 1; --j) {
if (!x[j - 1].isZero()) {
RS temp = x[j - 1];
for (int i = n; i >= j + 1; --i) {
x[i - 1] = x[i - 1].add(temp.multiply(a[j * a_dim1 + i + pointer_a]));
}
if (nounit) {
x[j - 1] = x[j - 1].multiply(a[j * a_dim1 + j + pointer_a]);
}
}
}
} else {
kx += (n - 1) * incx;
int jx = kx;
for (int j = n; j >= 1; --j) {
if (!x[jx - 1].isZero()) {
RS temp = x[jx - 1];
int ix = kx;
for (int i = n; i >= j + 1; --i) {
x[ix - 1] = x[ix - 1].add(temp.multiply(a[j * a_dim1 + i + pointer_a]));
ix -= incx;
}
if (nounit) {
x[jx - 1] = x[jx - 1].multiply(a[j * a_dim1 + j + pointer_a]);
}
}
jx -= incx;
}
}
}
} else {
/* Form x := A'x. */
if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
if (incx == 1) {
for (int j = n; j >= 1; --j) {
RS temp = x[j - 1];
if (nounit) {
temp = temp.multiply(a[j * a_dim1 + j + pointer_a]);
}
for (int i = j - 1; i >= 1; --i) {
temp = temp.add(a[j * a_dim1 + i + pointer_a].multiply(x[i - 1]));
}
x[j - 1] = temp;
}
} else {
int jx = kx + (n - 1) * incx;
for (int j = n; j >= 1; --j) {
RS temp = x[jx - 1];
int ix = jx;
if (nounit) {
temp = temp.multiply(a[j * a_dim1 + j + pointer_a]);
}
for (int i = j - 1; i >= 1; --i) {
ix -= incx;
temp = temp.add(a[j * a_dim1 + i + pointer_a].multiply(x[ix - 1]));
}
x[jx - 1] = temp;
jx -= incx;
}
}
} else {
if (incx == 1) {
for (int j = 1; j <= n; ++j) {
RS temp = x[j - 1];
if (nounit) {
temp = temp.multiply(a[j * a_dim1 + j + pointer_a]);
}
for (int i = j + 1; i <= n; ++i) {
temp = temp.add(a[j * a_dim1 + i + pointer_a].multiply(x[i - 1]));
}
x[j - 1] = temp;
}
} else {
int jx = kx;
for (int j = 1; j <= n; ++j) {
RS temp = x[jx - 1];
int ix = jx;
if (nounit) {
temp = temp.multiply(a[j * a_dim1 + j + pointer_a]);
}
for (int i = j + 1; i <= n; ++i) {
ix += incx;
temp = temp.add(a[j * a_dim1 + i + pointer_a].multiply(x[ix - 1]));
}
x[jx - 1] = temp;
jx += incx;
}
}
}
}
return 0;
}
}