[sdpb] 86/233: More fixes to cholesky stabilization
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:22 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 323148fdac6d9552fd935bf2f4baff665a16b919
Author: David Simmons-Duffin <dsd at athena.sns.ias.edu>
Date: Sun Oct 5 21:35:48 2014 -0400
More fixes to cholesky stabilization
---
src/BlockDiagonalMatrix.cpp | 9 +++++++++
src/BlockDiagonalMatrix.h | 4 ++++
src/Matrix.cpp | 8 ++++++++
src/SDPSolver.cpp | 20 ++++----------------
src/SDPSolver.h | 4 ++--
src/main.cpp | 3 +--
src/mpack/RpotrfStabilized.cpp | 2 ++
7 files changed, 30 insertions(+), 20 deletions(-)
diff --git a/src/BlockDiagonalMatrix.cpp b/src/BlockDiagonalMatrix.cpp
index 27846cd..8014680 100644
--- a/src/BlockDiagonalMatrix.cpp
+++ b/src/BlockDiagonalMatrix.cpp
@@ -97,6 +97,15 @@ void choleskyDecomposition(BlockDiagonalMatrix &A,
choleskyDecomposition(A.blocks[b], L.blocks[b]);
}
+void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
+ BlockDiagonalMatrix &L,
+ vector<vector<Integer> > &schurStabilizeIndices,
+ vector<vector<Real> > &schurStabilizeLambdas) {
+ #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]);
+}
+
void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
BlockDiagonalMatrix &X) {
#pragma omp parallel for schedule(dynamic)
diff --git a/src/BlockDiagonalMatrix.h b/src/BlockDiagonalMatrix.h
index 3bc99c3..d26e7c5 100644
--- a/src/BlockDiagonalMatrix.h
+++ b/src/BlockDiagonalMatrix.h
@@ -101,6 +101,10 @@ void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A, BlockDiagonalMatrix &B,
void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L);
Real minEigenvalue(BlockDiagonalMatrix &A, vector<Vector> &workspace, vector<Vector> &eigenvalues);
void choleskyDecomposition(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L);
+void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
+ BlockDiagonalMatrix &L,
+ vector<vector<Integer> > &schurStabilizeIndices,
+ vector<vector<Real> > &schurStabilizeLambdas);
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 34d1179..4cb7e8c 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -238,6 +238,14 @@ void choleskyDecomposition(Matrix &A, Matrix &L) {
L.elements[i + j*dim] = 0;
}
+// L (lower triangular) such that A = L L^T - \sum_i stabilizeLambdas_i^2 u_j u_j^T
+// where j = stabilizeIndices_i
+// Inputs:
+// - A : dim x dim symmetric matrix
+// - L : dim x dim lower-triangular matrix
+// - 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) {
int dim = A.rows;
assert(A.cols == dim);
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 7a40992..56286cd 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -37,7 +37,6 @@ SDPSolver::SDPSolver(const SDP &sdp):
basicKernelCoords(Q.rows),
schurStabilizeIndices(SchurBlocks.blocks.size()),
schurStabilizeLambdas(SchurBlocks.blocks.size()),
- schurStabilizeVectors(SchurBlocks.blocks.size()),
StepMatrixWorkspace(X)
{
// initialize bilinearPairingsWorkspace, eigenvaluesWorkspace, QRWorkspace
@@ -82,9 +81,6 @@ SDPSolver::SDPSolver(const SDP &sdp):
FreeVarMatrixBasicPivots,
&dualObjectiveReduced[0], 1,
dualObjectiveReduced.size());
-
- for (unsigned int b = 0; b < SchurBlocks.blocks.size(); b++)
- schurStabilizeVectors[b].resize(SchurBlocks.blocks[b].rows);
}
Vector SDPSolver::freeVariableSolution() {
@@ -483,14 +479,8 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
const BlockDiagonalMatrix &BilinearPairingsY) {
computeSchurBlocks(sdp, BilinearPairingsXInv, BilinearPairingsY, SchurBlocks);
- choleskyDecomposition(SchurBlocks, SchurBlocksCholesky);
+ choleskyDecompositionStabilized(SchurBlocks, SchurBlocksCholesky, schurStabilizeIndices, schurStabilizeLambdas);
- for (unsigned int b = 0; b < SchurBlocks.blocks.size(); b++)
- stabilizeCholesky(SchurBlocksCholesky.blocks[b],
- schurStabilizeVectors[b],
- schurStabilizeIndices[b],
- schurStabilizeLambdas[b]);
-
// SchurUpdateLowRank = {{- 1, 0}, {E, G}}
SchurUpdateLowRank.setCols(sdp.FreeVarMatrix.cols);
SchurUpdateLowRank.copyFrom(sdp.FreeVarMatrix);
@@ -498,9 +488,8 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
for (unsigned int i = 0; i < schurStabilizeIndices[b].size(); i++) {
int fullIndex = SchurBlocks.blockStartIndices[b] + schurStabilizeIndices[b][i];
SchurUpdateLowRank.addColumn();
- SchurUpdateLowRank.elt(fullIndex, SchurUpdateLowRank.cols - 1) = schurStabilizeLambdas[b];
+ SchurUpdateLowRank.elt(fullIndex, SchurUpdateLowRank.cols - 1) = schurStabilizeLambdas[b][i];
}
- schurStabilizeIndices[b].resize(0);
}
// SchurUpdateLowRank = SchurBlocksCholesky^{-1} {{- 1, 0}, {E, G}}
@@ -518,7 +507,6 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
}
void SDPSolver::solveSchurComplementEquation(Vector &dx) {
-
// dx = SchurBlocksCholesky^{-1} dx
blockMatrixLowerTriangularSolve(SchurBlocksCholesky, dx);
@@ -671,7 +659,7 @@ SDPSolverTerminateReason SDPSolver::run(const SDPSolverParameters ¶meters,
}
timers["Run solver"].stop();
- saveCheckpoint(checkpointFile);
- timers["Save checkpoint"].start();
+ //saveCheckpoint(checkpointFile);
+ //timers["Save checkpoint"].start();
return finished;
}
diff --git a/src/SDPSolver.h b/src/SDPSolver.h
index 9df8edb..a74893c 100644
--- a/src/SDPSolver.h
+++ b/src/SDPSolver.h
@@ -114,8 +114,8 @@ public:
vector<Integer> Qpivots;
Vector basicKernelCoords;
vector<vector<int> > schurStabilizeIndices;
- vector<Real> schurStabilizeLambdas;
- vector<Vector> schurStabilizeVectors;
+ vector<vector<Real> > schurStabilizeLambdas;
+ //vector<Vector> schurStabilizeVectors;
// additional workspace variables
BlockDiagonalMatrix StepMatrixWorkspace;
diff --git a/src/main.cpp b/src/main.cpp
index 28faac9..6271b18 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -184,6 +184,5 @@ int main(int argc, char** argv) {
return 1;
}
- //return solveSDP(sdpFile, outFile, checkpointFile, parameters);
- testCholeskyStabilize();
+ return solveSDP(sdpFile, outFile, checkpointFile, parameters);
}
diff --git a/src/mpack/RpotrfStabilized.cpp b/src/mpack/RpotrfStabilized.cpp
index cbf787b..e253aa2 100644
--- a/src/mpack/RpotrfStabilized.cpp
+++ b/src/mpack/RpotrfStabilized.cpp
@@ -76,6 +76,8 @@ void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGE
INTEGER j, jb, nb;
REAL One = 1.0;
double totalLogLambda = 0;
+ stabilizeIndices.resize(0);
+ stabilizeLambdas.resize(0);
//Test the input parameters.
*info = 0;
--
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