final int n = currentBest.getDimension();
final int npt = numberOfInterpolationPoints;
final int nptm = npt - n - 1;
// XXX Should probably be split into two arrays.
final ArrayRealVector work = new ArrayRealVector(npt + n);
double ztest = ZERO;
for (int k = 0; k < npt; k++) {
for (int j = 0; j < nptm; j++) {
// Computing MAX
ztest = Math.max(ztest, Math.abs(zMatrix.getEntry(k, j)));
}
}
ztest *= 1e-20;
// Apply the rotations that put zeros in the KNEW-th row of ZMAT.
for (int j = 1; j < nptm; j++) {
final double d1 = zMatrix.getEntry(knew, j);
if (Math.abs(d1) > ztest) {
// Computing 2nd power
final double d2 = zMatrix.getEntry(knew, 0);
// Computing 2nd power
final double d3 = zMatrix.getEntry(knew, j);
final double d4 = Math.sqrt(d2 * d2 + d3 * d3);
final double d5 = zMatrix.getEntry(knew, 0) / d4;
final double d6 = zMatrix.getEntry(knew, j) / d4;
for (int i = 0; i < npt; i++) {
final double d7 = d5 * zMatrix.getEntry(i, 0) + d6 * zMatrix.getEntry(i, j);
zMatrix.setEntry(i, j, d5 * zMatrix.getEntry(i, j) - d6 * zMatrix.getEntry(i, 0));
zMatrix.setEntry(i, 0, d7);
}
}
zMatrix.setEntry(knew, j, ZERO);
}
// Put the first NPT components of the KNEW-th column of HLAG into W,
// and calculate the parameters of the updating formula.
for (int i = 0; i < npt; i++) {
work.setEntry(i, zMatrix.getEntry(knew, 0) * zMatrix.getEntry(i, 0));
}
final double alpha = work.getEntry(knew);
final double tau = lagrangeValuesAtNewPoint.getEntry(knew);
lagrangeValuesAtNewPoint.setEntry(knew, lagrangeValuesAtNewPoint.getEntry(knew) - ONE);
// Complete the updating of ZMAT.
final double sqrtDenom = Math.sqrt(denom);
final double d1 = tau / sqrtDenom;
final double d2 = zMatrix.getEntry(knew, 0) / sqrtDenom;
for (int i = 0; i < npt; i++) {
zMatrix.setEntry(i, 0,
d1 * zMatrix.getEntry(i, 0) - d2 * lagrangeValuesAtNewPoint.getEntry(i));
}
// Finally, update the matrix BMAT.
for (int j = 0; j < n; j++) {
final int jp = npt + j;
work.setEntry(jp, bMatrix.getEntry(knew, j));
final double d3 = (alpha * lagrangeValuesAtNewPoint.getEntry(jp) - tau * work.getEntry(jp)) / denom;
final double d4 = (-beta * work.getEntry(jp) - tau * lagrangeValuesAtNewPoint.getEntry(jp)) / denom;
for (int i = 0; i <= jp; i++) {
bMatrix.setEntry(i, j,
bMatrix.getEntry(i, j) + d3 * lagrangeValuesAtNewPoint.getEntry(i) + d4 * work.getEntry(i));
if (i >= npt) {
bMatrix.setEntry(jp, (i - npt), bMatrix.getEntry(i, j));
}
}
}