[sdpb] 108/233: Ran cpplint on some files
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:26 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 ff6efd47baed988a3d5054ffa10f0b8eee64f3f1
Author: David Simmons-Duffin <dsd at minerva.sns.ias.edu>
Date: Tue Nov 18 17:27:49 2014 -0500
Ran cpplint on some files
---
src/BlockDiagonalMatrix.cpp | 52 ++++++++++++++++++--------------
src/BlockDiagonalMatrix.h | 51 +++++++++++++++++++++----------
src/Matrix.cpp | 73 ++++++++++++++++++++++++++-------------------
src/Matrix.h | 36 ++++++++++------------
src/SDP.cpp | 23 +++++++-------
src/parse.cpp | 7 +++--
6 files changed, 142 insertions(+), 100 deletions(-)
diff --git a/src/BlockDiagonalMatrix.cpp b/src/BlockDiagonalMatrix.cpp
index 38ea457..5a7e978 100644
--- a/src/BlockDiagonalMatrix.cpp
+++ b/src/BlockDiagonalMatrix.cpp
@@ -6,6 +6,8 @@
//=======================================================================
+#include <algorithm>
+#include <vector>
#include "BlockDiagonalMatrix.h"
ostream& operator<<(ostream& os, const BlockDiagonalMatrix& A) {
@@ -32,7 +34,7 @@ Real frobeniusProductSymmetric(const BlockDiagonalMatrix &A,
}
return result;
}
-
+
// (X + dX) . (Y + dY), where X, dX, Y, dY are symmetric
// BlockDiagonalMatrices and '.' is the Frobenius product.
//
@@ -43,7 +45,8 @@ Real frobeniusProductOfSums(const BlockDiagonalMatrix &X,
Real result = 0;
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < X.blocks.size(); b++) {
- Real f = frobeniusProductOfSums(X.blocks[b], dX.blocks[b], Y.blocks[b], dY.blocks[b]);
+ Real f = frobeniusProductOfSums(X.blocks[b], dX.blocks[b],
+ Y.blocks[b], dY.blocks[b]);
#pragma omp critical
{
result += f;
@@ -68,7 +71,8 @@ void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A,
blockDiagonalMatrixScaleMultiplyAdd(1, A, B, 0, C);
}
-void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L) {
+void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A,
+ BlockDiagonalMatrix &L) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < A.blocks.size(); b++)
lowerTriangularInverseCongruence(A.blocks[b], L.blocks[b]);
@@ -80,15 +84,19 @@ void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A, BlockDiagonalMatri
// eigenvalues : vector<Vector> of length A.blocks.size()
// workspace : vector<Vector> of length A.blocks.size()
//
-Real minEigenvalue(BlockDiagonalMatrix &A, vector<Vector> &workspace, vector<Vector> &eigenvalues) {
+Real minEigenvalue(BlockDiagonalMatrix &A,
+ vector<Vector> &workspace,
+ vector<Vector> &eigenvalues) {
assert(A.blocks.size() == eigenvalues.size());
assert(A.blocks.size() == workspace.size());
- // TODO: get rid of this hack
+ // TODO(davidsd): get rid of this hack
Real lambdaMin = 1e100;
#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]);
+ Real minBlockLambda = minEigenvalue(A.blocks[b],
+ workspace[b],
+ eigenvalues[b]);
#pragma omp critical
{
lambdaMin = min(lambdaMin, minBlockLambda);
@@ -107,14 +115,14 @@ void choleskyDecomposition(BlockDiagonalMatrix &A,
void choleskyDecompositionStabilized(BlockDiagonalMatrix &A,
BlockDiagonalMatrix &L,
- vector<vector<Integer> > &schurStabilizeIndices,
- vector<vector<Real> > &schurStabilizeLambdas,
+ vector<vector<Integer> > &stabilizeIndices,
+ vector<vector<Real> > &stabilizeLambdas,
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],
+ stabilizeIndices[b],
+ stabilizeLambdas[b],
stabilizeThreshold);
}
@@ -125,26 +133,26 @@ void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
matrixSolveWithCholesky(ACholesky.blocks[b], X.blocks[b]);
}
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, Matrix &B) {
- #pragma omp parallel for schedule(dynamic)
- for (unsigned int b = 0; b < L.blocks.size(); b++)
- lowerTriangularSolve(L.blocks[b], &B.elt(L.blockStartIndices[b], 0), B.cols, B.rows);
-}
-
-void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, Matrix &B) {
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
+ Matrix &B) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < L.blocks.size(); b++)
- lowerTriangularTransposeSolve(L.blocks[b], &B.elt(L.blockStartIndices[b], 0), B.cols, B.rows);
+ lowerTriangularSolve(L.blocks[b], &B.elt(L.blockStartIndices[b], 0),
+ B.cols, B.rows);
}
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, Vector &v) {
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
+ Vector &v) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < L.blocks.size(); b++)
- lowerTriangularSolve(L.blocks[b], &v[L.blockStartIndices[b]], 1, v.size());
+ lowerTriangularSolve(L.blocks[b], &v[L.blockStartIndices[b]],
+ 1, v.size());
}
-void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, Vector &v) {
+void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L,
+ Vector &v) {
#pragma omp parallel for schedule(dynamic)
for (unsigned int b = 0; b < L.blocks.size(); b++)
- lowerTriangularTransposeSolve(L.blocks[b], &v[L.blockStartIndices[b]], 1, v.size());
+ lowerTriangularTransposeSolve(L.blocks[b], &v[L.blockStartIndices[b]],
+ 1, v.size());
}
diff --git a/src/BlockDiagonalMatrix.h b/src/BlockDiagonalMatrix.h
index 691e648..b8a45a1 100644
--- a/src/BlockDiagonalMatrix.h
+++ b/src/BlockDiagonalMatrix.h
@@ -25,7 +25,7 @@ public:
vector<Matrix> blocks;
vector<int> blockStartIndices;
- BlockDiagonalMatrix(const vector<int> &blockSizes):
+ explicit BlockDiagonalMatrix(const vector<int> &blockSizes):
dim(0) {
for (unsigned int i = 0; i < blockSizes.size(); i++) {
blocks.push_back(Matrix(blockSizes[i], blockSizes[i]));
@@ -77,7 +77,9 @@ public:
Real maxAbs() const {
Real max = 0;
- for (vector<Matrix>::const_iterator b = blocks.begin(); b != blocks.end(); b++) {
+ for (vector<Matrix>::const_iterator b = blocks.begin();
+ b != blocks.end();
+ b++) {
const Real tmp = b->maxAbs();
if (tmp > max)
max = tmp;
@@ -86,12 +88,11 @@ public:
}
friend ostream& operator<<(ostream& os, const BlockDiagonalMatrix& A);
-
};
Real frobeniusProductSymmetric(const BlockDiagonalMatrix &A,
const BlockDiagonalMatrix &B);
-
+
// (X + dX) . (Y + dY), where X, dX, Y, dY are symmetric
// BlockDiagonalMatrices and '.' is the Frobenius product.
//
@@ -105,19 +106,37 @@ void blockDiagonalMatrixScaleMultiplyAdd(Real alpha,
BlockDiagonalMatrix &B,
Real beta,
BlockDiagonalMatrix &C);
-void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A, BlockDiagonalMatrix &B, BlockDiagonalMatrix &C);
-void lowerTriangularInverseCongruence(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L);
-Real minEigenvalue(BlockDiagonalMatrix &A, vector<Vector> &workspace, vector<Vector> &eigenvalues);
-void choleskyDecomposition(BlockDiagonalMatrix &A, BlockDiagonalMatrix &L);
+
+void blockDiagonalMatrixMultiply(BlockDiagonalMatrix &A,
+ BlockDiagonalMatrix &B,
+ BlockDiagonalMatrix &C);
+
+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,
+ vector<vector<Integer> > &stabilizeIndices,
+ vector<vector<Real> > &stabilizeLambdas,
const double stabilizeThreshold);
-void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky, BlockDiagonalMatrix &X);
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, Matrix &B);
-void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, Matrix &B);
-void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L, Vector &v);
-void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L, Vector &v);
-#endif // SDPB_MATRIX_H_
+void blockMatrixSolveWithCholesky(BlockDiagonalMatrix &ACholesky,
+ BlockDiagonalMatrix &X);
+
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
+ Matrix &B);
+
+void blockMatrixLowerTriangularSolve(BlockDiagonalMatrix &L,
+ Vector &v);
+
+void blockMatrixLowerTriangularTransposeSolve(BlockDiagonalMatrix &L,
+ Vector &v);
+
+#endif // SDPB_BLOCKDIAGONALMATRIX_H_
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 7023ce1..713486a 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -6,6 +6,7 @@
//=======================================================================
+#include <vector>
#include "Matrix.h"
ostream& operator<<(ostream& os, const Matrix& a) {
@@ -13,7 +14,7 @@ ostream& operator<<(ostream& os, const Matrix& a) {
for (int r = 0; r < a.rows; r++) {
os << "{";
for (int c = 0; c < a.cols; c++) {
- os << a.elt(r,c);
+ os << a.elt(r, c);
if (c < a.cols-1)
os << ", ";
}
@@ -32,12 +33,13 @@ void transpose(const Matrix &A, Matrix &B) {
for (int n = 0; n < A.cols; n++)
for (int m = 0; m < A.rows; m++)
- B.elt(n,m) = A.elt(m,n);
+ 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) {
+void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B,
+ Real beta, Matrix &C) {
assert(A.cols == B.rows);
assert(A.rows == C.rows);
assert(B.cols == C.cols);
@@ -59,13 +61,13 @@ void matrixMultiply(Matrix &A, Matrix &B, Matrix &C) {
void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol) {
assert(bRow + A.cols <= B.rows);
assert(bCol + A.cols <= B.cols);
-
+
#pragma omp parallel for schedule(dynamic)
for (int c = 0; c < A.cols; c++) {
for (int r = 0; r <= c; r++) {
Real tmp = 0;
for (int p = 0; p < A.rows; p++)
- tmp += A.elt(p,r) * A.elt(p,c);
+ tmp += A.elt(p, r) * A.elt(p, c);
B.elt(bRow + r, bCol + c) = tmp;
if (r != c)
B.elt(bRow + c, bCol + r) = tmp;
@@ -91,30 +93,32 @@ void lowerTriangularInverseCongruence(Matrix &A, Matrix &L) {
// y := alpha A x + beta y
//
-void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x, Real beta, Vector &y) {
- assert(A.cols <= (int)x.size());
- assert(A.rows <= (int)y.size());
+void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x,
+ Real beta, Vector &y) {
+ assert(A.cols <= static_cast<int>(x.size()));
+ assert(A.rows <= static_cast<int>(y.size()));
#pragma omp parallel for schedule(static)
for (int p = 0; p < A.rows; p++) {
Real tmp = 0;
for (int n = 0; n < A.cols; n++)
- tmp += A.elt(p,n)*x[n];
+ tmp += A.elt(p, n)*x[n];
y[p] = alpha*tmp + beta*y[p];
}
}
// y := alpha A^T x + beta y
//
-void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x, Real beta, Vector &y) {
- assert(A.cols <= (int)y.size());
- assert(A.rows <= (int)x.size());
+void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x,
+ Real beta, Vector &y) {
+ assert(A.cols <= static_cast<int>(y.size()));
+ assert(A.rows <= static_cast<int>(x.size()));
#pragma omp parallel for schedule(static)
for (int n = 0; n < A.cols; n++) {
Real tmp = 0;
for (int p = 0; p < A.rows; p++)
- tmp += A.elt(p,n)*x[p];
+ tmp += A.elt(p, n)*x[p];
y[n] = alpha*tmp + beta*y[n];
}
}
@@ -127,12 +131,12 @@ Real frobeniusProductSymmetric(const Matrix &A, const Matrix &B) {
Real result = 0;
for (int c = 0; c < A.cols; c++)
for (int r = 0; r < c ; r++)
- result += A.elt(r,c)*B.elt(r,c);
+ result += A.elt(r, c)*B.elt(r, c);
result *= 2;
for (int r = 0; r < A.rows; r++)
- result += A.elt(r,r)*B.elt(r,r);
-
+ result += A.elt(r, r)*B.elt(r, r);
+
return result;
}
@@ -145,11 +149,11 @@ Real frobeniusProductOfSums(const Matrix &X, const Matrix &dX,
for (int c = 0; c < X.cols; c++)
for (int r = 0; r < c; r++)
- result += (X.elt(r,c) + dX.elt(r,c)) * (Y.elt(r,c) + dY.elt(r,c));
+ result += (X.elt(r, c) + dX.elt(r, c)) * (Y.elt(r, c) + dY.elt(r, c));
result *= 2;
for (int r = 0; r < X.rows; r++)
- result += (X.elt(r,r) + dX.elt(r,r)) * (Y.elt(r,r) + dY.elt(r,r));
+ result += (X.elt(r, r) + dX.elt(r, r)) * (Y.elt(r, r) + dY.elt(r, r));
return result;
}
@@ -162,12 +166,13 @@ Real frobeniusProductOfSums(const Matrix &X, const Matrix &dX,
//
void matrixEigenvalues(Matrix &A, Vector &workspace, Vector &eigenvalues) {
assert(A.rows == A.cols);
- assert((int)eigenvalues.size() == A.rows);
- assert((int)workspace.size() == 3*A.rows - 1);
+ assert(static_cast<int>(eigenvalues.size()) == A.rows);
+ assert(static_cast<int>(workspace.size()) == 3*A.rows - 1);
Integer info;
Integer workSize = workspace.size();
- Rsyev("NoEigenvectors", "LowerTriangular", A.rows, &A.elements[0], A.rows, &eigenvalues[0], &workspace[0], workSize, &info);
+ Rsyev("NoEigenvectors", "LowerTriangular", A.rows, &A.elements[0],
+ A.rows, &eigenvalues[0], &workspace[0], workSize, &info);
assert(info == 0);
}
@@ -193,7 +198,8 @@ void LUDecomposition(Matrix &A, vector<Integer> &pivots) {
void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b) {
Integer info;
- Rgetrs("NoTranspose", LU.rows, 1, &LU.elements[0], LU.rows, &pivots[0], &b[0], b.size(), &info);
+ Rgetrs("NoTranspose", LU.rows, 1, &LU.elements[0], LU.rows,
+ &pivots[0], &b[0], b.size(), &info);
assert(info == 0);
}
@@ -220,15 +226,21 @@ 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
+// 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
+// - 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, const double stabilizeThreshold) {
+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);
@@ -237,7 +249,8 @@ 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, stabilizeThreshold);
+ RpotrfStabilized("Lower", dim, &L.elements[0], dim, &info,
+ stabilizeIndices, stabilizeLambdas, stabilizeThreshold);
assert(info == 0);
// Set the upper-triangular part of the L to zero
@@ -255,7 +268,7 @@ void lowerTriangularSolve(Matrix &L, Real *b, int bcols, int ldb) {
}
void lowerTriangularSolve(Matrix &L, Vector &b) {
- assert((int) b.size() == L.rows);
+ assert(static_cast<int>(b.size()) == L.rows);
lowerTriangularSolve(L, &b[0], 1, b.size());
}
@@ -268,7 +281,7 @@ void lowerTriangularTransposeSolve(Matrix &L, Real *b, int bcols, int ldb) {
}
void lowerTriangularTransposeSolve(Matrix &L, Vector &b) {
- assert((int) b.size() == L.rows);
+ assert(static_cast<int>(b.size()) == L.rows);
lowerTriangularTransposeSolve(L, &b[0], 1, b.size());
}
diff --git a/src/Matrix.h b/src/Matrix.h
index fe54130..f61b7bb 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -12,6 +12,7 @@
#include <assert.h>
#include <iostream>
#include <ostream>
+#include <vector>
#include "omp.h"
#include "types.h"
#include "Vector.h"
@@ -19,9 +20,6 @@
using std::ostream;
using std::vector;
-const Real CHOLESKY_STABILIZE_THRESHOLD = 1e-8;
-const Real BASIC_ROW_THRESHOLD = 1e-4;
-
class Matrix {
public:
int rows;
@@ -49,7 +47,7 @@ class Matrix {
assert(rows == cols);
for (int i = 0; i < rows; i++)
- elt(i,i) += c;
+ elt(i, i) += c;
}
void addColumn() {
@@ -73,9 +71,9 @@ class Matrix {
for (int r = 0; r < rows; r++) {
for (int c = 0; c < r; c++) {
- Real tmp = (elt(r,c) + elt(c,r))/2;
- elt(r,c) = tmp;
- elt(c,r) = tmp;
+ Real tmp = (elt(r, c) + elt(c, r))/2;
+ elt(r, c) = tmp;
+ elt(c, r) = tmp;
}
}
}
@@ -91,7 +89,7 @@ class Matrix {
void operator+=(const Matrix &A) {
for (unsigned int i = 0; i < elements.size(); i++)
elements[i] += A.elements[i];
- }
+ }
void operator-=(const Matrix &A) {
for (unsigned int i = 0; i < elements.size(); i++)
@@ -107,14 +105,6 @@ class Matrix {
return maxAbsVector(elements);
}
- void swapCols(int c1, int c2) {
- for (int r = 0; r < rows; r++) {
- Real tmp = elt(r, c1);
- elt(r,c1) = elt(r,c2);
- elt(r,c1) = tmp;
- }
- }
-
friend ostream& operator<<(ostream& os, const Matrix& a);
};
@@ -124,7 +114,8 @@ ostream& operator<<(ostream& os, const Matrix& a);
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);
+void matrixScaleMultiplyAdd(Real alpha, Matrix &A, Matrix &B,
+ Real beta, Matrix &C);
// C := A*B
void matrixMultiply(Matrix &A, Matrix &B, Matrix &C);
@@ -136,10 +127,12 @@ void matrixSquareIntoBlock(Matrix &A, Matrix &B, int bRow, int bCol);
void lowerTriangularInverseCongruence(Matrix &A, Matrix &L);
// y := alpha A x + beta y
-void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x, Real beta, Vector &y);
+void vectorScaleMatrixMultiplyAdd(Real alpha, Matrix &A, Vector &x,
+ Real beta, Vector &y);
// y := alpha A^T x + beta y
-void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x, Real beta, Vector &y);
+void vectorScaleMatrixMultiplyTransposeAdd(Real alpha, Matrix &A, Vector &x,
+ Real beta, Vector &y);
// Frobenius product Tr(A^T B) where A and B are symmetric matrices
Real frobeniusProductSymmetric(const Matrix &A, const Matrix &B);
@@ -175,7 +168,10 @@ 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, const double stabilizeThreshold);
+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/SDP.cpp b/src/SDP.cpp
index 476a1b4..687dabf 100644
--- a/src/SDP.cpp
+++ b/src/SDP.cpp
@@ -6,6 +6,7 @@
//=======================================================================
+#include <vector>
#include "SDP.h"
Matrix sampleBilinearBasis(const int maxDegree,
@@ -23,7 +24,8 @@ Matrix sampleBilinearBasis(const int maxDegree,
return b;
}
-SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m) {
+SampledMatrixPolynomial
+samplePolynomialVectorMatrix(const PolynomialVectorMatrix &m) {
SampledMatrixPolynomial s;
assert(m.rows == m.cols);
@@ -32,14 +34,14 @@ SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatri
int numSamples = s.degree + 1;
int numConstraints = numSamples * s.dim * (s.dim + 1)/2;
- int vectorDim = m.elt(0,0).size();
+ int vectorDim = m.elt(0, 0).size();
// The first element of each vector multiplies the constant 1
s.constraintConstants = Vector(numConstraints);
// The rest multiply decision variables
s.constraintMatrix = Matrix(numConstraints, vectorDim - 1);
-
+
int p = 0;
for (int c = 0; c < s.dim; c++) {
for (int r = c; r < s.dim; r++) {
@@ -47,7 +49,7 @@ SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatri
Real x = m.samplePoints[k];
Real scale = m.sampleScalings[k];
- s.constraintConstants[p] = scale*m.elt(r,c)[0](x);
+ s.constraintConstants[p] = scale*m.elt(r, c)[0](x);
for (int n = 1; n < vectorDim; n++)
s.constraintMatrix.elt(p, n-1) = -scale*m.elt(r, c)[n](x);
@@ -57,17 +59,19 @@ SampledMatrixPolynomial samplePolynomialVectorMatrix(const PolynomialVectorMatri
}
int delta1 = s.degree/2;
- s.bilinearBases.push_back(sampleBilinearBasis(delta1, numSamples, m.bilinearBasis,
+ s.bilinearBases.push_back(sampleBilinearBasis(delta1, numSamples,
+ m.bilinearBasis,
m.samplePoints,
m.sampleScalings));
int delta2 = (s.degree - 1)/2;
if (delta2 >= 0)
- s.bilinearBases.push_back(sampleBilinearBasis(delta2, numSamples, m.bilinearBasis,
+ s.bilinearBases.push_back(sampleBilinearBasis(delta2, numSamples,
+ m.bilinearBasis,
m.samplePoints,
- multiplyVectors(m.samplePoints, m.sampleScalings)));
+ multiplyVectors(m.samplePoints,
+ m.sampleScalings)));
return s;
-
}
SDP bootstrapSDP(const Vector &objective,
@@ -92,7 +96,6 @@ SDP bootstrapSDP(const Vector &objective,
for (vector<SampledMatrixPolynomial>::const_iterator s = sampledMatrixPols.begin();
s != sampledMatrixPols.end();
s++) {
-
vector<int> blocks;
for (vector<Matrix>::const_iterator b = s->bilinearBases.begin();
b != s->bilinearBases.end();
@@ -107,7 +110,7 @@ SDP bootstrapSDP(const Vector &objective,
for (int n = 0; n < s->constraintMatrix.cols; n++)
sdp.FreeVarMatrix.elt(p, n) = s->constraintMatrix.elt(k, n);
}
- assert(p == (int)sdp.primalObjective.size());
+ assert(p == static_cast<int>(sdp.primalObjective.size()));
sdp.initializeConstraintIndices();
return sdp;
diff --git a/src/parse.cpp b/src/parse.cpp
index 0c23dfc..c6505ef 100644
--- a/src/parse.cpp
+++ b/src/parse.cpp
@@ -19,7 +19,8 @@ using tinyxml2::XMLElement;
using boost::filesystem::path;
template <class T>
-vector<T> parseMany(const char *name, T(*parse)(XMLElement *), XMLElement *elt) {
+vector<T> parseMany(const char *name, T(*parse)(XMLElement *),
+ XMLElement *elt) {
XMLElement *e;
vector<T> v;
for (e = elt->FirstChildElement(name);
@@ -56,7 +57,9 @@ SampledMatrixPolynomial parseSampledMatrixPolynomial(XMLElement *xml) {
s.degree = parseInt(xml->FirstChildElement("degree"));
s.constraintMatrix = parseMatrix(xml->FirstChildElement("constraintMatrix"));
s.constraintConstants = parseVector(xml->FirstChildElement("constraintConstants"));
- s.bilinearBases = parseMany("bilinearBasisMatrix", parseMatrix, xml->FirstChildElement("bilinearBases"));
+ s.bilinearBases = parseMany("bilinearBasisMatrix",
+ parseMatrix,
+ xml->FirstChildElement("bilinearBases"));
return s;
}
--
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