Dpotf2.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 Dpotf2<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 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
=======
DPOTF2 computes the Cholesky factorization of a real symmetric
positive definite matrix A.
The factorization has the form
A = U' * U , if UPLO = 'U', or
A = L * L', if UPLO = 'L',
where U is an upper triangular matrix and L is lower triangular.
This is the unblocked version of the algorithm, calling Level 2 BLAS.
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 INFO = 0, the factor U or L from the Cholesky
factorization A = U'*U or A = L*L'.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -k, the k-th argument had an illegal value
> 0: if INFO = k, the leading minor of order k is not
positive definite, and the factorization could not be
completed.
=====================================================================
*/
/**
* @param uplo uplo
* @param n n
* @param a a
* @param lda lda
* @return result
*/
int execute(String uplo, int n, RS[] a, int lda) {
final RS unit = a[0].createUnit();
RS c_b10 = unit.create(-1);
RS c_b12 = unit.create(1);
int a_dim1 = lda;
int a_offset = 1 + a_dim1 * 1;
int pointer_a = -a_offset;
/* Function Body */
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("DPOTF2", -info); //$NON-NLS-1$
return info;
}
/* Quick return if possible */
if (n == 0) {
return info;
}
if (upper) {
/* Compute the Cholesky factorization A = U'U. */
for (int j = 1; j <= n; ++j) {
/* Compute U(J,J) and test for non-positive-definiteness. */
RS[] a_1jTemp = unit.createArray(a.length - (j * a_dim1 + 1 + pointer_a));
System.arraycopy(a, j * a_dim1 + 1 + pointer_a, a_1jTemp, 0, a_1jTemp.length);
RS ajj = a[j * a_dim1 + j + pointer_a].subtract(BLAS.ddot(j - 1, a_1jTemp, 1, a_1jTemp, 1));
System.arraycopy(a_1jTemp, 0, a, j * a_dim1 + 1 + pointer_a, a_1jTemp.length);
if (ajj.isLessThanOrEquals(0)) {
a[j * a_dim1 + j + pointer_a] = ajj;
info = j;
return info;
}
ajj = ajj.sqrt();
a[j * a_dim1 + j + pointer_a] = ajj;
/* Compute elements J+1:N of row J. */
if (j < n) {
RS[] a_1j1Temp = unit.createArray(a.length - ((j + 1) * a_dim1 + 1 + pointer_a));
RS[] a_jj1Temp = unit.createArray(a.length - ((j + 1) * a_dim1 + j + pointer_a));
System.arraycopy(a, j * a_dim1 + 1 + pointer_a, a_1jTemp, 0, a_1jTemp.length);
System.arraycopy(a, (j + 1) * a_dim1 + 1 + pointer_a, a_1j1Temp, 0, a_1j1Temp.length);
System.arraycopy(a, (j + 1) * a_dim1 + j + pointer_a, a_jj1Temp, 0, a_jj1Temp.length);
BLAS.dgemv("Transpose", j - 1, n - j, c_b10, a_1j1Temp, lda, a_1jTemp, 1, c_b12, a_jj1Temp, lda); //$NON-NLS-1$
RS d1 = unit.create(1).divide(ajj);
BLAS.dscal(n - j, d1, a_jj1Temp, lda);
System.arraycopy(a_1jTemp, 0, a, j * a_dim1 + 1 + pointer_a, a_1jTemp.length);
System.arraycopy(a_1j1Temp, 0, a, (j + 1) * a_dim1 + 1 + pointer_a, a_1j1Temp.length);
System.arraycopy(a_jj1Temp, 0, a, (j + 1) * a_dim1 + j + pointer_a, a_jj1Temp.length);
}
}
} else {
/* Compute the Cholesky factorization A = LL'. */
for (int j = 1; j <= n; ++j) {
/* Compute L(J,J) and test for non-positive-definiteness. */
RS[] a_j1Temp = unit.createArray(a.length - (a_dim1 + j + pointer_a));
System.arraycopy(a, a_dim1 + j + pointer_a, a_j1Temp, 0, a_j1Temp.length);
RS ajj = a[j * a_dim1 + j + pointer_a].subtract(BLAS.ddot(j - 1, a_j1Temp, lda, a_j1Temp, lda));
System.arraycopy(a_j1Temp, 0, a, a_dim1 + j + pointer_a, a_j1Temp.length);
if (ajj.isLessThanOrEquals(0)) {
a[j * a_dim1 + j + pointer_a] = ajj;
info = j;
return info;
}
ajj = ajj.sqrt();
a[j * a_dim1 + j + pointer_a] = ajj;
/* Compute elements J+1:N of column J. */
if (j < n) {
RS[] a_j11Temp = unit.createArray(a.length - (a_dim1 + j + 1 + pointer_a));
RS[] a_j1jTemp = unit.createArray(a.length - (j * a_dim1 + j + 1 + pointer_a));
a_j1Temp = unit.createArray(a.length - (a_dim1 + j + pointer_a));
System.arraycopy(a, a_dim1 + j + 1 + pointer_a, a_j11Temp, 0, a_j11Temp.length);
System.arraycopy(a, j * a_dim1 + j + 1 + pointer_a, a_j1jTemp, 0, a_j1jTemp.length);
System.arraycopy(a, a_dim1 + j + pointer_a, a_j1Temp, 0, a_j1Temp.length);
BLAS.dgemv("No transpose", n - j, j - 1, c_b10, a_j11Temp, lda, a_j1Temp, lda, c_b12, a_j1jTemp, 1); //$NON-NLS-1$
System.arraycopy(a_j1jTemp, 0, a, j * a_dim1 + j + 1 + pointer_a, a_j1jTemp.length);
System.arraycopy(a, a_dim1 + j + 1 + pointer_a, a_j11Temp, 0, a_j11Temp.length);
System.arraycopy(a, j * a_dim1 + j + 1 + pointer_a, a_j1jTemp, 0, a_j1jTemp.length);
System.arraycopy(a, a_dim1 + j + pointer_a, a_j1Temp, 0, a_j1Temp.length);
RS d1 = unit.create(1).divide(ajj);// 1. / ajj;
BLAS.dscal(n - j, d1, a_j1jTemp, 1);
System.arraycopy(a_j1jTemp, 0, a, j * a_dim1 + j + 1 + pointer_a, a_j1jTemp.length);
}
}
}
return info;
}
}