Dlansy.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 Dlansy<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>> {
/** */
private RS[] a;
/** */
private RS scale;
/** */
private RS sum;
/* -- 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
=======
DLANSY returns the value of the one norm, or the Frobenius norm, or
the infinity norm, or the element of largest absolute value of a
real symmetric matrix A.
Description
===========
DLANSY returns the value
DLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm'
(
( norm1(A), NORM = '1', 'O' or 'o'
(
( normI(A), NORM = 'I' or 'i'
(
( normF(A), NORM = 'F', 'f', 'E' or 'e'
where norm1 denotes the one norm of a matrix (maximum column sum),
normI denotes the infinity norm of a matrix (maximum row sum) and
normF denotes the Frobenius norm of a matrix (square root of sum of
squares). Note that max(abs(A(i,j))) is not a matrix norm.
Arguments
=========
NORM (input) CHARACTER*1
Specifies the value to be returned in DLANSY as described
above.
UPLO (input) CHARACTER*1
Specifies whether the upper or lower triangular part of the
symmetric matrix A is to be referenced.
= 'U': Upper triangular part of A is referenced
= 'L': Lower triangular part of A is referenced
N (input) INTEGER
The order of the matrix A. N >= 0. When N = 0, DLANSY is
set to zero.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
The symmetric matrix A. If UPLO = 'U', the leading n by n
upper triangular part of A contains the upper triangular part
of the matrix A, and the strictly lower triangular part of A
is not referenced. If UPLO = 'L', the leading n by n lower
triangular part of A contains the lower triangular part of
the matrix A, and the strictly upper triangular part of A is
not referenced.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(N,1).
WORK (workspace) DOUBLE PRECISION array, dimension (LWORK),
where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
WORK is not referenced.
=====================================================================
*/
/**
* サブルーチン
*
* @param norm norm
* @param uplo uplo
* @param n n
* @param ain ain
* @param lda lda
* @param workin workin
* @return result
*/
RS execute(String norm, String uplo, int n, RS[] ain, int lda, RS[] workin) {
final RS unit = ain[0].create(1);
this.a = ain;
RS[] work = workin;
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int pointer_work = -1;
RS value = null;
if (n == 0) {
value = unit.createZero();
} else if (BLAS.lsame(norm, "M")) { //$NON-NLS-1$
/* Find max(abs(A(i,j))). */
value = unit.create(0);
if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
for (int j = 1; j <= n; ++j) {
for (int i = 1; i <= j; ++i) {
RS d1 = this.a[(j) * a_dim1 + i + pointer_a];
value = value.max(d1.abs());
}
}
} else {
for (int j = 1; j <= n; ++j) {
for (int i = j; i <= n; ++i) {
RS d1 = this.a[(j) * a_dim1 + i + pointer_a];
value = value.max(d1.abs());
}
}
}
} else if (BLAS.lsame(norm, "I") || BLAS.lsame(norm, "O") || norm.charAt(0) == '1') { //$NON-NLS-1$ //$NON-NLS-2$
/* Find normI(A) ( = norm1(A), since A is symmetric). */
value = unit.create(0);
if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
for (int j = 1; j <= n; ++j) {
this.sum = unit.create(0);
for (int i = 1; i <= j - 1; ++i) {
RS d1 = this.a[(j) * a_dim1 + i + pointer_a];
RS absa = d1.abs();
this.sum = this.sum.add(absa);
work[i + pointer_work] = work[i + pointer_work].add(absa);
}
RS d1 = this.a[(j) * a_dim1 + j + pointer_a];
work[j + pointer_work] = this.sum.add(d1.abs());
}
for (int i = 1; i <= n; ++i) {
RS d2 = work[i + pointer_work];
value = value.max(d2);
}
} else {
for (int i = 1; i <= n; ++i) {
work[i + pointer_work] = unit.create(0);
}
for (int j = 1; j <= n; ++j) {
RS d1 = this.a[(j) * a_dim1 + j + pointer_a];
this.sum = work[j + pointer_work].add(d1.abs());
for (int i = j + 1; i <= n; ++i) {
d1 = this.a[(j) * a_dim1 + i + pointer_a];
RS absa = d1.abs();
this.sum = this.sum.add(absa);
work[i + pointer_work] = work[i + pointer_work].add(absa);
}
value = value.max(this.sum);
}
}
} else if (BLAS.lsame(norm, "F") || BLAS.lsame(norm, "E")) { //$NON-NLS-1$ //$NON-NLS-2$
/* Find normF(A). */
this.scale = unit.create(0);
this.sum = unit.create(1);
if (BLAS.lsame(uplo, "U")) { //$NON-NLS-1$
for (int j = 2; j <= n; ++j) {
dlassq(j - 1, (j) * a_dim1 + 1 + pointer_a, 1);
}
} else {
for (int j = 1; j <= n - 1; ++j) {
dlassq(n - j, (j) * a_dim1 + j + 1 + pointer_a, 1);
}
}
this.sum = this.sum.multiply(2);
dlassq(n, a_offset + pointer_a, lda + 1);
value = this.scale.multiply(this.sum.sqrt());
}
RS ret_val = value;
return ret_val;
}
/* -- LAPACK auxiliary routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
June 30, 1999
Purpose
=======
DLASSQ returns the values scl and smsq such that
( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is
assumed to be non-negative and scl returns the value
scl = max( scale, abs( x( i ) ) ).
scale and sumsq must be supplied in SCALE and SUMSQ and
scl and smsq are overwritten on SCALE and SUMSQ respectively.
The routine makes only one pass through the vector x.
Arguments
=========
N (input) INTEGER
The number of elements to be used from the vector X.
X (input) DOUBLE PRECISION array, dimension (N)
The vector for which a scaled sum of squares is computed.
x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
INCX (input) INTEGER
The increment between successive values of the vector X.
INCX > 0.
SCALE (input/output) DOUBLE PRECISION
On entry, the value scale in the equation above.
On exit, SCALE is overwritten with scl , the scaling factor
for the sum of squares.
SUMSQ (input/output) DOUBLE PRECISION
On entry, the value sumsq in the equation above.
On exit, SUMSQ is overwritten with smsq , the basic sum of
squares from which scl has been factored out.
=====================================================================
*/
/**
* @param n n
* @param xIndex xIndex
* @param incx incx
* @return result
*/
private int dlassq(int n, int xIndex, int incx) {
RS[] x = this.a;
int pointer_x = xIndex;
--pointer_x;
if (n > 0) {
int i1 = (n - 1) * incx + 1;
int i2 = incx;
for (int ix = 1; i2 < 0 ? ix >= i1 : ix <= i1; ix += i2) {
if (!x[ix + pointer_x].isZero()) {
RS d = x[ix + pointer_x];
RS absxi = d.abs();
if (this.scale.isLessThan(absxi)) {
d = this.scale.divide(absxi);
this.sum = this.sum.multiply((d.multiply(d))).add(1);
this.scale = absxi;
} else {
d = absxi.divide(this.scale);
this.sum = this.sum.add(d.multiply(d));
}
}
/* L10: */
}
}
return 0;
}
}