Dsterf.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.gpack.f2clibs.LibF77;
/**
* 無平方根QL法またはQR法を使って、対称三重対角行列のすべての固有値を求める
*
* @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 Dsterf<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[] d;
/** */
private RS[] e;
/* -- LAPACK 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
=======
DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
using the Pal-Walker-Kahan variant of the QL or QR algorithm.
Arguments
=========
N (input) INTEGER
The order of the matrix. N >= 0.
D (input/output) DOUBLE PRECISION array, dimension (N)
On entry, the n diagonal elements of the tridiagonal matrix.
On exit, if INFO = 0, the eigenvalues in ascending order.
RS (input/output) DOUBLE PRECISION array, dimension (N-1)
On entry, the (n-1) subdiagonal elements of the tridiagonal
matrix.
On exit, RS has been destroyed.
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: the algorithm failed to find all of the eigenvalues in
a total of 30*N iterations; if INFO = i, then i
elements of RS have not converged to zero.
=====================================================================
*/
/**
* @param n n
* @param di di
* @param ei ei
* @return result
*/
public int execute(int n, RS[] di, RS[] ei) {
final RS unit = di[0].createUnit();
this.d = di;
this.e = ei;
//E rt1 = unit.create(0);
//E rt2 = unit.create(0);
int pointer_e = -1;
int pointer_d = -1;
int info = 0;
/* Quick return if possible */
if (n < 0) {
info = -1;
BLAS.xerbla("DSTERF", -info); //$NON-NLS-1$
return info;
}
if (n <= 1) {
return info;
}
/* Determine the unit roundoff for this environment. */
RS eps = Clapack.dlamch("E",unit); //$NON-NLS-1$
RS eps2 = eps.multiply(eps);
RS safmin = Clapack.dlamch("S",unit); //$NON-NLS-1$
RS safmax = unit.create(1).divide(safmin);
RS ssfmax = (safmax.sqrt()).divide(3);
RS ssfmin = (safmin.sqrt()).divide(eps2);
//System.out.println(eps + " \n" + eps2 + "\n" + safmin + "\n" + safmax + "\n" + ssfmax + "\n" + ssfmin + "\n"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$
/* Compute the eigenvalues of the tridiagonal matrix. */
int nmaxit = n * 30;
RS sigma = unit.create(0);
int jtot = 0;
/*
* Determine where the matrix splits and choose QL or QR iteration for each
* block, according to whether top or bottom diagonal element is smaller.
*/
int l1 = 1;
int l;
int lend;
int iscale;
int lendsv;
int lsv;
RS anorm;
// L10:
do {
do {
if (l1 > n) {
info = dlasrt("I", n, 1 + pointer_d); //$NON-NLS-1$
return info;
}
boolean isThrowGotoL30Statement = false;
if (l1 > 1) {
this.e[l1 - 1 + pointer_e] = unit.create(0);
}
int m;
for (m = l1; m <= n - 1; ++m) {
RS d1 = this.d[m + pointer_d];
RS d2 = this.d[m + 1 + pointer_d];
RS d3 = this.e[m + pointer_e];
if (d3.abs().isLessThanOrEquals((d1.abs()).sqrt().multiply(d2.abs().sqrt().multiply(eps)))) {
this.e[m + pointer_e] = unit.create(0);
// goto L30;
isThrowGotoL30Statement = true;
break;
}
/* L20: */
}
if (!isThrowGotoL30Statement) {
m = n;
}
// L30:
l = l1;
lsv = l;
lend = m;
lendsv = lend;
l1 = m + 1;
} while (lend == l);
/* Scale submatrix in rows and columns L to LEND */
anorm = dlanst("I", lend - l + 1, l + pointer_d, l + pointer_e); //$NON-NLS-1$
iscale = 0;
if (anorm.isGreaterThan(ssfmax)) {
iscale = 1;
info = dlascl_d("G", 0, 0, anorm, ssfmax, lend - l + 1, 1, l + pointer_d, n); //$NON-NLS-1$
info = dlascl_e("G", 0, 0, anorm, ssfmax, lend - l, 1, l + pointer_e, n); //$NON-NLS-1$
} else if (anorm.isLessThan(ssfmin)) {
iscale = 2;
info = dlascl_d("G", 0, 0, anorm, ssfmin, lend - l + 1, 1, l + pointer_d, n); //$NON-NLS-1$
info = dlascl_e("G", 0, 0, anorm, ssfmin, lend - l, 1, l + pointer_e, n); //$NON-NLS-1$
}
for (int i = l; i <= lend - 1; ++i) {
RS d1 = this.e[i + pointer_e];
this.e[i + pointer_e] = d1.multiply(d1);
/* L40: */
}
/* Choose between QL and QR iteration */
RS d1 = this.d[lend + pointer_d];
RS d2 = this.d[l + pointer_d];
if (d1.abs().isLessThan(d2.abs())) {
lend = lsv;
l = lendsv;
}
if (lend >= l) {
/*
* QL Iteration
* Look for small subdiagonal element.
*/
boolean isGotoL50 = false;
do {
isGotoL50 = false;
// L50:
int m = 0;
boolean isThrowGotoL70Statement = false;
if (l != lend) {
for (m = l; m <= lend - 1; ++m) {
d1 = this.d[m + pointer_d];
d2 = this.e[m + pointer_e];
if ((d2).abs().isLessThanOrEquals(eps2.multiply((d1.multiply(this.d[m + 1 + pointer_d])).abs()))) {
isThrowGotoL70Statement = true;
break;
}
/* L60: */
}
}
if (!isThrowGotoL70Statement) {
m = lend;
}
// L70:
if (m < lend) {
this.e[m + pointer_e] = unit.create(0);
}
RS p = this.d[l + pointer_d];
if (m == l) {
// L90:
this.d[l + pointer_d] = p;
++l;
if (l <= lend) {
isGotoL50 = true;
continue;
} else {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
}
/*
* If remaining matrix is 2 by 2, use DLAE2 to compute its eigenvalues.
*/
if (m == l + 1) {
RS rte = this.e[l + pointer_e].sqrt();
RS[] temp = Clapack.dlae2(this.d[l + pointer_d], rte, this.d[l + 1 + pointer_d]);
RS rt1 = temp[0];
RS rt2 = temp[1];
this.d[l + pointer_d] = rt1;
this.d[l + 1 + pointer_d] = rt2;
this.e[l + pointer_e] = unit.create(0);
l += 2;
if (l <= lend) {
isGotoL50 = true;
continue;
} else {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
}
if (jtot == nmaxit) {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
++jtot;
/* Form shift. */
RS rte = this.e[l + pointer_e].sqrt();
sigma = (this.d[l + 1 + pointer_d].subtract(p)).divide(rte.multiply(2));
RS c_b32 = unit.create(1);
RS r = Clapack.dlapy2(sigma, c_b32);
sigma = p.subtract(rte.divide((sigma.add(LibF77.d_sign(r, sigma)))));
RS c = unit.create(1);
RS s = unit.create(0);
RS gamma = this.d[m + pointer_d].subtract(sigma);
p = gamma.multiply(gamma);
/* Inner loop */
for (int i = m - 1; i >= l; --i) {
RS bb = this.e[i + pointer_e];
r = p.add(bb);
if (i != m - 1) {
this.e[i + 1 + pointer_e] = s.multiply(r);
}
RS oldc = c;
c = p.divide(r);
s = bb.divide(r);
RS oldgam = gamma;
RS alpha = this.d[i + pointer_d];
gamma = (c.multiply(alpha.subtract(sigma))).subtract(s.multiply(oldgam));
this.d[i + 1 + pointer_d] = oldgam.add(alpha.subtract(gamma));
if (!c.equals(unit.create(0))) {
p = gamma.multiply(gamma).divide(c);
} else {
p = oldc.multiply(bb);
}
/* L80: */
}
this.e[l + pointer_e] = s.multiply(p);
this.d[l + pointer_d] = sigma.add(gamma);
isGotoL50 = true;
} while (isGotoL50);
/* Eigenvalue found. */
// L90:
// d__[l+pointer_d__] = p;
// ++l;
// if (l <= lend) {
// goto L50;
// }
// goto L150;
// -------------------たぶんここまでQL Iteration--------
} else {
/*
* QR Iteration
* Look for small superdiagonal element.
*/
boolean isGotoL100 = false;
do {
isGotoL100 = false;
// L100:
boolean isGotoL120 = false;
int m;
for (m = l; m >= lend + 1; --m) {
d1 = this.d[m + pointer_d];
d2 = this.e[m - 1 + pointer_e];
if ((d2).abs().isLessThanOrEquals(eps2.multiply((d1.multiply(this.d[m - 1 + pointer_d])).abs()))) {
isGotoL120 = true;
break;
}
/* L110: */
}
if (!isGotoL120) {
m = lend;
}
// L120:
if (m > lend) {
this.e[m - 1 + pointer_e] = unit.create(0);
}
RS p = this.d[l + pointer_d];
if (m == l) {
// L140:
this.d[l + pointer_d] = p;
--l;
if (l >= lend) {
isGotoL100 = true;
continue;
} else {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
}
/*
* If remaining matrix is 2 by 2, use DLAE2 to compute its eigenvalues.
*/
if (m == l - 1) {
RS rte = this.e[l - 1 + pointer_e].sqrt();
RS[] temp = Clapack.dlae2(this.d[l + pointer_d], rte, this.d[l - 1 + pointer_d]);
RS rt1 = temp[0];
RS rt2 = temp[1];
this.d[l + pointer_d] = rt1;
this.d[l - 1 + pointer_d] = rt2;
this.e[l - 1 + pointer_e] = unit.create(0);
l += -2;
if (l >= lend) {
isGotoL100 = true;
continue;
} else {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
}
if (jtot == nmaxit) {
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
break;
}
++jtot;
/* Form shift. */
RS rte = this.e[l - 1 + pointer_e].sqrt();
sigma = (this.d[l - 1 + pointer_d].subtract(p)).divide(rte.multiply(2d));
RS c_b32 = unit.create(1);
RS r = Clapack.dlapy2(sigma, c_b32);
sigma = (p.subtract(rte)).divide((sigma.add(LibF77.d_sign(r, sigma))));
RS c = unit.create(1);
RS s = unit.create(0);
RS gamma = this.d[m + pointer_d].subtract(sigma);
p = gamma.multiply(gamma);
/* Inner loop */
for (int i = m; i <= l - 1; ++i) {
RS bb = this.e[i + pointer_e];
r = p.add(bb);
if (i != m) {
this.e[i - 1 + pointer_e] = s.multiply(r);
}
RS oldc = c;
c = p.divide(r);
s = bb.divide(r);
RS oldgam = gamma;
RS alpha = this.d[i + 1 + pointer_d];
gamma = (c.multiply(alpha.subtract(sigma))).subtract(s.multiply(oldgam));
this.d[i + pointer_d] = oldgam.add(alpha.subtract(gamma));
if (!c.equals(unit.create(0))) {
p = gamma.multiply(gamma).divide(c);
} else {
p = oldc.multiply(bb);
}
/* L130: */
}
this.e[l - 1 + pointer_e] = s.multiply(p);
this.d[l + pointer_d] = sigma.add(gamma);
isGotoL100 = true;
} while (isGotoL100);
/* Eigenvalue found. */
// L140:
// d__[l+pointer_d__] = p;
//
// --l;
// if (l >= lend) {
// goto L100;
// }
// goto L150;
// ---------ここまっでQR Iteration--------------
}
/* Undo scaling if necessary */
info = gotoL150(pointer_d, n, iscale, lendsv, lsv, anorm, ssfmin, ssfmax);
// L150:if (iscale == 1) {
// i__1 = lendsv - lsv + 1;
// dlascl_d("G", c__0, c__0, ssfmax, anorm, i__1, c__1, lsv + pointer_d__,
// n, info);
// }
// if (iscale == 2) {
// i__1 = lendsv - lsv + 1;
// dlascl_d("G", c__0, c__0, ssfmin, anorm, i__1, c__1, lsv + pointer_d__,
// n, info);
// }
/*
* Check for no convergence to an eigenvalue after a total of NMAXIT iterations.
*/
} while (jtot < nmaxit);
for (int i = 1; i <= n - 1; ++i) {
if (!this.e[i + pointer_e].equals(unit.create(0))) {
info = info + 1;
}
/* L160: */
}
return info;
/* Sort eigenvalues in increasing order. */
// L170:
// dlasrt("I", n, 1+pointer_d__, info);
// L180:
// return 0;
}
/**
* @param pointer_d pointer_d
* @param n n
* @param iscale iscale
* @param lendsv lendsv
* @param lsv lsv
* @param anorm anorm
* @param ssfmin ssfmin
* @param ssfmax ssfmax
* @return result
*/
private int gotoL150(int pointer_d, int n, int iscale, int lendsv, int lsv, RS anorm, RS ssfmin, RS ssfmax) {
// L150:
int info = 0;
if (iscale == 1) {
info = dlascl_d("G", 0, 0, ssfmax, anorm, lendsv - lsv + 1, 1, lsv + pointer_d, n); //$NON-NLS-1$
}
if (iscale == 2) {
info = dlascl_d("G", 0, 0, ssfmin, anorm, lendsv - lsv + 1, 1, lsv + pointer_d, n); //$NON-NLS-1$
}
return info;
}
/**
* @param type type
* @param kl kl
* @param ku ku
* @param cfrom cfrom
* @param cto cto
* @param m m
* @param n2 n2
* @param aIndex aIndex
* @param lda lda
* @return result
*/
private int dlascl_d(String type, int kl, int ku, RS cfrom, RS cto, int m, int n2, int aIndex, int lda) {
final RS unit = cfrom.createUnit();
RS[] aTemp = unit.createArray(this.d.length - aIndex);
System.arraycopy(this.d, aIndex, aTemp, 0, aTemp.length);
int info = Clapack.dlascl(type, kl, ku, cfrom, cto, m, n2, aTemp, lda);
System.arraycopy(aTemp, 0, this.d, aIndex, aTemp.length);
return info;
}
/**
* @param type type
* @param kl kl
* @param ku ku
* @param cfrom cfrom
* @param cto cto
* @param m m
* @param n2 n2
* @param aIndex aIndex
* @param lda lda
* @return result
*/
private int dlascl_e(String type, int kl, int ku, RS cfrom, RS cto, int m, int n2, int aIndex, int lda) {
final RS unit = cfrom.createUnit();
RS[] aTemp = unit.createArray(this.e.length - aIndex);
System.arraycopy(this.e, aIndex, aTemp, 0, aTemp.length);
int info = Clapack.dlascl(type, kl, ku, cfrom, cto, m, n2, aTemp, lda);
System.arraycopy(aTemp, 0, this.e, aIndex, aTemp.length);
return info;
}
/**
* @param id id
* @param n n
* @param dIndex d input/output
* @return result
*/
private int dlasrt(String id, int n, int dIndex) {
final RS unit = this.d[0].createUnit();
RS[] dTemp = unit.createArray(this.d.length - dIndex);
System.arraycopy(this.d, dIndex, dTemp, 0, dTemp.length);
int info = Clapack.dlasrt(id, n, dTemp);
System.arraycopy(dTemp, 0, this.d, dIndex, dTemp.length);
return info;
}
/**
* @param norm norm
* @param n n
* @param dIndex dIndex
* @param eIndex eIndex
* @return result
*/
private RS dlanst(String norm, int n, int dIndex, int eIndex) {
final RS unit = this.d[0].createUnit();
RS[] dTemp = unit.createArray(this.d.length - dIndex);
RS[] eTemp = unit.createArray(this.e.length - eIndex);
System.arraycopy(this.d, dIndex, dTemp, 0, dTemp.length);
System.arraycopy(this.e, eIndex, eTemp, 0, eTemp.length);
return Clapack.dlanst(norm, n, this.d, this.e);
}
}