[sdpb] 101/233: Made choleskyStabilizeThreshold a parameter; added (temporary) timers for a bunch of things
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:25 UTC 2017
This is an automated email from the git hooks/post-receive script.
thansen pushed a commit to branch master
in repository sdpb.
commit 3af2357837cf16ca5819638b2200d22eb24342ae
Author: David Simmons-Duffin <dsd at minerva.sns.ias.edu>
Date: Wed Oct 29 15:58:49 2014 -0400
Made choleskyStabilizeThreshold a parameter; added (temporary) timers for a bunch of things
---
src/BlockDiagonalMatrix.cpp | 8 ++++++--
src/BlockDiagonalMatrix.h | 3 ++-
src/Matrix.cpp | 4 ++--
src/Matrix.h | 2 +-
src/SDPSolver.cpp | 43 +++++++++++++++++++++++++++---------------
src/SDPSolver.h | 9 +++++++--
src/SDPSolverIO.cpp | 15 ++++++++++-----
src/main.cpp | 9 ++++++++-
src/mpack/Rgetrf.cpp | 12 ++++++++++++
src/mpack/Rpotf2Stabilized.cpp | 16 +++++++---------
src/mpack/RpotrfStabilized.cpp | 8 ++++----
src/mpack/mlapack_gmp.h | 4 ++--
src/tests.cpp | 3 ++-
13 files changed, 91 insertions(+), 45 deletions(-)
diff --git a/src/BlockDiagonalMatrix.cpp b/src/BlockDiagonalMatrix.cpp
index 8014680..b0d3760 100644
--- a/src/BlockDiagonalMatrix.cpp
+++ b/src/BlockDiagonalMatrix.cpp
@@ -100,10 +100,14 @@ void choleskyDecomposition(BlockDiagonalMatrix &A,
void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L,
vector<vector<Integer> > &schurStabilizeIndices,
- vector<vector<Real> > &schurStabilizeLambdas) {
+ vector<vector<Real> > &schurStabilizeLambdas,
+ const double stabilizeThreshold) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < A.blocks.size(); b++)
- choleskyDecompositionStabilized(A.blocks[b], L.blocks[b], schurStabilizeIndices[b], schurStabilizeLambdas[b]);
+ choleskyDecompositionStabilized(A.blocks[b], L.blocks[b],
+ schurStabilizeIndices[b],
+ schurStabilizeLambdas[b],
+ stabilizeThreshold);
}
void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
diff --git a/src/BlockDiagonalMatrix.h b/src/BlockDiagonalMatrix.h
index d26e7c5..9ea258e 100644
--- a/src/BlockDiagonalMatrix.h
+++ b/src/BlockDiagonalMatrix.h
@@ -104,7 +104,8 @@ void choleskyDecomposition(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L);
void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L,
vector<vector<Integer> > &schurStabilizeIndices,
- vector<vector<Real> > &schurStabilizeLambdas);
+ vector<vector<Real> > &schurStabilizeLambdas,
+ const double stabilizeThreshold);
void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky, BlockDiagonalMatrix &X);
void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, Matrix &B);
void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, Matrix &B);
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 0353bdd..ae108be 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -220,7 +220,7 @@ void choleskyDecomposition(Matrix &A, Matrix &L) {
// - stabilizeIndices: vector to hold indices where small diagonal elements were compensated
// - stabilizeLambdas: vector to hold values used to compensate small diagonal elements
//
-void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas) {
+void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas, const double stabilizeThreshold) {
int dim = A.rows;
assert(A.cols == dim);
assert(L.rows == dim);
@@ -229,7 +229,7 @@ void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stab
// Set lower-triangular part of L to cholesky decomposition
L.copyFrom(A);
Integer info;
- RpotrfStabilized("Lower", dim, &L.elements[0], dim, &info, stabilizeIndices, stabilizeLambdas);
+ RpotrfStabilized("Lower", dim, &L.elements[0], dim, &info, stabilizeIndices, stabilizeLambdas, stabilizeThreshold);
assert(info == 0);
// Set the upper-triangular part of the L to zero
diff --git a/src/Matrix.h b/src/Matrix.h
index 2b2cff8..6d206ac 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -167,7 +167,7 @@ void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b);
// - L : dim x dim lower-triangular matrix
void choleskyDecomposition(Matrix &A, Matrix &L);
-void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas);
+void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas, const double stabilizeThreshold);
void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb);
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 0f9de3e..e29ad4e 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -400,13 +400,17 @@ Real stepLength(BlockDiagonalMatrix &XCholesky,
}
void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
- const BlockDiagonalMatrix &BilinearPairingsY) {
-
+ const BlockDiagonalMatrix &BilinearPairingsY,
+ const Real &choleskyStabilizeThreshold) {
+ timers["schurblocks/cholesky"].resume();
computeSchurBlocks(sdp, BilinearPairingsXInv, BilinearPairingsY, SchurBlocks);
- choleskyDecompositionStabilized(SchurBlocks, SchurBlocksCholesky, schurStabilizeIndices, schurStabilizeLambdas);
-
+ choleskyDecompositionStabilized(SchurBlocks, SchurBlocksCholesky,
+ schurStabilizeIndices,
+ schurStabilizeLambdas,
+ cast2double(choleskyStabilizeThreshold));
+ timers["schurblocks/cholesky"].stop();
// SchurUpdateLowRank = {{- 1, 0}, {E, G}}
-
+ timers["make schur update"].resume();
SchurUpdateLowRank.copyFrom(sdp.FreeVarMatrix);
blockMatrixLowerTriangularSolve(SchurBlocksCholesky, SchurUpdateLowRank);
int updateColumns = SchurUpdateLowRank.cols;
@@ -447,7 +451,8 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
&stabilizeBlocks[b].elt(0, 0),
stabilizeBlocks[b].rows);
}
-
+ timers["make schur update"].stop();
+ timers["make Q"].resume();
// Q = SchurUpdateLowRank^T SchurUpdateLowRank - {{0,0},{0,1}}
Q.setRowsCols(updateColumns, updateColumns);
Q.setZero();
@@ -488,9 +493,12 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
for (int c = 0; c < SchurUpdateLowRank.cols; c++)
for (int r = SchurUpdateLowRank.cols; r < Q.rows; r++)
Q.elt(c,r) = Q.elt(r,c);
+ timers["make Q"].stop();
+ timers["LU of Q"].resume();
Qpivots.resize(Q.rows);
LUDecomposition(Q, Qpivots);
+ timers["LU of Q"].stop();
}
@@ -600,12 +608,14 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
if (timers["Solver runtime"].elapsed().wall >= parameters.maxRuntime * 1000000000LL)
return MaxRuntimeExceeded;
+ timers["cholesky"].resume();
choleskyDecomposition(X, XCholesky);
choleskyDecomposition(Y, YCholesky);
-
+ timers["cholesky"].stop();
+ timers["bilinear pairings"].resume();
computeInvBilinearPairingsWithCholesky(XCholesky, sdp.bilinearBases, bilinearPairingsWorkspace, BilinearPairingsXInv);
computeBilinearPairings(Y, sdp.bilinearBases, bilinearPairingsWorkspace, BilinearPairingsY);
-
+ timers["bilinear pairings"].stop();
// d_k = c_k - Tr(F_k Y) - (D y)_k
computeDualResidues(sdp, y, BilinearPairingsY, dualResidues);
@@ -631,33 +641,36 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
// functions for the current point
else if (iteration > parameters.maxIterations)
return MaxIterationsExceeded;
-
- initializeSchurComplementSolver(BilinearPairingsXInv, BilinearPairingsY);
+ timers["initialize schur solver"].resume();
+ initializeSchurComplementSolver(BilinearPairingsXInv, BilinearPairingsY, parameters.choleskyStabilizeThreshold);
+ timers["initialize schur solver"].stop();
Real mu = frobeniusProductSymmetric(X, Y)/X.dim;
if (mu > parameters.maxComplementarity)
return MaxComplementarityExceeded;
-
+ timers["predictor"].resume();
// Mehrotra predictor solution for (dx, dX, dY)
Real betaPredictor = predictorCenteringParameter(parameters, isPrimalFeasible && isDualFeasible);
computeSearchDirection(betaPredictor, mu, false);
-
+ timers["predictor"].stop();
+ timers["corrector"].resume();
// Mehrotra corrector solution for (dx, dX, dY)
Real betaCorrector = correctorCenteringParameter(parameters, X, dX, Y, dY, mu, isPrimalFeasible && isDualFeasible);
computeSearchDirection(betaCorrector, mu, true);
-
+ timers["corrector"].stop();
+ timers["step length"].resume();
// Step length to preserve positive definiteness
primalStepLength = stepLength(XCholesky, dX, StepMatrixWorkspace,
eigenvaluesWorkspace, QRWorkspace, parameters);
dualStepLength = stepLength(YCholesky, dY, StepMatrixWorkspace,
eigenvaluesWorkspace, QRWorkspace, parameters);
-
+ timers["step length"].stop();
if (isPrimalFeasible && isDualFeasible) {
primalStepLength = min(primalStepLength, dualStepLength);
dualStepLength = primalStepLength;
}
- printSolverInfo(iteration, mu, status, primalStepLength, dualStepLength, betaCorrector);
+ printSolverInfo(iteration, mu, status, primalStepLength, dualStepLength, betaCorrector, sdp.dualObjective.size(), Q.rows);
// Update current point
scaleVector(dx, primalStepLength);
diff --git a/src/SDPSolver.h b/src/SDPSolver.h
index 67c1b7e..b81af6c 100644
--- a/src/SDPSolver.h
+++ b/src/SDPSolver.h
@@ -34,6 +34,7 @@ public:
Real feasibleCenteringParameter;
Real infeasibleCenteringParameter;
Real stepLengthReduction;
+ Real choleskyStabilizeThreshold;
Real maxComplementarity;
void resetPrecision() {
@@ -45,6 +46,7 @@ public:
setPrecision(feasibleCenteringParameter, precision);
setPrecision(infeasibleCenteringParameter, precision);
setPrecision(stepLengthReduction, precision);
+ setPrecision(choleskyStabilizeThreshold, precision);
setPrecision(maxComplementarity, precision);
}
@@ -130,7 +132,8 @@ public:
void initialize(const SDPSolverParameters ¶meters);
SDPSolverTerminateReason run(const SDPSolverParameters ¶meters, const path checkpointFile);
void initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
- const BlockDiagonalMatrix &BilinearPairingsY);
+ const BlockDiagonalMatrix &BilinearPairingsY,
+ const Real &choleskyStabilizeThreshold);
void solveSchurComplementEquation(Vector &dx, Vector &dz);
void computeSearchDirection(const Real &beta, const Real &mu, const bool correctorPhase);
Vector freeVariableSolution();
@@ -145,6 +148,8 @@ void printSolverInfo(int iteration,
SDPSolverStatus status,
Real primalStepLength,
Real dualStepLength,
- Real betaCorrector);
+ Real betaCorrector,
+ int dualObjectiveSize,
+ int Qrows);
#endif // SDP_BOOTSTRAP_SDPSOLVER_H_
diff --git a/src/SDPSolverIO.cpp b/src/SDPSolverIO.cpp
index 045d78f..0a0fee0 100644
--- a/src/SDPSolverIO.cpp
+++ b/src/SDPSolverIO.cpp
@@ -17,8 +17,8 @@ using boost::posix_time::microseconds;
using std::cout;
void printSolverHeader() {
- cout << "\n time mu P-obj D-obj gap P-err D-err P-step D-step beta\n";
- cout << "---------------------------------------------------------------------------------------------------------\n";
+ cout << "\n time mu P-obj D-obj gap P-err D-err P-step D-step beta dim/stabilized\n";
+ cout << "-------------------------------------------------------------------------------------------------------------------------\n";
}
void printSolverInfo(int iteration,
@@ -26,13 +26,15 @@ void printSolverInfo(int iteration,
SDPSolverStatus status,
Real primalStepLength,
Real dualStepLength,
- Real betaCorrector) {
+ Real betaCorrector,
+ int dualObjectiveSize,
+ int Qrows) {
Real time = Real(timers["Solver runtime"].elapsed().wall)/1000000000;
time_duration td(microseconds(timers["Solver runtime"].elapsed().wall)/1000);
std::stringstream ss;
ss << td;
gmp_fprintf(stdout,
- "%3d %s %-8.1Fe %-+11.2Fe %-+11.2Fe %-9.2Fe %-+10.2Fe %-+10.2Fe %-8.3Fg %-8.3Fg %-4.2Fg",
+ "%3d %s %-8.1Fe %-+11.2Fe %-+11.2Fe %-9.2Fe %-+10.2Fe %-+10.2Fe %-8.3Fg %-8.3Fg %-4.2Fg %d/%d",
iteration,
ss.str().substr(0,8).c_str(),
mu.get_mpf_t(),
@@ -43,7 +45,9 @@ void printSolverInfo(int iteration,
status.dualError.get_mpf_t(),
primalStepLength.get_mpf_t(),
dualStepLength.get_mpf_t(),
- betaCorrector.get_mpf_t());
+ betaCorrector.get_mpf_t(),
+ dualObjectiveSize,
+ Qrows);
cout << endl;
}
@@ -65,6 +69,7 @@ ostream& operator<<(ostream& os, const SDPSolverParameters& p) {
os << "feasibleCenteringParameter = " << p.feasibleCenteringParameter << endl;
os << "infeasibleCenteringParameter = " << p.infeasibleCenteringParameter << endl;
os << "stepLengthReduction = " << p.stepLengthReduction << endl;
+ os << "choleskyStabilizeThreshold = " << p.choleskyStabilizeThreshold << endl;
os << "maxComplementarity = " << p.maxComplementarity << endl;
return os;
}
diff --git a/src/main.cpp b/src/main.cpp
index 4ec7336..4275bf0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -59,7 +59,7 @@ int solveSDP(const path &sdpFile,
SDPSolverTerminateReason reason = solver.run(parameters, checkpointFile);
timers["Solver runtime"].stop();
- cout << "-----" << setfill('-') << setw(100) << std::left << reason << endl << endl;
+ cout << "-----" << setfill('-') << setw(116) << std::left << reason << endl << endl;
cout << solver.status << endl;
if (!parameters.noFinalCheckpoint)
@@ -164,6 +164,13 @@ int main(int argc, char** argv) {
po::value<Real>(¶meters.stepLengthReduction)->default_value(Real("0.7")),
"Shrink each newton step by this factor (smaller means slower, more "
"stable convergence). Corresponds to SDPA's gammaStar.")
+ ("choleskyStabilizeThreshold",
+ po::value<Real>(¶meters.choleskyStabilizeThreshold)->default_value(Real("1e-40")),
+ "Adds stabilizing terms to the cholesky decomposition of the schur complement "
+ "matrix for diagonal entries which are smaller than this threshold times the "
+ "geometric mean of other diagonal entries. Somewhat higher choleskyStabilizeThreshold "
+ "can improve numerical stability but if the threshold is large enough that a high "
+ "proportion of eigenvalues are being stabilized, the computation will slow substantially.")
("maxComplementarity",
po::value<Real>(¶meters.maxComplementarity)->default_value(Real("1e100")),
"Terminate if the complementarity mu = Tr(X Y)/dim(X) exceeds this value.")
diff --git a/src/mpack/Rgetrf.cpp b/src/mpack/Rgetrf.cpp
index 6cf6257..958252d 100644
--- a/src/mpack/Rgetrf.cpp
+++ b/src/mpack/Rgetrf.cpp
@@ -68,6 +68,8 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <mblas.h>
#include <mlapack.h>
+#include "../Timers.h"
+
void Rgetrf(INTEGER m, INTEGER n, REAL * A, INTEGER lda, INTEGER * ipiv, INTEGER * info)
{
INTEGER i, j, jb, nb, iinfo;
@@ -101,7 +103,9 @@ void Rgetrf(INTEGER m, INTEGER n, REAL * A, INTEGER lda, INTEGER * ipiv, INTEGER
jb = min(min(m, n) - j + 1, nb);
//Factor diagonal and subdiagonal blocks and test for exact
//singularity.
+ timers["Rgetf2"].resume();
Rgetf2(m - j + 1, jb, &A[(j - 1) + (j - 1) * lda], lda, &ipiv[j - 1], &iinfo);
+ timers["Rgetf2"].stop();
//Adjust INFO and the pivot indices.
if (*info == 0 && iinfo > 0) {
*info = iinfo + j - 1;
@@ -110,16 +114,24 @@ void Rgetrf(INTEGER m, INTEGER n, REAL * A, INTEGER lda, INTEGER * ipiv, INTEGER
ipiv[i - 1] = j - 1 + ipiv[i - 1];
}
//Apply interchanges to columns 1:J-one
+ timers["Rlaswp"].resume();
Rlaswp(j - 1, A, lda, j, j + jb - 1, ipiv, 1);
+ timers["Rlaswp"].stop();
if (j + jb <= n) {
//Apply interchanges to columns J+JB:N.
+ timers["Rlaswp2"].resume();
Rlaswp(n - j - jb + 1, &A[0 + (j + jb - 1) * lda], lda, j, j + jb - 1, ipiv, 1);
+ timers["Rlaswp2"].stop();
//Compute block row of U.
+ timers["Rtrsm"].resume();
Rtrsm("Left", "Lower", "No transpose", "Unit", jb, n - j - jb + 1, One, &A[(j - 1) + (j - 1) * lda], lda, &A[(j - 1) + (j + jb - 1) * lda], lda);
+ timers["Rtrsm"].stop();
if (j + jb <= m) {
//Update trailing submatrix.
+ timers["RgemmParallel"].resume();
RgemmParallel("No transpose", "No transpose", m - j - jb + 1,
n - j - jb + 1, jb, -One, &A[(j + jb - 1) + (j - 1) * lda], lda, &A[(j - 1) + (j + jb - 1) * lda], lda, One, &A[(j + jb - 1) + (j + jb - 1) * lda], lda);
+ timers["RgemmParallel"].stop();
}
}
}
diff --git a/src/mpack/Rpotf2Stabilized.cpp b/src/mpack/Rpotf2Stabilized.cpp
index 065f0ed..e22f6fa 100644
--- a/src/mpack/Rpotf2Stabilized.cpp
+++ b/src/mpack/Rpotf2Stabilized.cpp
@@ -71,26 +71,24 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <vector>
#include <assert.h>
-const double CHOLESKY_STABILIZE_THRESHOLD = 1e-60;
-
-REAL lambdaGeometricMean(double totalLogLambda, INTEGER numLambdas) {
+REAL lambdaGeometricMean(const double totalLogLambda, const INTEGER numLambdas) {
assert(numLambdas != 0);
return REAL(exp(totalLogLambda/numLambdas));
}
-bool isSmallDiagonal(const REAL &ajj, double totalLogLambda, INTEGER numLambdas) {
+bool isSmallDiagonal(const REAL &ajj, const double totalLogLambda, const INTEGER numLambdas, const double stabilizeThreshold) {
double d = cast2double(ajj);
- return (d <= 0) || (numLambdas != 0 && log(d)/2 < totalLogLambda/numLambdas + log(CHOLESKY_STABILIZE_THRESHOLD));
+ return (d <= 0) || (numLambdas != 0 && log(d)/2 < totalLogLambda/numLambdas + log(stabilizeThreshold));
}
-void correctDiagonal(REAL &ajj, INTEGER jIndex, double totalLogLambda, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas) {
+void correctDiagonal(REAL &ajj, const INTEGER jIndex, const double totalLogLambda, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas) {
REAL lambda = lambdaGeometricMean(totalLogLambda, jIndex);
ajj += lambda * lambda;
stabilizeIndices.push_back(jIndex);
stabilizeLambdas.push_back(lambda);
}
-void Rpotf2Stabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, INTEGER indexStart, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas, double &totalLogLambda)
+void Rpotf2Stabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, const INTEGER indexStart, const double stabilizeThreshold, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas, double &totalLogLambda)
{
INTEGER j, upper, success = 1;
REAL ajj;
@@ -119,7 +117,7 @@ void Rpotf2Stabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
//Compute U(J,J) and test for non-positive-definiteness.
ajj = A[j + j * lda] - Rdot(j, &A[j * lda], 1, &A[j * lda], 1);
- if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j))
+ if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j, stabilizeThreshold))
correctDiagonal(ajj, indexStart + j, totalLogLambda, stabilizeIndices, stabilizeLambdas);
ajj = sqrt(ajj);
@@ -138,7 +136,7 @@ void Rpotf2Stabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
//Compute L(J,J) and test for non-positive-definiteness.
ajj = A[j + j * lda] - Rdot(j, &A[j], lda, &A[j], lda);
- if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j))
+ if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j, stabilizeThreshold))
correctDiagonal(ajj, indexStart + j, totalLogLambda, stabilizeIndices, stabilizeLambdas);
ajj = sqrt(ajj);
diff --git a/src/mpack/RpotrfStabilized.cpp b/src/mpack/RpotrfStabilized.cpp
index e253aa2..c3cdd9c 100644
--- a/src/mpack/RpotrfStabilized.cpp
+++ b/src/mpack/RpotrfStabilized.cpp
@@ -70,7 +70,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <vector>
-void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas)
+void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas, const double stabilizeThreshold)
{
INTEGER upper;
INTEGER j, jb, nb;
@@ -101,7 +101,7 @@ void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
nb = iMlaenv(1, "Rpotrf", uplo, n, -1, -1, -1);
if (nb <= 1 || nb >= n) {
//Use unblocked code.
- Rpotf2Stabilized(uplo, n, A, lda, info, 0, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+ Rpotf2Stabilized(uplo, n, A, lda, info, 0, stabilizeThreshold, stabilizeIndices, stabilizeLambdas, totalLogLambda);
} else {
//Use blocked code.
if (upper) {
@@ -111,7 +111,7 @@ void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
//for non-positive-definiteness.
jb = min(nb, n - j + 1);
Rsyrk("Upper", "Transpose", jb, j - 1, -One, &A[0 + (j - 1) * lda], lda, One, &A[(j - 1) + (j - 1) * lda], lda);
- Rpotf2Stabilized("Upper", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+ Rpotf2Stabilized("Upper", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeThreshold, stabilizeIndices, stabilizeLambdas, totalLogLambda);
if (*info != 0) {
goto L30;
}
@@ -129,7 +129,7 @@ void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
//for non-positive-definiteness.
jb = min(nb, n - j + 1);
Rsyrk("Lower", "No transpose", jb, j - 1, -One, &A[(j - 1) + 0 * lda], lda, One, &A[(j - 1) + (j - 1) * lda], lda);
- Rpotf2Stabilized("Lower", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+ Rpotf2Stabilized("Lower", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeThreshold, stabilizeIndices, stabilizeLambdas, totalLogLambda);
if (*info != 0) {
goto L30;
}
diff --git a/src/mpack/mlapack_gmp.h b/src/mpack/mlapack_gmp.h
index fdab543..a7bc622 100644
--- a/src/mpack/mlapack_gmp.h
+++ b/src/mpack/mlapack_gmp.h
@@ -44,7 +44,7 @@ void
Rsyev(const char *jobz, const char *uplo, mpackint n, mpf_class * A,
mpackint lda, mpf_class * w, mpf_class * work, mpackint lwork, mpackint *info);
void Rpotrf(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint *info);
-void RpotrfStabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas);
+void RpotrfStabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas, const double stabilizeThreshold);
mpackint iMlaenv_gmp(mpackint ispec, const char *name, const char *opts, mpackint n1, mpackint n2,
mpackint n3, mpackint n4);
mpf_class Rlamch_gmp(const char *cmach);
@@ -82,7 +82,7 @@ void Rorg2r(mpackint m, mpackint n, mpackint k, mpf_class * A, mpackint lda, mpf
void Rlarf(const char *side, mpackint m, mpackint n, mpf_class * v, mpackint incv,
mpf_class tau, mpf_class * C, mpackint ldc, mpf_class * work);
void Rpotf2(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint *info);
-void Rpotf2Stabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, mpackint indexStart, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas, double &totalLogLambda);
+void Rpotf2Stabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, const mpackint indexStart, const double stabilizeThreshold, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas, double &totalLogLambda);
void Rlaset(const char *uplo, mpackint m, mpackint n, mpf_class alpha, mpf_class beta,
mpf_class * A, mpackint lda);
void Rlaev2(mpf_class a, mpf_class b, mpf_class c, mpf_class * rt1,
diff --git a/src/tests.cpp b/src/tests.cpp
index a0de85d..1643971 100644
--- a/src/tests.cpp
+++ b/src/tests.cpp
@@ -30,7 +30,8 @@ void testCholeskyStabilize() {
vector<Integer> stabilizeIndices;
vector<Real> stabilizeLambdas;
- choleskyDecompositionStabilized(A, L, stabilizeIndices, stabilizeLambdas);
+ double stabilizeThreshold = 1e-10;
+ choleskyDecompositionStabilized(A, L, stabilizeIndices, stabilizeLambdas, stabilizeThreshold);
cout << "A = " << A << ";\n";
cout << "L = " << L << ";\n";
cout << "stabilizeIndices = " << stabilizeIndices << ";\n";
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/sdpb.git
More information about the debian-science-commits
mailing list