[sdpb] 127/233: Commented Matrix, BlockDiagonalMatrix, and main; some minor fixes and cleanup
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:28 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 44201eb57e0d1a377cf006fc4402c7ce44dea006
Author: David Simmons-Duffin <dsd at minerva.sns.ias.edu>
Date: Sat Jan 10 03:05:07 2015 -0500
Commented Matrix, BlockDiagonalMatrix, and main; some minor fixes and cleanup
---
src/BlockDiagonalMatrix.cpp | 17 ++++++-
src/BlockDiagonalMatrix.h | 106 +++++++++++++++++++++++++++++++++++---------
src/Matrix.cpp | 77 ++++++++++++++------------------
src/Matrix.h | 99 +++++++++++++++++++++++++++++++++--------
src/SDPSolver.cpp | 17 ++++---
src/main.cpp | 11 ++++-
6 files changed, 232 insertions(+), 95 deletions(-)
diff --git a/src/BlockDiagonalMatrix.cpp b/src/BlockDiagonalMatrix.cpp
index 7ee9dbb..4897aea 100644
--- a/src/BlockDiagonalMatrix.cpp
+++ b/src/BlockDiagonalMatrix.cpp
@@ -21,12 +21,16 @@ ostream& operator<<(ostream& os, const BlockDiagonalMatrix& A) {
return os;
}
+// Tr(A B), where A and B are symmetric
Real frobeniusProductSymmetric(const BlockDiagonalMatrix &A,
const BlockDiagonalMatrix &B) {
Real result = 0;
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < A.blocks.size(); b++) {
Real f = frobeniusProductSymmetric(A.blocks[b], B.blocks[b]);
+ // this pragma means that other threads should be stopped while
+ // this operation is performed (this prevents result from being
+ // modified by multiple threads simultaneously)
#pragma omp critical
{
result += f;
@@ -55,6 +59,7 @@ Real frobeniusProductOfSums(const BlockDiagonalMatrix &X,
return result;
}
+// C := alpha*A*B + beta*C
void blockDiagonalMatrixScaleMultiplyAdd(Real alpha,
BlockDiagonalMatrix &A,
BlockDiagonalMatrix &B,
@@ -65,12 +70,14 @@ void blockDiagonalMatrixScaleMultiplyAdd(Real alpha,
matrixScaleMultiplyAdd(alpha, A.blocks[b], B.blocks[b], beta, C.blocks[b]);
}
+// C := A*B
void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &B,
BlockDiagonalMatrix &C) {
blockDiagonalMatrixScaleMultiplyAdd(1, A, B, 0, C);
}
+// A := L^{-1} A L^{-T}
void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L) {
#pragma omp parallel for schedule(dynamic)
@@ -91,12 +98,13 @@ Real minEigenvalue(BlockDiagonalMatrix &A,
assert(A.blocks.size() == workspace.size());
// TODO(davidsd): get rid of this hack
- Real lambdaMin = 1e100;
+ Real lambdaMin = 1e100; // we really want lambdaMin = infinity here
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < A.blocks.size(); b++) {
Real minBlockLambda = minEigenvalue(A.blocks[b],
workspace[b],
eigenvalues[b]);
+ // ensure only one thread modifies lambdaMin at a time
#pragma omp critical
{
lambdaMin = min(lambdaMin, minBlockLambda);
@@ -106,6 +114,7 @@ Real minEigenvalue(BlockDiagonalMatrix &A,
return lambdaMin;
}
+// Compute L (lower triangular) such that A = L L^T
void choleskyDecomposition(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L) {
#pragma omp parallel for schedule(dynamic)
@@ -113,6 +122,8 @@ void choleskyDecomposition(BlockDiagonalMatrix &A,
choleskyDecomposition(A.blocks[b], L.blocks[b]);
}
+// Compute L (lower triangular) such that A + U U^T = L L^T, where U
+// is a low-rank update matrix.
void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L,
vector<vector<Integer> > &stabilizeIndices,
@@ -126,6 +137,7 @@ void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
stabilizeThreshold);
}
+// X := ACholesky^{-T} ACholesky^{-1} X = A^{-1} X
void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
BlockDiagonalMatrix &X) {
#pragma omp parallel for schedule(dynamic)
@@ -133,6 +145,7 @@ void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
matrixSolveWithCholesky(ACholesky.blocks[b], X.blocks[b]);
}
+// B := L^{-1} B, where L is lower-triangular
void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
Matrix &B) {
#pragma omp parallel for schedule(dynamic)
@@ -141,6 +154,7 @@ void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
B.cols, B.rows);
}
+// v := L^{-1} v, where L is lower-triangular
void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
Vector &v) {
#pragma omp parallel for schedule(dynamic)
@@ -149,6 +163,7 @@ void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
1, v.size());
}
+// v := L^{-T} v, where L is lower-triangular
void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L,
Vector &v) {
#pragma omp parallel for schedule(dynamic)
diff --git a/src/BlockDiagonalMatrix.h b/src/BlockDiagonalMatrix.h
index f11eb9b..23f877d 100644
--- a/src/BlockDiagonalMatrix.h
+++ b/src/BlockDiagonalMatrix.h
@@ -19,12 +19,26 @@
using std::vector;
using std::ostream;
+// A block-diagonal square matrix
+//
+// M = Diagonal(M_0, M_1, ..., M_{bMax-1})
+//
+// where each block M_b is a square-matrix (of possibly different
+// sizes).
class BlockDiagonalMatrix {
public:
+ // The total number of rows (or columns) in M
int dim;
+
+ // The blocks M_b for 0 <= b < bMax
vector<Matrix> blocks;
+
+ // The rows (or columns) of M corresponding to the top-left entry of
+ // each block M_b
vector<int> blockStartIndices;
+ // Construct a BlockDiagonalMatrix from a vector of dimensions {s_0,
+ // ..., s_{bMax-1}} for each block.
explicit BlockDiagonalMatrix(const vector<int> &blockSizes):
dim(0) {
for (unsigned int i = 0; i < blockSizes.size(); i++) {
@@ -34,47 +48,55 @@ public:
}
}
+ // M = 0
void setZero() {
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b].setZero();
}
+ // Add a constant c to each diagonal element
void addDiagonal(const Real &c) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b].addDiagonal(c);
}
+ // M = M + A
void operator+=(const BlockDiagonalMatrix &A) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b] += A.blocks[b];
}
+ // M = M - A
void operator-=(const BlockDiagonalMatrix &A) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b] -= A.blocks[b];
}
+ // M = c*M, where c is a constant
void operator*=(const Real &c) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b] *= c;
}
+ // M = A
void copyFrom(const BlockDiagonalMatrix &A) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b].copyFrom(A.blocks[b]);
}
+ // Symmetrize M in place
void symmetrize() {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < blocks.size(); b++)
blocks[b].symmetrize();
}
+ // The maximal absolute value of the elements of M
Real maxAbs() const {
Real max = 0;
for (vector<Matrix>::const_iterator b = blocks.begin();
@@ -90,6 +112,7 @@ public:
friend ostream& operator<<(ostream& os, const BlockDiagonalMatrix& A);
};
+// Tr(A B), where A and B are symmetric
Real frobeniusProductSymmetric(const BlockDiagonalMatrix &A,
const BlockDiagonalMatrix &B);
@@ -101,42 +124,81 @@ Real frobeniusProductOfSums(const BlockDiagonalMatrix &X,
const BlockDiagonalMatrix &Y,
const BlockDiagonalMatrix &dY);
+// C := alpha*A*B + beta*C
void blockDiagonalMatrixScaleMultiplyAdd(Real alpha,
- BlockDiagonalMatrix &A,
- BlockDiagonalMatrix &B,
+ BlockDiagonalMatrix &A, // constant
+ BlockDiagonalMatrix &B, // constant
Real beta,
- BlockDiagonalMatrix &C);
-
-void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A,
- BlockDiagonalMatrix &B,
- BlockDiagonalMatrix &C);
-
-void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A,
- BlockDiagonalMatrix &L);
-
+ BlockDiagonalMatrix &C); // overwritten
+
+// C := A B
+void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A, // constant
+ BlockDiagonalMatrix &B, // constant
+ BlockDiagonalMatrix &C); // overwritten
+
+// A := L^{-1} A L^{-T}
+void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A, // overwritten
+ BlockDiagonalMatrix &L); // constant
+
+// Minimum eigenvalue of A, via the QR method
+// Inputs:
+// - A : BlockDiagonalMatrix with blocks of size n_b x n_b (will be
+// overwritten)
+// - eigenvalues : vector of Vectors of length n_b (0 <= b < bMax)
+// - workspace : vector of Vectors of lenfth 3*n_b-1 (0 <= b < bMax)
+// (temporary workspace)
+//
Real minEigenvalue(BlockDiagonalMatrix &A,
vector<Vector> &workspace,
vector<Vector> &eigenvalues);
+// Compute L (lower triangular) such that A = L L^T
+// Inputs:
+// - A : dim x dim symmetric matrix (constant)
+// - L : dim x dim lower-triangular matrix (overwritten)
+//
void choleskyDecomposition(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L);
+// Compute L (lower triangular) such that A + U U^T = L L^T, where U
+// is a low-rank update matrix. (See Matrix.h for details).
+//
+// Input:
+// - A
+// - stabilizeThreshold: parameter for deciding which directions to
+// stabilize
+// Output:
+// - L (overwritten)
+// - stabilizeIndices: a list of directions which have been
+// stabilized, for each block of A (overwritten)
+// - stabilizeLambdas: a list of corresponding Lambdas, for each block
+// of A (overwritten)
+//
void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L,
vector<vector<Integer> > &stabilizeIndices,
vector<vector<Real> > &stabilizeLambdas,
const double stabilizeThreshold);
-void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
- BlockDiagonalMatrix &X);
-
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
- Matrix &B);
-
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
- Vector &v);
-
-void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L,
- Vector &v);
+// X := ACholesky^{-T} ACholesky^{-1} X = A^{-1} X
+// Input:
+// - ACholesky, a lower triangular BlockDiagonalMatrix such that A =
+// ACholesky ACholesky^T (constant)
+// Output:
+// - X (overwritten)
+void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky, // constant
+ BlockDiagonalMatrix &X); // overwritten
+
+// B := L^{-1} B, where L is lower-triangular
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, // constant
+ Matrix &B); // overwritten
+
+// v := L^{-1} v, where L is lower-triangular
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, // constant
+ Vector &v); // overwritten
+
+// v := L^{-T} v, where L is lower-triangular
+void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, // constant
+ Vector &v); // overwritten
#endif // SDPB_BLOCKDIAGONALMATRIX_H_
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index bbdd510..62cc036 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -9,6 +9,20 @@
#include <vector>
#include "Matrix.h"
+// Most of the routines below are helpfully-named wrappers for
+// functions in MBLAS or MLAPACK. See Matrix.h for a more detailed
+// description of the various input/output parameters.
+//
+// For a list of MBLAS routines with documentation, see
+// http://mplapack.sourceforge.net/mblas_routines.html
+//
+// For a list of MLAPACK routines with documentation, see
+// http://mplapack.sourceforge.net/mlapack_routines.html
+//
+// We have chosen not to parallelize the operations here that are used
+// within BlockDiagonalMatrices, since there parallelism can be
+// achieved by parallelizing loops over blocks.
+
ostream& operator<<(ostream& os, const Matrix& a) {
os << "{";
for (int r = 0; r < a.rows; r++) {
@@ -26,18 +40,7 @@ ostream& operator<<(ostream& os, const Matrix& a) {
return os;
}
-// B := A^T
-void transpose(const Matrix &A, Matrix &B) {
- assert(A.cols == B.rows);
- assert(A.rows == B.cols);
-
- for (int n = 0; n < A.cols; n++)
- for (int m = 0; m < A.rows; m++)
- B.elt(n, m) = A.elt(m, n);
-}
-
// C := alpha*A*B + beta*C
-//
void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B,
Real beta, Matrix &C) {
assert(A.cols == B.rows);
@@ -52,7 +55,6 @@ void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B,
}
// C := A*B
-//
void matrixMultiply(Matrix &A, Matrix &B, Matrix &C) {
matrixScaleMultiplyAdd(1, A, B, 0, C);
}
@@ -62,6 +64,11 @@ void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol) {
assert(bRow + A.cols <= B.rows);
assert(bCol + A.cols <= B.cols);
+ // This operation is not used within a BlockDiagonalMatrix, so it is
+ // worthwhile to parallelize. In fact, this function, used in
+ // computing TopLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal (see
+ // SDPSolver.cpp) is one of the main performance bottlenecks in the
+ // solver.
#pragma omp parallel for schedule(dynamic)
for (int c = 0; c < A.cols; c++) {
for (int r = 0; r <= c; r++) {
@@ -82,17 +89,18 @@ void lowerTriangularInverseCongruence(Matrix &A, Matrix &L) {
assert(L.rows == dim);
assert(L.cols == dim);
+ // A := A L^{-T}
Rtrsm("Right", "Lower", "Transpose", "NonUnitDiagonal", dim, dim, 1,
&L.elements[0], dim,
&A.elements[0], dim);
+ // A := L^{-1} A
Rtrsm("Left", "Lower", "NoTranspose", "NonUnitDiagonal", dim, dim, 1,
&L.elements[0], dim,
&A.elements[0], dim);
}
// y := alpha A x + beta y
-//
void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x,
Real beta, Vector &y) {
assert(A.cols <= static_cast<int>(x.size()));
@@ -108,7 +116,6 @@ void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x,
}
// y := alpha A^T x + beta y
-//
void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x,
Real beta, Vector &y) {
assert(A.cols <= static_cast<int>(y.size()));
@@ -123,6 +130,7 @@ void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x,
}
}
+// Tr(A B), where A and B are symmetric
Real frobeniusProductSymmetric(const Matrix &A, const Matrix &B) {
assert(A.rows == B.rows);
assert(A.cols == B.cols);
@@ -158,12 +166,7 @@ Real frobeniusProductOfSums(const Matrix &X, const Matrix &dX,
return result;
}
-// Eigenvalues of A, via the QR method
-// Inputs:
-// A : n x n Matrix (will be overwritten)
-// eigenvalues : Vector of length n
-// workspace : Vector of lenfth 3*n-1 (temporary workspace)
-//
+// Compute the eigenvalues of A, via the QR method
void matrixEigenvalues(Matrix &A, Vector &workspace, Vector &eigenvalues) {
assert(A.rows == A.cols);
assert(static_cast<int>(eigenvalues.size()) == A.rows);
@@ -177,16 +180,13 @@ void matrixEigenvalues(Matrix &A, Vector &workspace, Vector &eigenvalues) {
}
// Minimum eigenvalue of A, via the QR method
-// Inputs:
-// A : n x n Matrix (will be overwritten)
-// eigenvalues : Vector of length n
-// workspace : Vector of lenfth 3*n-1 (temporary workspace)
-//
Real minEigenvalue(Matrix &A, Vector &workspace, Vector &eigenvalues) {
matrixEigenvalues(A, workspace, eigenvalues);
return eigenvalues[0];
}
+// Compute an in-place LU decomposition of A, with pivots, suitable
+// for use with 'solveWithLUDecomposition'
void LUDecomposition(Matrix &A, vector<Integer> &pivots) {
int dim = A.rows;
assert(A.cols == dim);
@@ -196,6 +196,7 @@ void LUDecomposition(Matrix &A, vector<Integer> &pivots) {
assert(info == 0);
}
+// b := A^{-1} b, where LU and pivots encode the LU decomposition of A
void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b) {
Integer info;
Rgetrs("NoTranspose", LU.rows, 1, &LU.elements[0], LU.rows,
@@ -204,10 +205,6 @@ void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b) {
}
// L (lower triangular) such that A = L L^T
-// Inputs:
-// - A : dim x dim symmetric matrix
-// - L : dim x dim lower-triangular matrix
-//
void choleskyDecomposition(Matrix &A, Matrix &L) {
int dim = A.rows;
assert(A.cols == dim);
@@ -226,17 +223,7 @@ 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
-//
+// Compute L (lower triangular) such that A + U U^T = L L^T.
void choleskyDecompositionStabilized(Matrix &A, Matrix &L,
vector<Integer> &stabilizeIndices,
vector<Real> &stabilizeLambdas,
@@ -259,6 +246,8 @@ void choleskyDecompositionStabilized(Matrix &A, Matrix &L,
L.elements[i + j*dim] = 0;
}
+// B := L^{-1} B, where L is lower-triangular and B is a matrix
+// pointed to by b
void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb) {
int dim = L.rows;
assert(L.cols == dim);
@@ -267,11 +256,14 @@ void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb) {
dim, bcols, 1, &L.elements[0], dim, b, ldb);
}
+// b := L^{-1} b, where L is lower-triangular
void lowerTriangularSolve(Matrix &L, Vector &b) {
assert(static_cast<int>(b.size()) == L.rows);
lowerTriangularSolve(L, &b[0], 1, b.size());
}
+// B := L^{-T} B, where L is lower-triangular and B is a matrix
+// pointed to by b
void lowerTriangularTransposeSolve(Matrix &L, Real *b, int bcols, int ldb) {
int dim = L.rows;
assert(L.cols == dim);
@@ -280,16 +272,13 @@ void lowerTriangularTransposeSolve(Matrix &L, Real *b, int bcols, int ldb) {
dim, bcols, 1, &L.elements[0], dim, b, ldb);
}
+// b := L^{-T} b, where L is lower-triangular
void lowerTriangularTransposeSolve(Matrix &L, Vector &b) {
assert(static_cast<int>(b.size()) == L.rows);
lowerTriangularTransposeSolve(L, &b[0], 1, b.size());
}
// X := ACholesky^{-1 T} ACholesky^{-1} X = A^{-1} X
-// Inputs:
-// - ACholesky : dim x dim lower triangular matrix
-// - X : dim x dim matrix
-//
void matrixSolveWithCholesky(Matrix &ACholesky, Matrix &X) {
int dim = X.rows;
assert(X.cols == dim);
diff --git a/src/Matrix.h b/src/Matrix.h
index 68c5a7b..fd1fe59 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -20,10 +20,18 @@
using std::ostream;
using std::vector;
+// A matrix M with Real entries
class Matrix {
public:
int rows;
int cols;
+ // Elements of M in row-major order:
+ //
+ // elements = { M_{0,0}, ..., M_{cols-1,0},
+ // M_{0,1}, ..., M_{cols-1,1},
+ // ...
+ // M_{0,rows-1}, ..., M_{cols-1,rows-1} }
+ //
Vector elements;
Matrix(int rows = 0, int cols = 0):
@@ -39,28 +47,21 @@ class Matrix {
return elements[r + c*rows];
}
+ // M := 0
void setZero() {
fillVector(elements, 0);
}
+ // M += c*I, where I is the identity and c is a constant
void addDiagonal(const Real &c) {
+ // ensure M is square
assert(rows == cols);
for (int i = 0; i < rows; i++)
elt(i, i) += c;
}
- void addColumn() {
- elements.resize(elements.size() + rows);
- cols++;
- }
-
- void setCols(int c) {
- elements.resize(rows*c);
- cols = c;
- }
-
- void setRowsCols(int r, int c) {
+ void resize(int r, int c) {
elements.resize(r*c);
rows = r;
cols = c;
@@ -78,6 +79,7 @@ class Matrix {
}
}
+ // M := A
void copyFrom(const Matrix &A) {
assert(rows == A.rows);
assert(cols == A.cols);
@@ -86,21 +88,25 @@ class Matrix {
elements[i] = A.elements[i];
}
+ // M := M + A
void operator+=(const Matrix &A) {
for (unsigned int i = 0; i < elements.size(); i++)
elements[i] += A.elements[i];
}
+ // M := M - A
void operator-=(const Matrix &A) {
for (unsigned int i = 0; i < elements.size(); i++)
elements[i] -= A.elements[i];
}
+ // M := c*M, where c is a constant
void operator*=(const Real &c) {
for (unsigned int i = 0; i < elements.size(); i++)
elements[i] *= c;
}
+ // The maximum absolute value of the elemnts of M
Real maxAbs() const {
return maxAbsVector(elements);
}
@@ -110,9 +116,6 @@ class Matrix {
ostream& operator<<(ostream& os, const Matrix& a);
-// B := A^T
-void transpose(const Matrix &A, Matrix &B);
-
// C := alpha*A*B + beta*C
void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B,
Real beta, Matrix &C);
@@ -151,7 +154,7 @@ void matrixEigenvalues(Matrix &A, Vector &workspace, Vector &eigenvalues);
// Minimum eigenvalue of A, via the QR method
// Inputs:
-// A : n x n Matrix (will be overwritten)
+// A : n x n Matrix (overwritten)
// eigenvalues : Vector of length n
// workspace : Vector of lenfth 3*n-1 (temporary workspace)
Real minEigenvalue(Matrix &A, Vector &workspace, Vector &eigenvalues);
@@ -160,31 +163,89 @@ Real minEigenvalue(Matrix &A, Vector &workspace, Vector &eigenvalues);
// for use with 'solveWithLUDecomposition'
void LUDecomposition(Matrix &A, vector<Integer> &pivots);
+// b := A^{-1} b, where LU and pivots encode the LU decomposition of A
void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b);
// L (lower triangular) such that A = L L^T
-// Inputs:
+// Input:
// - A : dim x dim symmetric matrix
-// - L : dim x dim lower-triangular matrix
+// Output:
+// - L : dim x dim lower-triangular matrix (overwritten)
void choleskyDecomposition(Matrix &A, Matrix &L);
+// Compute L (lower triangular) such that A + U U^T = L L^T. Here,
+// the 'update' matrix U has columns given by
+//
+// U = ( Lambda_{p_1} e_{p_1}, ..., Lambda_{p_M} e_{p_M} )
+//
+// where e_p is a unit vector in the p-th direction and the
+// Lambda_{p_m} are constants. If p_m appears above, we say the
+// direction p_m has been `stabilized.'
+//
+// We choose which direction to stabilize by comparing diagonal
+// entries A_{ii} encountered during the Cholesky decomposition to
+// stabilizeThreshold*Lambda_GM, where Lambda_GM is the geometric mean
+// of the diagonal entries L computed so far. Smaller
+// stabilizeThreshold means fewer directions will be stabilized.
+// Larger stabilizeThreshold means more directions will be stabilized.
+//
+// Input:
+// - A
+// - stabilizeThreshold: parameter for deciding which directions to
+// stabilize
+// Output:
+// - L (overwritten)
+// - stabilizeIndices: a list of directions which have been
+// stabilized (overwritten)
+// - stabilizeLambdas: a list of corresponding Lambdas (overwritten)
+//
void choleskyDecompositionStabilized(Matrix &A, Matrix &L,
vector<Integer> &stabilizeIndices,
vector<Real> &stabilizeLambdas,
const double stabilizeThreshold);
+// B := L^{-1} B, where L is lower-triangular and B is a matrix
+// pointed to by b
+//
+// Input:
+// - L
+// - b, a pointer to the top-left element of B
+// - bcols, the number of columns of B
+// - ldb, the distance in pointer-space between the first elements of
+// each row of B
+// Output:
+// - elements of B, which are values pointed to by b are overwritten
+//
+// (The use of pointers and ldb is to allow using this function for
+// submatrices of a larger matrix.)
void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb);
+// b := L^{-1} b, where L is lower-triangular
void lowerTriangularSolve(Matrix &L, Vector &b);
+// B := L^{-T} B, where L is lower-triangular and B is a matrix
+// pointed to by b
+//
+// Input:
+// - L
+// - b, a pointer to the top-left element of B
+// - bcols, the number of columns of B
+// - ldb, the distance in pointer-space between the first elements of
+// each row of B
+// Output:
+// - elements of B, which are values pointed to by b are overwritten
+//
void lowerTriangularTransposeSolve(Matrix &L, Real *b, int bcols, int ldb);
+// b := L^{-T} b, where L is lower-triangular
void lowerTriangularTransposeSolve(Matrix &L, Vector &b);
// X := ACholesky^{-1 T} ACholesky^{-1} X = A^{-1} X
// Inputs:
-// - ACholesky : dim x dim lower triangular matrix
-// - X : dim x dim matrix
+// - ACholesky : dim x dim lower triangular matrix (constant)
+// Output:
+// - X : dim x dim matrix (overwritten)
+//
void matrixSolveWithCholesky(Matrix &ACholesky, Matrix &X);
#endif // SDPB_MATRIX_H_
diff --git a/src/SDPSolver.cpp b/src/SDPSolver.cpp
index 210e0c0..5a8542a 100644
--- a/src/SDPSolver.cpp
+++ b/src/SDPSolver.cpp
@@ -565,7 +565,7 @@ Real correctorCenteringParameter(const SDPSolverParameters ¶meters,
// - MCholesky = L, the Cholesky decomposition of M (M itself is not needed)
// - dM, a BlockDiagonalMatrix with the same structure as M
// Workspace:
-// - MInvDM
+// - MInvDM (NB: overwritten when computing minEigenvalue)
// - eigenvalues, a Vector of eigenvalues for each block of M
// - workspace, a vector of Vectors needed by the minEigenvalue function
// Output:
@@ -686,8 +686,8 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
int startIndex = schurStabilizeIndices[j][0];
// set the dimensions of U_j
- stabilizeBlocks[j].setRowsCols(SchurComplement.blocks[j].rows - startIndex,
- schurStabilizeIndices[j].size());
+ stabilizeBlocks[j].resize(SchurComplement.blocks[j].rows - startIndex,
+ schurStabilizeIndices[j].size());
// set U_j = 0
stabilizeBlocks[j].setZero();
// for each column of U_j add Lambda_p in the row (p - startIndex)
@@ -727,19 +727,19 @@ void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &Bilin
//
// Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
//
- // Where B' = (B U). We think of Q as containing four-blocks called
+ // Where B' = (B U). We think of Q as containing four blocks called
// Upper/Lower-Left/Right.
// Set the dimensions of Q
- Q.setRowsCols(offDiagonalColumns, offDiagonalColumns);
+ Q.resize(offDiagonalColumns, offDiagonalColumns);
Q.setZero();
- // At this point, SchurOffDiagonal = L'^{-1} B.
+ // Here, SchurOffDiagonal = L'^{-1} B.
//
// UpperLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal
matrixSquareIntoBlock(SchurOffDiagonal, Q, 0, 0);
- // At this point, stabilizeBlocks the blocks of V = L'^{-1} U.
+ // Here, stabilizeBlocks contains the blocks of V = L'^{-1} U.
//
// LowerRight(Q) = V^T V - 1
for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
@@ -894,6 +894,9 @@ void SDPSolver::computeSearchDirection(const Real &beta,
dY *= -1;
}
+/***********************************************************************/
+// The main solver loop
+
SDPSolverTerminateReason SDPSolver::run(const path checkpointFile) {
Real primalStepLength;
Real dualStepLength;
diff --git a/src/main.cpp b/src/main.cpp
index a11b3e7..95aae9f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -39,10 +39,17 @@ int solveSDP(const path &sdpFile,
const path &outFile,
const path &checkpointFile,
SDPSolverParameters parameters) {
+ // Set the default precision of all Real numbers to that specified
+ // by the 'precision' parameter.
setDefaultPrecision(parameters.precision);
+
+ // Set cout to print the appropriate number of digits
cout.precision(min(static_cast<int>(parameters.precision * 0.31 + 5), 30));
+
// Ensure all the Real parameters have the appropriate precision
parameters.resetPrecision();
+
+ // Set the maximum number of threads for parallel loops
omp_set_num_threads(parameters.maxThreads);
cout << "SDPB started at " << second_clock::local_time() << endl;
@@ -53,8 +60,8 @@ int solveSDP(const path &sdpFile,
cout << "\nParameters:\n";
cout << parameters << endl;
- const SDP sdp = readBootstrapSDP(sdpFile);
- SDPSolver solver(sdp, parameters);
+ // Read an SDP from sdpFile and create a solver for it
+ SDPSolver solver(readBootstrapSDP(sdpFile), parameters);
if (exists(checkpointFile))
solver.loadCheckpoint(checkpointFile);
--
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