Dsyev.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;
import org.mklab.sdpj.tool.Tools;
/**
* @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 Dsyev<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 driver 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
=======
DSYEV computes all eigenvalues and, optionally, eigenvectors of a
real symmetric matrix A.
Arguments
=========
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
On entry, 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. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
W (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK (input) INTEGER
The length of the array WORK. LWORK >= max(1,3*N-1).
For optimal efficiency, LWORK >= (NB+2)*N,
where NB is the blocksize for DSYTRD returned by ILAENV.
If LWORK = -1, then a workspace query is assumed; the routine
only calculates the optimal size of the WORK array, returns
this value as the first entry of the WORK array, and no error
message related to LWORK is issued by XERBLA.
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, the algorithm failed to converge; i
off-diagonal elements of an intermediate tridiagonal
form did not converge to zero.
=====================================================================
*/
//work[3]を監視
/**
* @param jobz jobz
* @param uplo uplo
* @param n n
* @param a a
* @param lda lda
* @param w w
* @param work work
* @param lwork lwork
* @return result
*/
int execute(String jobz, String uplo, int n, RS[] a, int lda, RS[] w, RS[] work, int lwork) {
final RS unit = a[0].createUnit();
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int pointer_w = -1;
int pointer_work = -1;
boolean wantz = BLAS.lsame(jobz, "V"); //$NON-NLS-1$
boolean lower = BLAS.lsame(uplo, "L"); //$NON-NLS-1$
boolean lquery = lwork == -1;
int info = 0;
if (!(wantz || BLAS.lsame(jobz, "N"))) { //$NON-NLS-1$
info = -1;
} else if (!(lower || BLAS.lsame(uplo, "U"))) { //$NON-NLS-1$
info = -2;
} else if (n < 0) {
info = -3;
} else if (lda < Math.max(1, n)) {
info = -5;
} else /* if(complicated condition) */{
if (lwork < Math.max(1, n * 3 - 1) && !lquery) {
info = -8;
}
}
int lwkopt = 0;
if (info == 0) {
int nb = Clapack.ilaenv(1, "DSYTRD", uplo, n, -1, -1, -1, 6, 1, unit); //$NON-NLS-1$
lwkopt = Math.max(1, (nb + 2) * n);
work[1 + pointer_work] = unit.create(lwkopt);
}
if (info != 0) {
BLAS.xerbla("DSYEV ", -info); //$NON-NLS-1$
return info;
} else if (lquery) {
return info;
}
/* Quick return if possible */
if (n == 0) {
work[1 + pointer_work] = unit.create(1);
return info;
}
if (n == 1) {
w[1 + pointer_w] = a[a_dim1 + 1 + pointer_a];
work[1 + pointer_work] = unit.create(3);
if (wantz) {
a[a_dim1 + 1 + pointer_a] = unit.create(1);// 1.;
}
return info;
}
/* Get machine constants. */
//RS safmin = Clapack.dlamch("Safe minimum",unit); //$NON-NLS-1$
RS safmin = Clapack.dlamch("S",unit); //$NON-NLS-1$
//RS eps = Clapack.dlamch("Precision",unit); //$NON-NLS-1$
RS eps = Clapack.dlamch("E",unit); //$NON-NLS-1$
RS smlnum = safmin.divide(eps);
RS bignum = unit.create(1).divide(smlnum);
RS rmin = smlnum.sqrt();
RS rmax = bignum.sqrt();
/* Scale matrix to allowable range, if necessary. */
RS[] aTemp = unit.createArray(a.length - (a_offset + pointer_a));
RS[] workTemp = unit.createArray(work.length - (1 + pointer_work));
System.arraycopy(a, a_offset + pointer_a, aTemp, 0, aTemp.length);
System.arraycopy(work, 1 + pointer_work, workTemp, 0, workTemp.length);
RS anrm = Clapack.dlansy("M", uplo, n, aTemp, lda, workTemp); //$NON-NLS-1$
System.arraycopy(aTemp, 0, a, a_offset + pointer_a, aTemp.length);
System.arraycopy(workTemp, 0, work, 1 + pointer_work, workTemp.length);
int iscale = 0;
RS sigma = null;
if (anrm.isGreaterThan(0) && anrm.isLessThan(rmin)) {
iscale = 1;
sigma = rmin.divide(anrm);
} else if (anrm.isGreaterThan(rmax)) {
iscale = 1;
sigma = rmax.divide(anrm);
}
if (iscale == 1) {
RS[] a_temp = unit.createArray(a.length - (a_offset + pointer_a));
System.arraycopy(a, a_offset + pointer_a, a_temp, 0, a_temp.length);
RS c_b17 = unit.create(1);
info = Clapack.dlascl(uplo, 0, 0, c_b17, sigma, n, n, a_temp, lda);
System.arraycopy(a_temp, 0, a, a_offset + pointer_a, a_temp.length);
}
/* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */
int inde = 1;
int indtau = inde + n;
int indwrk = indtau + n;
int llwork = lwork - indwrk + 1;
RS[] wTemp = unit.createArray(w.length - (1 + pointer_w));
RS[] workTemp1 = unit.createArray(work.length - (inde + pointer_work));
RS[] workTemp2 = unit.createArray(work.length - (indtau + pointer_work));
RS[] workTemp3 = unit.createArray(work.length - (indwrk + pointer_work));
aTemp = unit.createArray(a.length - (a_offset + pointer_a));
System.arraycopy(a, a_offset + pointer_a, aTemp, 0, aTemp.length);
System.arraycopy(w, 1 + pointer_w, wTemp, 0, wTemp.length);
System.arraycopy(work, inde + pointer_work, workTemp1, 0, workTemp1.length);
System.arraycopy(work, indtau + pointer_work, workTemp2, 0, workTemp2.length);
System.arraycopy(work, indwrk + pointer_work, workTemp3, 0, workTemp3.length);
info = Clapack.dsytrd(uplo, n, aTemp, lda, wTemp, workTemp1, workTemp2, workTemp3, llwork);
System.arraycopy(aTemp, 0, a, a_offset + pointer_a, aTemp.length);
System.arraycopy(wTemp, 0, w, 1 + pointer_w, wTemp.length);
System.arraycopy(workTemp1, 0, work, inde + pointer_work, workTemp1.length);
System.arraycopy(workTemp2, 0, work, indtau + pointer_work, workTemp2.length);
System.arraycopy(workTemp3, 0, work, indwrk + pointer_work, workTemp3.length);
int lopt = Tools.NumericalScalarToInt(work[indwrk + pointer_work].add(n << 1));
/*
* For eigenvalues only, call DSTERF. For eigen vectors, first call DORGTR to
* generate the orthogonal matrix, then call DSTEQR.
*/
if (!wantz) {
wTemp = unit.createArray(w.length - (1 + pointer_w));
workTemp = unit.createArray(work.length - (inde + pointer_work));
System.arraycopy(w, 1 + pointer_w, wTemp, 0, wTemp.length);
System.arraycopy(work, inde + pointer_work, workTemp, 0, workTemp.length);
info = Clapack.dsterf(n, wTemp, workTemp);
System.arraycopy(wTemp, 0, w, 1 + pointer_w, wTemp.length);
System.arraycopy(workTemp, 0, work, inde + pointer_work, workTemp.length);
} else {
RS[] a_temp = unit.createArray(a.length - (a_offset + pointer_a));
RS[] work_indtau = unit.createArray(work.length - (indtau + pointer_work));// [a.length-indtau];
RS[] work_indwrk = unit.createArray(work.length - (indwrk + pointer_work));// [a.length- indwrk ];
RS[] work_inde = unit.createArray(work.length - (inde + pointer_work));// [a.length - inde ];
RS[] w1 = unit.createArray(w.length - (1 + pointer_w));// [a.length - 1 ];
System.arraycopy(a, a_offset + pointer_a, a_temp, 0, a_temp.length);
System.arraycopy(work, indtau + pointer_work, work_indtau, 0, work_indtau.length);
System.arraycopy(work, indwrk + pointer_work, work_indwrk, 0, work_indwrk.length);
System.arraycopy(w, 1 + pointer_w, w1, 0, w1.length);
System.arraycopy(work, inde + pointer_work, work_inde, 0, work_inde.length);
info = Clapack.dorgtr(uplo, n, a_temp, lda, work_indtau, work_indwrk, llwork);
//debug
info = Clapack.dsteqr(jobz, n, w1, work_inde, a_temp, lda, work_indtau);
System.arraycopy(a_temp, 0, a, a_offset + pointer_a, a_temp.length);
System.arraycopy(work_inde, 0, work, inde + pointer_work, work_inde.length);
System.arraycopy(work_indtau, 0, work, indtau + pointer_work, work_indtau.length);
System.arraycopy(work_indwrk, 0, work, indwrk + pointer_work, work_indwrk.length);
System.arraycopy(w1, 0, w, 1 + pointer_w, w1.length);
}
/* If matrix was scaled, then rescale eigenvalues appropriately. */
if (iscale == 1) {
int imax;
if (info == 0) {
imax = n;
} else {
imax = info - 1;
}
RS d1 = unit.create(1).divide(sigma);
RS[] w1 = unit.createArray(w.length - (1 + pointer_w));
System.arraycopy(w, 1 + pointer_w, w1, 0, w1.length);
BLAS.dscal(imax, d1, w1, 1);
System.arraycopy(w1, 0, w, 1 + pointer_w, w1.length);
}
/* Set WORK(1) to optimal workspace size. */
work[1 + pointer_work] = unit.create(lwkopt);
return info;
}
}