Dsytd2.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 Dsytd2<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[] tau;
/* -- LAPACK 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
=======
DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
form T by an orthogonal similarity transformation: Q' * A * Q = T.
Arguments
=========
UPLO (input) CHARACTER*1
Specifies whether the upper or lower triangular part of the
symmetric matrix A is stored:
= 'U': Upper triangular
= 'L': Lower triangular
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, 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.
On exit, if UPLO = 'U', the diagonal and first superdiagonal
of A are overwritten by the corresponding elements of the
tridiagonal matrix T, and the elements above the first
superdiagonal, with the array TAU, represent the orthogonal
matrix Q as a product of elementary reflectors; if UPLO
= 'L', the diagonal and first subdiagonal of A are over-
written by the corresponding elements of the tridiagonal
matrix T, and the elements below the first subdiagonal, with
the array TAU, represent the orthogonal matrix Q as a product
of elementary reflectors. See Further Details.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
D (output) DOUBLE PRECISION array, dimension (N)
The diagonal elements of the tridiagonal matrix T:
D(i) = A(i,i).
RS (output) DOUBLE PRECISION array, dimension (N-1)
The off-diagonal elements of the tridiagonal matrix T:
E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
TAU (output) DOUBLE PRECISION array, dimension (N-1)
The scalar factors of the elementary reflectors (see Further
Details).
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value.
Further Details
===============
If UPLO = 'U', the matrix Q is represented as a product of elementary
reflectors
Q = H(n-1) . . . H(2) H(1).
Each H(i) has the form
H(i) = I - tau * v * v'
where tau is a real scalar, and v is a real vector with
v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
A(1:i-1,i+1), and tau in TAU(i).
If UPLO = 'L', the matrix Q is represented as a product of elementary
reflectors
Q = H(1) H(2) . . . H(n-1).
Each H(i) has the form
H(i) = I - tau * v * v'
where tau is a real scalar, and v is a real vector with
v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
and tau in TAU(i).
The contents of A on exit are illustrated by the following examples
with n = 5:
if UPLO = 'U': if UPLO = 'L':
( d e v2 v3 v4 ) ( d )
( d e v3 v4 ) ( e d )
( d e v4 ) ( v1 e d )
( d e ) ( v1 v2 e d )
( d ) ( v1 v2 v3 e d )
where d and e denote diagonal and off-diagonal elements of T, and vi
denotes an element of the vector defining H(i).
=====================================================================
*/
/**
* @param uplo uplo
* @param n n
* @param ai ai
* @param lda lda
* @param d d
* @param e e
* @param tauii tauii
* @return result
*/
int execute(String uplo, int n, RS[] ai, int lda, RS[] d, RS[] e, RS[] tauii) {
final RS unit = ai[0].createUnit();
this.a = ai;
this.tau = tauii;
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
int pointer_d = -1;
int pointer_e = -1;
int pointer_tau = -1;
int info = 0;
boolean upper = BLAS.lsame(uplo, "U"); //$NON-NLS-1$
if (!upper && !BLAS.lsame(uplo, "L")) { //$NON-NLS-1$
info = -1;
} else if (n < 0) {
info = -2;
} else if (lda < Math.max(1, n)) {
info = -4;
}
if (info != 0) {
BLAS.xerbla("DSYTD2", info); //$NON-NLS-1$
return info;
}
/* Quick return if possible */
if (n <= 0) {
return info;
}
RS taui = unit.create(0);
if (upper) {
/* Reduce the upper triangle of A */
for (int i = n - 1; i >= 1; --i) {
/* Generate elementary reflector H(i) = I - tau * v * v' to annihilate A(1:i-1,i+1) */
RS alpha1 = this.a[(i + 1) * a_dim1 + i + pointer_a];
RS[] ans = dlarfg(i, alpha1, (i + 1) * a_dim1 + 1 + pointer_a, 1, taui);
this.a[(i + 1) * a_dim1 + i + pointer_a] = ans[0];
taui = ans[1];
e[pointer_e + i] = this.a[(i + 1) * a_dim1 + i + pointer_a];
//e[3]の値
if (!taui.isZero()) {
/* Apply H(i) from both sides to A(1:i,1:i) */
this.a[(i + 1) * a_dim1 + i + pointer_a] = unit.create(1);
/* Compute x := tau * A * v storing x in TAU(1:i) */
RS c_b8 = unit.create(0);
dsymv(uplo, i, taui, pointer_a + a_offset, lda, (i + 1) * a_dim1 + 1 + pointer_a, 1, c_b8, pointer_tau + 1, 1);
/* Compute w := x - 1/2 * tau * (x'*v) * v */
RS alpha = taui.divide(2).unaryMinus().multiply(ddot(i, pointer_tau + 1, 1, (i + 1) * a_dim1 + 1 + pointer_a, 1));
daxpy(i, alpha, (i + 1) * a_dim1 + 1 + pointer_a, 1, pointer_tau + 1, 1);
/* Apply the transformation as a rank-2 update:A := A - v * w' - w * v' */
RS c_b14 = unit.create(-1);
dsyr2(uplo, i, c_b14, (i + 1) * a_dim1 + 1 + pointer_a, 1, pointer_tau + 1, 1, pointer_a + a_offset, lda);
this.a[(i + 1) * a_dim1 + i + pointer_a] = e[pointer_e + i];
}
d[pointer_d + i + 1] = this.a[(i + 1) * a_dim1 + i + 1 + pointer_a];
this.tau[pointer_tau + i] = taui;
/* L10: */
}
d[pointer_d + 1] = this.a[a_dim1 + 1 + pointer_a];
} else {
/* Reduce the lower triangle of A */
for (int i = 1; i <= n - 1; ++i) {
/* Generate elementary reflector H(i) = I - tau * v * v' to annihilate A(i+2:n,i) */
RS alpha1 = this.a[i * a_dim1 + i + 1 + pointer_a];
RS[] ans = dlarfg(n - i, alpha1, i * a_dim1 + Math.min(i + 2, n) + pointer_a, 1, taui);
this.a[i * a_dim1 + i + 1 + pointer_a] = ans[0];
taui = ans[1];
e[pointer_e + i] = this.a[i * a_dim1 + i + 1 + pointer_a];
if (!taui.isZero()) {
// Apply H(i) from both sides to A(i+1:n,i+1:n)
this.a[i * a_dim1 + i + 1 + pointer_a] = unit.create(1);
// Compute x := tau * A * v storing y in TAU(i:n-1)
RS c_b8 = unit.create(0);
dsymv(uplo, n - i, taui, (i + 1) * a_dim1 + i + 1 + pointer_a, lda, i * a_dim1 + i + 1 + pointer_a, 1, c_b8, pointer_tau + i, 1);
// Compute w := x - 1/2 * tau * (x'*v) * v
RS alpha = taui.divide(2).unaryMinus().multiply(ddot(n - i, pointer_tau + i, 1, i * a_dim1 + i + 1 + pointer_a, 1));
daxpy(n - i, alpha, i * a_dim1 + i + 1 + pointer_a, 1, pointer_tau + i, 1);
// Apply the transformation as a rank-2 update: A := A - v * w' - w * v'
RS c_b14 = unit.create(-1);
dsyr2(uplo, n - i, c_b14, i * a_dim1 + i + 1 + pointer_a, 1, pointer_tau + i, 1, (i + 1) * a_dim1 + i + 1 + pointer_a, lda);
this.a[i * a_dim1 + i + 1 + pointer_a] = e[pointer_e + i];
}
d[pointer_d + i] = this.a[i * a_dim1 + i + pointer_a];
this.tau[pointer_tau + i] = taui;
/* L20: */
}
d[pointer_d + n] = this.a[n * a_dim1 + n + pointer_a];
}
return info;
}
/**
* @param n n
* @param dxIndex tau
* @param incx indx
* @param dyIndex a
* @param incy incy
* @return result
*/
private RS ddot(int n, int dxIndex, int incx, int dyIndex, int incy) {
final RS unit = this.a[0].createUnit();
RS[] dx = unit.createArray(this.tau.length - dxIndex);
RS[] dy = unit.createArray(this.a.length - dyIndex);
System.arraycopy(this.tau, dxIndex, dx, 0, dx.length);
System.arraycopy(this.a, dyIndex, dy, 0, dy.length);
return BLAS.ddot(n, dx, incx, dy, incy);
}
/**
* @param uplo uplo
* @param n n
* @param alpha alpha
* @param xIndex a
* @param incx indx
* @param yIndex tau
* @param incy incy
* @param aIndex a
* @param lda lda
*/
private void dsyr2(String uplo, int n, RS alpha, int xIndex, int incx, int yIndex, int incy, int aIndex, int lda) {
final RS unit = alpha.createUnit();
RS[] x = unit.createArray(this.a.length - xIndex);
RS[] y = unit.createArray(this.tau.length - yIndex);
RS[] aTemp = unit.createArray(this.a.length - aIndex);
System.arraycopy(this.a, xIndex, x, 0, x.length);
System.arraycopy(this.tau, yIndex, y, 0, y.length);
System.arraycopy(this.a, aIndex, aTemp, 0, aTemp.length);
BLAS.dsyr2(uplo, n, alpha, x, incx, y, incy, aTemp, lda);
System.arraycopy(aTemp, 0, this.a, aIndex, aTemp.length);
}
/**
* @param n n
* @param da da
* @param dxIndex dxIndex
* @param incx incx
* @param dyIndex dyIndex
* @param incy incy
*/
private void daxpy(int n, RS da, int dxIndex, int incx, int dyIndex, int incy) {
final RS unit = da.createUnit();
RS[] dx = unit.createArray(this.a.length - dxIndex);
RS[] dy = unit.createArray(this.tau.length - dyIndex);
System.arraycopy(this.a, dxIndex, dx, 0, dx.length);
System.arraycopy(this.tau, dyIndex, dy, 0, dy.length);
BLAS.daxpy(n, da, dx, incx, dy, incy);
System.arraycopy(dy, 0, this.tau, dyIndex, dy.length);
}
/**
* @param uplo uplo
* @param n n
* @param alpha alpha
* @param aIndex aIndex
* @param lda lda
* @param xIndex xIndex
* @param incx incx
* @param beta beta
* @param yIndex tau
* @param incy incy
*/
private void dsymv(String uplo, int n, RS alpha, int aIndex, int lda, int xIndex, int incx, RS beta, int yIndex, int incy) {
final RS unit = alpha.createUnit();
RS[] aTemp = unit.createArray(this.a.length - aIndex);
RS[] x = unit.createArray(this.a.length - xIndex);
RS[] y = unit.createArray(this.tau.length - yIndex);
System.arraycopy(this.a, aIndex, aTemp, 0, aTemp.length);
System.arraycopy(this.a, xIndex, x, 0, x.length);
System.arraycopy(this.tau, yIndex, y, 0, y.length);
BLAS.dsymv(uplo, n, alpha, aTemp, lda, x, incx, beta, y, incy);
System.arraycopy(y, 0, this.tau, yIndex, y.length);
}
/**
* @param n input
* @param alpha input/output
* @param xIndex input/output array"a"
* @param incx input
* @param tau output
* @return result
*/
private RS[] dlarfg(int n, RS alpha, int xIndex, int incx, RS tau) {
final RS unit = alpha.createUnit();
RS[] x = unit.createArray(this.a.length - xIndex);
System.arraycopy(this.a, xIndex, x, 0, x.length);
RS[] ans = Clapack.dlarfg(n, alpha, x, incx, tau);
System.arraycopy(x, 0, this.a, xIndex, x.length);
return ans;
}
}