Dpotrf.java
package org.mklab.sdpj.algebra;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
/**
* @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 Dpotrf<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 T POTRF_NONZERO;
// /** */
// private T POTRF_ASSIGN;
// /** */
// private T POTRF_LIMIT;
//
// /**
// * Set parameters.
// * @param unit 単位
// */
// private void setParameters(T unit) {
// if (unit instanceof MPFloat) {
//
// this.POTRF_NONZERO = unit.create(10).power(-14); //-100
// this.POTRF_ASSIGN = unit.create(10).power(100); //100
// this.POTRF_LIMIT = unit.create(10).power(-6).unaryMinus(); //-100
// }
// if (unit instanceof DoubleNumber) {
// this.POTRF_NONZERO = unit.create(10).power(-14);//-10 //-7
// this.POTRF_ASSIGN = unit.create(10).power(52);//52 //50
// this.POTRF_LIMIT = unit.create(10).power(-7).unaryMinus(); //-5 //-6
// }
// }
//
// /**
// * @param n
// * @param A
// * @param lda
// * @return
// */
// int rATL_dpotrfL(int n, T[] A, int lda) {
// final T unit = A[0].createUnit();
// setParameters(unit);
//
// if (n > 4) {
// int nLeft = n >> 1;
// int nRight = n - nLeft;
// int ierr = rATL_dpotrfL(nLeft, A, lda);
//
// /*
// * The value of ierr is maybe 0,1,2,3 or 4.
// * condition : if(ierr != 1) ,which is NO WARRANTY.
// * memo is made by wu,2010/09/05
// */
// if (ierr != 1) {
// int pointerAr = 0 + nLeft;
// T[] Ar = unit.createArray(A.length - pointerAr);
// System.arraycopy(A, pointerAr, Ar, 0, A.length - pointerAr); //A[ArPointer];
//
// int pointerAn = pointerAr + lda * nLeft;
// T[] An = unit.createArray(A.length - pointerAn);
//
// BLAS.dtrsm("R", "L", "T", "N", nRight, nLeft, unit.create(1), A, lda, Ar, lda); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
// System.arraycopy(Ar, 0, A, pointerAr, Ar.length);
// System.arraycopy(A, pointerAn, An, 0, A.length - pointerAn);
// BLAS.dsyrk("L", "N", nRight, nLeft, unit.create(-1), Ar, lda, unit.create(1), An, lda); //$NON-NLS-1$ //$NON-NLS-2$
// ierr = rATL_dpotrfL(nRight, An, lda);
// System.arraycopy(An, 0, A, pointerAn, An.length);
// if (ierr == 1) {
// return (ierr + nLeft);
// }
// } else {
// return (ierr);
// }
// } else if (n == 4) {
// return (potrf4(A, lda));
// } else if (n == 3) {
// return (potrf3(A, lda));
// } else if (n == 2) {
// return (potrf2(A, lda));
// } else if (n == 1) {
// if (A[0].isLessThan(this.POTRF_LIMIT)) {
// return 1;
// }
// if (A[0].isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// A[0] = this.POTRF_ASSIGN;
// }
// A[0] = A[0].sqrt();
// }
// return 0;
// }
//
// /**
// * @param A
// * @param n
// * @return
// */
// private int potrf4(T[] A, final int n) {
// final T unit = A[0].createUnit();
//
// int pointerA1 = 0 + n + 1;
// T[] A1 = unit.createArray(A.length - pointerA1);// = A+n+1;
// System.arraycopy(A, pointerA1, A1, 0, A.length - pointerA1);
//
// int pointerA2 = pointerA1 + n + 1;
// T[] A2 = unit.createArray(A.length - pointerA2);// =
// // A1+n+1;
// System.arraycopy(A, pointerA2, A2, 0, A.length - pointerA2);
//
// int pointerA3 = pointerA2 + n + 1;
// T[] A3 = unit.createArray(A.length - pointerA3);
// System.arraycopy(A, pointerA3, A3, 0, A.length - pointerA3);
//
// T L11 = A[0];
// T L21 = A[1];
// T L22 = A1[0];
// T L31 = A[2];
// T L32 = A1[1];
// T L33 = A2[0];
// T L41 = A[3];
// T L42 = A1[2];
// T L43 = A2[1];
// T L44 = A3[0];
//
// if (L11.isLessThan(this.POTRF_LIMIT)) {
// return 1;
// }
// if (L11.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L11 = this.POTRF_ASSIGN;
// }
//
// L11 = L11.sqrt();
// A[0] = L11;
// L11 = unit.create(1).divide(L11);
// L21 = L21.multiply(L11);
// L31 = L31.multiply(L11);
// L41 = L41.multiply(L11);
//
// L22 = L22.subtract(L21.multiply(L21));
// if (L22.isLessThan(this.POTRF_LIMIT)) {
// return 2;
// }
// if (L22.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L22 = this.POTRF_ASSIGN;
// }
// L22 = L22.sqrt();
// A1[0] = L22;
// L22 = unit.create(1).divide(L22);
// L32 = (L32.subtract(L31.multiply(L21))).multiply(L22);
// L42 = (L42.subtract(L41.multiply(L21))).multiply(L22);
//
// L33 = L33.subtract((L31.multiply(L31)).add(L32.multiply(L32)));
// // L33 -= L31*L31 + L32*L32;
//
// if (L33.isLessThan(this.POTRF_LIMIT)) {
// return 3;
// }
// if (L33.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L33 = this.POTRF_ASSIGN;
// }
//
// L33 = L33.sqrt();
// A2[0] = L33;
// L43 = (L43.subtract(L41.multiply(L31)).subtract(L42.multiply(L32))).divide(L33);
// // L43 = (L43-L41*L31-L42*L32)/L33;
//
// L44 = L44.subtract(L41.power(2).add(L42.power(2)).add(L43.power(2)));
// // L44 -= L41*L41 + L42*L42 + L43*L43;
//
// if (L44.isLessThan(this.POTRF_LIMIT)) {
// return 4;
// }
//
// if (L44.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L44 = this.POTRF_ASSIGN;
// }
//
// A3[0] = L44.sqrt();
// A[1] = L21;
// A[2] = L31;
// A1[1] = L32;
// A[3] = L41;
// A1[2] = L42;
// A2[1] = L43;
//
// System.arraycopy(A1, 0, A, pointerA1, A1.length);
// System.arraycopy(A2, 0, A, pointerA2, A2.length);
// System.arraycopy(A3, 0, A, pointerA3, A3.length);
//
// return 0;
// }
//
// /**
// * @param A
// * @param n
// * @return
// */
// private int potrf3(T[] A, final int n) {
// final T unit = A[0].createUnit();
//
// int pointerA1 = 0 + n + 1;
// T[] A1 = unit.createArray(A.length - pointerA1);// A+n+1;
// System.arraycopy(A, pointerA1, A1, 0, A.length - pointerA1);
//
// int pointerA2 = pointerA1 + n + 1;
// T[] A2 = unit.createArray(A.length - pointerA2);// A1+n+1;
// System.arraycopy(A, pointerA2, A2, 0, A.length - pointerA2);
//
// T L11 = A[0];
// T L21 = A[1];
// T L22 = A1[0];
// T L31 = A[2];
// T L32 = A1[1];
// T L33 = A2[0];
//
// if (L11.isLessThan(this.POTRF_LIMIT)) {
// return 1;
// }
// if (L11.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L11 = this.POTRF_ASSIGN;
// }
//
// L11 = L11.sqrt();
// A[0] = L11;
// L11 = unit.create(1).divide(L11);
// L21 = L21.multiply(L11);
// L31 = L31.multiply(L11);
// L22 = L22.subtract(L21.power(2));
//
// if (L22.isLessThan(this.POTRF_LIMIT)) {
// return 2;
// }
//
// if (L22.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L22 = this.POTRF_ASSIGN;
// }
//
// L22 = L22.sqrt();
// L32 = (L32.subtract(L31.multiply(L21))).divide(L22);
// L33 = L33.subtract(L31.power(2).add(L32.power(2)));
//
// if (L33.isLessThan(this.POTRF_LIMIT)) {
// return 3;
// }
//
// if (L33.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L33 = this.POTRF_ASSIGN;
// }
//
// A2[0] = L33.sqrt();
// A[1] = L21;
// A1[0] = L22;
// A[2] = L31;
// A1[1] = L32;
//
// System.arraycopy(A1, 0, A, pointerA1, A1.length);
// System.arraycopy(A2, 0, A, pointerA2, A2.length);
//
// return 0;
// }
//
// /**
// * @param A
// * @param n
// * @return
// */
// private int potrf2(T[] A, final int n) {
// final T unit = A[0].createUnit();
//
// int pointerA1 = 0 + n + 1;
// T[] A1 = unit.createArray(A.length - pointerA1);// A+n+1;
// System.arraycopy(A, pointerA1, A1, 0, A.length - pointerA1);
//
// T L11 = A[0];
// T L21 = A[1];
// T L22 = A1[0];
//
// if (L11.isLessThan(this.POTRF_LIMIT)) {
// return 1;
// }
//
// if (L11.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L11 = this.POTRF_ASSIGN;
// }
//
// L11 = L11.sqrt();
// A[0] = L11;
// L21 = L21.divide(L11);
// L22 = L22.subtract(L21.power(2));
//
// if (L22.isLessThan(this.POTRF_LIMIT)) {
// return 2;
// }
//
// if (L22.isLessThan(this.POTRF_NONZERO)) {
// CHOLESKY_ADJUST();
// L22 = this.POTRF_ASSIGN;
// }
//
// A[0] = L11;
// A[1] = L21;
// A1[0] = L22.sqrt();
// System.arraycopy(A1, 0, A, pointerA1, A1.length);
// return 0;
// }
//
// /**
// *
// */
// private void CHOLESKY_ADJUST() {
// // nothing to do
// }
}