Dlascl.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 takafumi
* @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 Dlascl<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
February 29, 1992
Purpose
=======
DLASCL multiplies the M by N real matrix A by the real scalar
CTO/CFROM. This is done without over/underflow as long as the final
result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
A may be full, upper triangular, lower triangular, upper Hessenberg,
or banded.
Arguments
=========
TYPE (input) CHARACTER*1
TYPE indices the storage type of the input matrix.
= 'G': A is a full matrix.
= 'L': A is a lower triangular matrix.
= 'U': A is an upper triangular matrix.
= 'H': A is an upper Hessenberg matrix.
= 'B': A is a symmetric band matrix with lower bandwidth KL
and upper bandwidth KU and with the only the lower
half stored.
= 'Q': A is a symmetric band matrix with lower bandwidth KL
and upper bandwidth KU and with the only the upper
half stored.
= 'Z': A is a band matrix with lower bandwidth KL and upper
bandwidth KU.
KL (input) INTEGER
The lower bandwidth of A. Referenced only if TYPE = 'B',
'Q' or 'Z'.
KU (input) INTEGER
The upper bandwidth of A. Referenced only if TYPE = 'B',
'Q' or 'Z'.
CFROM (input) DOUBLE PRECISION
CTO (input) DOUBLE PRECISION
The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
without over/underflow if the final result CTO*A(I,J)/CFROM
can be represented without over/underflow. CFROM must be
nonzero.
M (input) INTEGER
The number of rows of the matrix A. M >= 0.
N (input) INTEGER
The number of columns of the matrix A. N >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,M)
The matrix to be multiplied by CTO/CFROM. See TYPE for the
storage type.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
INFO (output) INTEGER
0 - successful exit
<0 - if INFO = -i, the i-th argument had an illegal value.
===================================================================== */
/**
* @param type typye
* @param kl kl
* @param ku ku
* @param cfrom cfrom
* @param cto cto
* @param m m
* @param n n
* @param a a
* @param lda lda
* @return result
*/
int execute(String type, int kl, int ku, RS cfrom, RS cto, int m, int n, RS[] a, int lda) {
final RS unit = a[0].createUnit();
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int info = 0;
int itype;
if (BLAS.lsame(type, "G")) { //$NON-NLS-1$
itype = 0;
} else if (BLAS.lsame(type, "L")) { //$NON-NLS-1$
itype = 1;
} else if (BLAS.lsame(type, "U")) { //$NON-NLS-1$
itype = 2;
} else if (BLAS.lsame(type, "H")) { //$NON-NLS-1$
itype = 3;
} else if (BLAS.lsame(type, "B")) { //$NON-NLS-1$
itype = 4;
} else if (BLAS.lsame(type, "Q")) { //$NON-NLS-1$
itype = 5;
} else if (BLAS.lsame(type, "Z")) { //$NON-NLS-1$
itype = 6;
} else {
itype = -1;
}
if (itype == -1) {
info = -1;
} else if (cfrom.isZero()) {
info = -4;
} else if (m < 0) {
info = -6;
} else if (n < 0 || itype == 4 && n != m || itype == 5 && n != m) {
info = -7;
} else if (itype <= 3 && lda < Math.max(1, m)) {
info = -9;
} else if (itype >= 4) {
/* Computing MAX */
if (kl < 0 || kl > Math.max(m - 1, 0)) {
info = -2;
} else /* if(complicated condition) */{
if (ku < 0 || ku > Math.max(n - 1, 0) || (itype == 4 || itype == 5) && kl != ku) {
info = -3;
} else if (itype == 4 && lda < kl + 1 || itype == 5 && lda < ku + 1 || itype == 6 && lda < (kl << 1) + ku + 1) {
info = -9;
}
}
}
if (info != 0) {
BLAS.xerbla("DLASCL", -(info)); //$NON-NLS-1$
return info;
}
/* Quick return if possible */
if (n == 0 || m == 0) {
return info;
}
/* Get machine parameters */
RS smlnum = Clapack.dlamch("S", unit); //$NON-NLS-1$
RS bignum = unit.create(1).divide(smlnum);
RS cfromc = cfrom;
RS ctoc = cto;
boolean done;
do {
RS cfrom1 = cfromc.multiply(smlnum);
RS cto1 = ctoc.divide(bignum);
RS mul;
if (cfrom1.abs().isGreaterThan(ctoc.abs()) && !ctoc.isZero()) {
mul = smlnum;
done = false;
cfromc = cfrom1;
} else if (cto1.abs().isGreaterThan(cfromc.abs())) {
mul = bignum;
done = false;
ctoc = cto1;
} else {
mul = ctoc.divide(cfromc);
done = true;
}
if (itype == 0) {
// Full matrix
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= m; ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 1) {
// Lower triangular matrix
for (int j = 1; j <= n; ++j) {
for (int i = j; i <= m; ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 2) {
// Upper triangular matrix
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= Math.min(j, m); ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 3) {
// Upper Hessenberg matrix
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= Math.min(j + 1, m); ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 4) {
// Lower half of a symmetric band matrix
for (int j = 1; j <= n; ++j) {
int iMax = Math.min(kl + 1, (n + 1) - j);
for (int i = 1; i <= iMax; ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 5) {
// Upper half of a symmetric band matrix
for (int j = 1; j <= n; ++j) {
int iMin = Math.max((ku + 2) - j, 1);
for (int i = iMin; i <= ku + 1; ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
} else if (itype == 6) {
// Band matrix
for (int j = 1; j <= n; ++j) {
int iMin = Math.max((kl + ku + 2) - j, kl + 1);
int iMax = Math.min((kl << 1) + ku + 1, (kl + ku + 1 + m) - j);
for (int i = iMin; i <= iMax; ++i) {
int p = j * a_dim1 + i + pointer_a;
a[p] = a[p].multiply(mul);
}
}
}
} while (!done);
return info;
}
}