Dlasr.java
package org.mklab.sdpj.gpack.lapackwrap;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.gpack.blaswrap.BLAS;
/**
* @author koga
* @version $Revision$, 2009/04/25
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
*/
public class Dlasr<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>> {
/* -- LAPACK auxiliary routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
October 31, 1992
Purpose
=======
DLASR performs the transformation
A := P*A, when SIDE = 'L' or 'l' ( Left-hand side )
A := A*P', when SIDE = 'R' or 'r' ( Right-hand side )
where A is an m by n real matrix and P is an orthogonal matrix,
consisting of a sequence of plane rotations determined by the
parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'
and z = n when SIDE = 'R' or 'r' ):
When DIRECT = 'F' or 'f' ( Forward sequence ) then
P = P( z - 1 )*...*P( 2 )*P( 1 ),
and when DIRECT = 'B' or 'b' ( Backward sequence ) then
P = P( 1 )*P( 2 )*...*P( z - 1 ),
where P( k ) is a plane rotation matrix for the following planes:
when PIVOT = 'V' or 'v' ( Variable pivot ),
the plane ( k, k + 1 )
when PIVOT = 'T' or 't' ( Top pivot ),
the plane ( 1, k + 1 )
when PIVOT = 'B' or 'b' ( Bottom pivot ),
the plane ( k, z )
c( k ) and s( k ) must contain the cosine and sine that define the
matrix P( k ). The two by two plane rotation part of the matrix
P( k ), R( k ), is assumed to be of the form
R( k ) = ( c( k ) s( k ) ).
( -s( k ) c( k ) )
This version vectorises across rows of the array A when SIDE = 'L'.
Arguments
=========
SIDE (input) CHARACTER*1
Specifies whether the plane rotation matrix P is applied to
A on the left or the right.
= 'L': Left, compute A := P*A
= 'R': Right, compute A:= A*P'
DIRECT (input) CHARACTER*1
Specifies whether P is a forward or backward sequence of
plane rotations.
= 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )
= 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )
PIVOT (input) CHARACTER*1
Specifies the plane for which P(k) is a plane rotation
matrix.
= 'V': Variable pivot, the plane (k,k+1)
= 'T': Top pivot, the plane (1,k+1)
= 'B': Bottom pivot, the plane (k,z)
M (input) INTEGER
The number of rows of the matrix A. If m <= 1, an immediate
return is effected.
N (input) INTEGER
The number of columns of the matrix A. If n <= 1, an
immediate return is effected.
C, S (input) DOUBLE PRECISION arrays, dimension
(M-1) if SIDE = 'L'
(N-1) if SIDE = 'R'
c(k) and s(k) contain the cosine and sine that define the
matrix P(k). The two by two plane rotation part of the
matrix P(k), R(k), is assumed to be of the form
R( k ) = ( c( k ) s( k ) ).
( -s( k ) c( k ) )
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
The m by n matrix A. On exit, A is overwritten by P*A if
SIDE = 'R' or by A*P' if SIDE = 'L'.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
*/
/**
* @param side side
* @param pivot pivot
* @param direct direct
* @param m m
* @param n n
* @param c input
* @param s input
* @param a input/output
* @param lda lda
* @return result
*/
int execute(String side, String pivot, String direct, int m, int n, RS[] c, RS[] s, RS[] a, int lda) {
int pointer_c = -1;
int pointer_s = -1;
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int info = 0;
if (!(BLAS.lsame(side, "L") || BLAS.lsame(side, "R"))) { //$NON-NLS-1$ //$NON-NLS-2$
info = 1;
} else if (!(BLAS.lsame(pivot, "V") || BLAS.lsame(pivot, "T") || BLAS.lsame(pivot, "B"))) { //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
info = 2;
} else if (!(BLAS.lsame(direct, "F") || BLAS.lsame(direct, "B"))) { //$NON-NLS-1$ //$NON-NLS-2$
info = 3;
} else if (m < 0) {
info = 4;
} else if (n < 0) {
info = 5;
} else if (lda < Math.max(1, m)) {
info = 9;
}
if (info != 0) {
BLAS.xerbla("DLASR ", info); //$NON-NLS-1$
return info;
}
/* Quick return if possible */
if (m == 0 || n == 0) {
return info;
}
if (BLAS.lsame(side, "L")) { //$NON-NLS-1$
/* Form P * A */
if (BLAS.lsame(pivot, "V")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 1; j <= m - 1; ++j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + j + pointer_a;
int p2 = i * a_dim1 + j + 1 + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = m - 1; j >= 1; --j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + j + pointer_a;
int p2 = i * a_dim1 + j + 1 + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
}
} else if (BLAS.lsame(pivot, "T")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 2; j <= m; ++j) {
RS ctemp = c[pointer_c + j - 1];
RS stemp = s[pointer_s + j - 1];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + 1 + pointer_a;
int p2 = i * a_dim1 + j + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = m; j >= 2; --j) {
RS ctemp = c[pointer_c + j - 1];
RS stemp = s[pointer_s + j - 1];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + 1 + pointer_a;
int p2 = i * a_dim1 + j + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
}
} else if (BLAS.lsame(pivot, "B")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 1; j <= m - 1; ++j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + m + pointer_a;
int p2 = i * a_dim1 + j + pointer_a;
RS temp = a[p2];
a[p2] = (stemp.multiply(a[p1])).add(ctemp.multiply(temp));
a[p1] = (ctemp.multiply(a[p1])).subtract(stemp.multiply(temp));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = m - 1; j >= 1; --j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= n; ++i) {
int p1 = i * a_dim1 + m + pointer_a;
int p2 = i * a_dim1 + j + pointer_a;
RS temp = a[p2];
a[p2] = (stemp.multiply(a[p1])).add(ctemp.multiply(temp));
a[p1] = (ctemp.multiply(a[p1])).subtract(stemp.multiply(temp));
}
}
}
}
}
} else if (BLAS.lsame(side, "R")) { //$NON-NLS-1$
/* Form A * P' */
if (BLAS.lsame(pivot, "V")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 1; j <= n - 1; ++j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = (j) * a_dim1 + i + pointer_a;
int p2 = (j + 1) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = n - 1; j >= 1; --j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = (j) * a_dim1 + i + pointer_a;
int p2 = (j + 1) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
}
} else if (BLAS.lsame(pivot, "T")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 2; j <= n; ++j) {
RS ctemp = c[pointer_c + j - 1];
RS stemp = s[pointer_s + j - 1];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = a_dim1 + i + pointer_a;
int p2 = (j) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = n; j >= 2; --j) {
RS ctemp = c[pointer_c + j - 1];
RS stemp = s[pointer_s + j - 1];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = a_dim1 + i + pointer_a;
int p2 = (j) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (ctemp.multiply(temp)).subtract(stemp.multiply(a[p1]));
a[p1] = (stemp.multiply(temp)).add(ctemp.multiply(a[p1]));
}
}
}
}
} else if (BLAS.lsame(pivot, "B")) { //$NON-NLS-1$
if (BLAS.lsame(direct, "F")) { //$NON-NLS-1$
for (int j = 1; j <= n - 1; ++j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = n * a_dim1 + i + pointer_a;
int p2 = (j) * a_dim1 + i + pointer_a;
int p3 = (n) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (stemp.multiply(a[p3])).add(ctemp.multiply(temp));
a[p1] = (ctemp.multiply(a[p3])).subtract(stemp.multiply(temp));
}
}
}
} else if (BLAS.lsame(direct, "B")) { //$NON-NLS-1$
for (int j = n - 1; j >= 1; --j) {
RS ctemp = c[pointer_c + j];
RS stemp = s[pointer_s + j];
if (!ctemp.isUnit() || !stemp.isZero()) {
for (int i = 1; i <= m; ++i) {
int p1 = (n) * a_dim1 + i + pointer_a;
int p2 = (j) * a_dim1 + i + pointer_a;
RS temp = a[p2];
a[p2] = (stemp.multiply(a[p1])).add(ctemp.multiply(temp));
a[p1] = (ctemp.multiply(a[p1])).subtract(stemp.multiply(temp));
}
}
}
}
}
}
return info;
}
}