[sdpb] 02/233: Consolidated the code into main for now; wrote an initial schur complement routine; testing for block congruence, cholesky decomp.
Tobias Hansen
thansen at moszumanska.debian.org
Thu Mar 9 04:06:11 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 72f8c90ab4c98bf13726973fedad88a281ead8f8
Author: David Simmons-Duffin <dsd at romulus.sns.ias.edu>
Date: Fri Jun 20 00:24:36 2014 -0400
Consolidated the code into main for now; wrote an initial schur complement routine; testing for block congruence, cholesky decomp.
---
BlockDenseMatrix.h | 39 -----
DenseMatrix.h | 90 ----------
LowRankConstraintMatrix.h | 29 ----
Makefile | 2 +-
SquareMatrix.h | 56 ++++++
main.cpp | 426 ++++++++++++++++++++++++++++++++++++++++++----
6 files changed, 447 insertions(+), 195 deletions(-)
diff --git a/BlockDenseMatrix.h b/BlockDenseMatrix.h
deleted file mode 100644
index 8ab8438..0000000
--- a/BlockDenseMatrix.h
+++ /dev/null
@@ -1,39 +0,0 @@
-#ifndef SDP_BOOTSTRAP_BLOCKDENSEMATRIX_H_
-#define SDP_BOOTSTRAP_BLOCKDENSEMATRIX_H_
-
-#include <vector>
-#include "types.h"
-#include "DenseMatrix.h"
-
-using std::vector;
-
-class BlockDenseMatrix {
-public:
- vector<Real> diagonalPart;
- vector<DenseMatrix> blocks;
-
- BlockDenseMatrix(const vector<Real> &diagonalPart, const vector<int> &blockSizes):
- diagonalPart(diagonalPart) {
- for (unsigned int i = 0; i < blockSizes.size(); i++) {
- blocks.push_back(DenseMatrix(blockSizes[i]));
- }
- }
-
- void inverseCholeskyAndInverse(BlockDenseMatrix &invCholesky,
- BlockDenseMatrix &inverse,
- BlockDenseMatrix &workSpace) {
-
- for (unsigned int i = 0; i < diagonalPart.size(); i++) {
- Real d = diagonalPart[i];
- invCholesky.diagonalPart[i] = 1/sqrt(d);
- inverse.diagonalPart[i] = 1/d;
- }
-
- for (unsigned int b = 0; b < blocks.size(); b++) {
- blocks[b].inverseCholeskyAndInverse(invCholesky.blocks[b], inverse.blocks[b], workSpace.blocks[b]);
- }
- }
-
-};
-
-#endif // SDP_BOOTSTRAP_BLOCKDENSEMATRIX_H_
diff --git a/DenseMatrix.h b/DenseMatrix.h
deleted file mode 100644
index 3bc927e..0000000
--- a/DenseMatrix.h
+++ /dev/null
@@ -1,90 +0,0 @@
-#ifndef SDP_BOOTSTRAP_DENSEMATRIX_H_
-#define SDP_BOOTSTRAP_DENSEMATRIX_H_
-
-#include <vector>
-#include "types.h"
-
-using std::vector;
-
-class DenseMatrix {
-public:
- int dim;
-
- // The elements of the matrix, in column-major order, for
- // consistency with lapack's conventions.
- vector<Real> elements;
-
- DenseMatrix(int dim):
- dim(dim),
- elements(vector<Real>(dim*dim, 0)) {}
-
- inline Real elem(int i, int j) const {
- return elements[i + j*dim];
- }
-
- // thinking of our matrix M as a block matrix where each block has
- // size v.size() == u.size(), compute the pairing v^T M[ij] u, where
- // M[ij] denotes the i,j-th block of M
- Real blockBilinearForm(int i, int j, const vector<Real> &v, const vector<Real> &u) const {
- int blockDim = v.size();
-
- Real total = 0;
- for (int k = 0; k < blockDim; k++) {
- Real uk = u[k];
- for (int l = 0; l < blockDim; l++) {
- total += v[l] * uk * elem(i*blockDim+l, j*blockDim + k);
- }
- }
-
- return total;
- }
-
- void setToZero() {
- std::fill(elements.begin(), elements.end(), 0);
- }
-
- void setToIdentity() {
- setToZero();
- for (int i = 0; i < dim; i++)
- elements[i * (dim + 1)] = 1;
- }
-
- void choleskyDecomposition(DenseMatrix &result) {
- mpackint info;
- Real *resultArray = &result.elements[0];
- Rcopy(dim*dim, &elements[0], 1, resultArray, 1);
-
- // The lower-triangular part of result is now our cholesky matrix
- Rpotrf("Lower", dim, resultArray, dim, &info);
-
- // Set the upper-triangular part of the result to zero
- for (int j = 0; j < dim; j++) {
- for (int i = 0; i < j; i++) {
- result.elements[i + j*dim] = 0;
- }
- }
- }
-
- void inverseLowerTriangular(DenseMatrix &result) {
- result.setToIdentity();
- Rtrsm("Left", "Lower", "NoTranspose", "NonUnitDiagonal",
- dim, dim, 1, &elements[0], dim, &result.elements[0], dim);
- }
-
- void inverseCholesky(DenseMatrix &result, DenseMatrix &workSpace) {
- choleskyDecomposition(workSpace);
- workSpace.inverseLowerTriangular(result);
- }
-
- void inverseCholeskyAndInverse(DenseMatrix &invCholesky, DenseMatrix &inverse, DenseMatrix &workSpace) {
- inverseCholesky(invCholesky, workSpace);
-
- inverse.elements = invCholesky.elements;
- Rtrmm("Left", "Lower", "Transpose", "NonUnitDiag", dim, dim, 1,
- &invCholesky.elements[0], dim,
- &inverse.elements[0], dim);
- }
-
-};
-
-#endif // SDP_BOOTSTRAP_DENSEMATRIX_H_
diff --git a/LowRankConstraintMatrix.h b/LowRankConstraintMatrix.h
deleted file mode 100644
index b74406d..0000000
--- a/LowRankConstraintMatrix.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef SDP_BOOTSTRAP_LOWRANKCONSTRAINTMATRIX_H_
-#define SDP_BOOTSTRAP_LOWRANKCONSTRAINTMATRIX_H_
-
-#include <vector>
-#include "types.h"
-
-using std::pair;
-using std::vector;
-
-// A matrix of the form E_{ij} \otimes v v^T, where E_{ij} is the
-// matrix with (i,j)-th entry 1, and the rest zero, and v is a vector.
-class RankOneSingleBlockMatrix {
-public:
- int i;
- int j;
- vector<Real> vec;
-
- RankOneSingleBlockMatrix(int i, int j, const vector<Real> &vec):
- i(i), j(j), vec(vec) {}
-};
-
-class LowRankConstraintMatrix {
-public:
- vector<Real> diagonalPart;
- vector<pair<int, RankOneSingleBlockMatrix> > lowRankBlocks;
- vector<pair<int, DenseMatrix> > denseBlocks;
-};
-
-#endif // SDP_BOOTSTRAP_LOWRANKCONSTRAINTMATRIX_H_
diff --git a/Makefile b/Makefile
index 5a1fddc..315da65 100755
--- a/Makefile
+++ b/Makefile
@@ -16,7 +16,7 @@ OBJECTS = main.o mpack/Rdot.o \
mpack/Rlaev2.o mpack/Rlasr.o mpack/Rlartg.o \
mpack/Rswap.o mpack/Rsyev.o mpack/Rlansy.o \
mpack/Mutils.o
-HEADERS = types.h mpack/mblas_gmp.h mpack/mlapack_gmp.h mpack/mpack_config.h mpack/mutils_gmp.h
+HEADERS = types.h mpack/mblas_gmp.h mpack/mlapack_gmp.h mpack/mpack_config.h mpack/mutils_gmp.h SquareMatrix.h
SOURCES = $(OJBECTS:.o=.cpp)
RESULT = sdp-bootstrap
diff --git a/SquareMatrix.h b/SquareMatrix.h
new file mode 100644
index 0000000..9026201
--- /dev/null
+++ b/SquareMatrix.h
@@ -0,0 +1,56 @@
+#ifndef SDP_BOOTSTRAP_SQUAREMATRIX_H_
+#define SDP_BOOTSTRAP_SQUAREMATRIX_H_
+
+#include <vector>
+#include <assert.h>
+#include <ostream>
+#include "types.h"
+
+using std::vector;
+using std::ostream;
+
+class Matrix {
+ public:
+ int rows;
+ int cols;
+ vector<Real> elements;
+
+ Matrix(int rows, int cols):
+ rows(rows),
+ cols(cols),
+ elements(vector<Real>(rows*cols, 0)) {}
+
+ inline Real get(int r, int c) const {
+ return elements[r + c*rows];
+ }
+
+ inline void set(int r, int c, const Real &a) {
+ elements[r + c*rows] = a;
+ }
+
+ void setZero() {
+ std::fill(elements.begin(), elements.end(), 0);
+ }
+
+ void setIdentity() {
+ assert(rows == cols);
+
+ setZero();
+ for (int i = 0; i < rows; i++)
+ elements[i * (rows + 1)] = 1;
+ }
+
+ friend ostream& operator<<(ostream& os, const Matrix& a);
+};
+
+void blockMatrixCongruence(const Matrix &a, const Matrix &b, Matrix &work, Matrix &result);
+
+void choleskyDecomposition(Matrix &a, Matrix &result);
+
+void inverseLowerTriangular(Matrix &a, Matrix &result);
+
+void inverseCholesky(Matrix &a, Matrix &work, Matrix &result);
+
+void inverseCholeskyAndInverse(Matrix &a, Matrix &work, Matrix &invCholesky, Matrix &inverse);
+
+#endif // SDP_BOOTSTRAP_SQUAREMATRIX_H_
diff --git a/main.cpp b/main.cpp
index 7c2ab8e..c14cfaa 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,61 +1,415 @@
+#include <iterator>
#include <iostream>
+#include <ostream>
#include <vector>
+#include <assert.h>
#include "types.h"
-#include "DenseMatrix.h"
-#include "BlockDenseMatrix.h"
-#include "LowRankConstraintMatrix.h"
+using std::vector;
using std::cout;
using std::endl;
+using std::ostream;
-Real schurComplementPairing(const LowRankConstraintMatrix &F1,
- const LowRankConstraintMatrix &F2,
- const BlockDenseMatrix &X,
- const BlockDenseMatrix &ZInv) {
- Real total = 0;
+class Matrix {
+ public:
+ int rows;
+ int cols;
+ vector<Real> elements;
+
+ Matrix(int rows, int cols):
+ rows(rows),
+ cols(cols),
+ elements(vector<Real>(rows*cols, 0)) {}
- for (unsigned int k1 = 0; k1 < F1.lowRankBlocks.size(); k1++) {
- int b1 = F1.lowRankBlocks[k1].first;
+ inline Real get(int r, int c) const {
+ return elements[r + c*rows];
+ }
+
+ inline void set(int r, int c, const Real &a) {
+ elements[r + c*rows] = a;
+ }
+
+ void setZero() {
+ std::fill(elements.begin(), elements.end(), 0);
+ }
+
+ void setIdentity() {
+ assert(rows == cols);
+
+ setZero();
+ for (int i = 0; i < rows; i++)
+ elements[i * (rows + 1)] = 1;
+ }
+
+ friend ostream& operator<<(ostream& os, const Matrix& a);
+};
+
+ostream& operator<<(ostream& os, const Matrix& a) {
+ os << "{";
+ for (int r = 0; r < a.rows; r++) {
+ os << "{";
+ for (int c = 0; c < a.cols; c++) {
+ os << a.get(r,c);
+ if (c < a.cols-1)
+ os << ", ";
+ }
+ os << "}";
+ if (r < a.rows-1)
+ os << ", ";
+ }
+ os << "}";
+ return os;
+}
+
+// result = choleskyDecomposition(a) (lower triangular)
+// Inputs:
+// - a : dim x dim symmetric matrix
+// - result : dim x dim lower-triangular matrix
+//
+void choleskyDecomposition(Matrix &a, Matrix &result) {
+ int dim = a.rows;
+ assert(a.cols == dim);
+ assert(result.rows == dim);
+ assert(result.cols == dim);
+
+ mpackint info;
+ Real *resultArray = &result.elements[0];
+
+ Rcopy(dim*dim, &a.elements[0], 1, resultArray, 1);
+
+ // The lower-triangular part of result is now our cholesky matrix
+ Rpotrf("Lower", dim, resultArray, dim, &info);
+
+ // Set the upper-triangular part of the result to zero
+ for (int j = 0; j < dim; j++) {
+ for (int i = 0; i < j; i++) {
+ result.elements[i + j*dim] = 0;
+ }
+ }
+}
+
+// result = a^-1
+// Inputs:
+// - a : dim x dim lower-triangular matrix
+// - result : dim x dim lower-triangular matrix
+//
+void inverseLowerTriangular(Matrix &a, Matrix &result) {
+ int dim = a.rows;
+ assert(a.cols == dim);
+ assert(result.rows == dim);
+ assert(result.cols == dim);
+
+ result.setIdentity();
+ Rtrsm("Left", "Lower", "NoTranspose", "NonUnitDiagonal",
+ dim, dim, 1, &a.elements[0], dim, &result.elements[0], dim);
+}
+
+// result = choleskyDecomposition(a)^-1
+// Inputs:
+// - a : dim x dim symmetric matrix
+// - work : dim x dim matrix
+// - result : dim x dim lower-triangular matrix
+//
+void inverseCholesky(Matrix &a, Matrix &work, Matrix &result) {
+ choleskyDecomposition(a, work);
+ inverseLowerTriangular(work, result);
+}
+
+// invCholesky = choleskyDecomposition(a)^-1
+// inverse = a^-1
+// Inputs:
+// - a : dim x dim symmetric matrix
+// - work : dim x dim matrix
+// - invCholesky : dim x dim lower-triangular matrix
+// - inverse : dim x dim symmetric matrix
+//
+void inverseCholeskyAndInverse(Matrix &a, Matrix &work, Matrix &invCholesky, Matrix &inverse) {
+ int dim = a.rows;
+ assert(a.cols == dim);
+ assert(work.rows == dim && work.cols == dim);
+ assert(invCholesky.rows == dim && invCholesky.cols == dim);
+ assert(inverse.rows == dim && inverse.cols == dim);
- for (unsigned int k2 = 0; k2 < F2.lowRankBlocks.size(); k2++) {
- int b2 = F2.lowRankBlocks[k2].first;
+ inverseCholesky(a, work, invCholesky);
- if (b1 == b2) {
- int i1 = F1.lowRankBlocks[k1].second.i;
- int j1 = F1.lowRankBlocks[k1].second.j;
+ inverse.elements = invCholesky.elements;
+ Rtrmm("Left", "Lower", "Transpose", "NonUnitDiag", dim, dim, 1,
+ &invCholesky.elements[0], dim,
+ &inverse.elements[0], dim);
+}
+
+// result = b'^T a b', where b' = b \otimes 1
+// Inputs:
+// - a : d x d symmetric matrix
+// - b : k x n matrix, where d = numBlocks*k, with numBlocks an integer
+// - work : d x (n*numBlocks) matrix
+// - result : (n*numBlocks) x (n*numBlocks) symmetric matrix
+//
+void blockMatrixCongruence(const Matrix &a, const Matrix &b, Matrix &work, Matrix &result) {
+ int numBlocks = a.rows / b.rows;
+
+ assert(result.rows == b.cols * numBlocks);
+ assert(result.cols == b.cols * numBlocks);
+
+ assert(work.rows == a.rows);
+ assert(work.cols == result.cols);
+
+ // work = a b'
+ for (int c = 0; c < work.cols; c++) {
+ int bCol = c % b.cols;
+ int aColOffset = (c / b.cols) * b.rows;
+
+ for (int r = 0; r < work.rows; r++) {
+
+ Real tmp = 0;
+ for (int k = 0; k < b.rows; k++) {
+ tmp += a.get(r, aColOffset + k) * b.get(k, bCol);
+ }
- int i2 = F2.lowRankBlocks[k2].second.i;
- int j2 = F2.lowRankBlocks[k2].second.j;
+ work.set(r, c, tmp);
+ }
+ }
- const vector<Real> v = F1.lowRankBlocks[k1].second.vec;
- const vector<Real> u = F2.lowRankBlocks[k2].second.vec;
+ // result = b'^T work
+ for (int c = 0; c < result.cols; c++) {
- total += (X.blocks[b1].blockBilinearForm(j2, i1, u, v) *
- ZInv.blocks[b1].blockBilinearForm(j1, i2, v, u));
+ // since result is symmetric, only compute its upper triangle
+ for (int r = 0; r <= c; r++) {
+ int bCol = r % b.cols;
+ int workRowOffset = (r / b.cols) * b.rows;
+ Real tmp = 0;
+ for (int k = 0; k < b.rows; k++) {
+ tmp += b.get(k, bCol) * work.get(workRowOffset + k, c);
+ }
+
+ result.set(r, c, tmp);
+
+ // lower triangle is the same as upper triangle
+ if (c != r) {
+ result.set(c, r, tmp);
}
}
}
- return total;
-
}
-int main(int argc, char** argv) {
+class BlockDiagonalMatrix {
+public:
+ vector<Real> diagonalPart;
+ vector<Matrix> blocks;
- Real x = 1;
- Real y = 2;
- DenseMatrix m(2);
- DenseMatrix n(2);
- m.choleskyDecomposition(n);
- vector<int> sizes;
- sizes.push_back(3);
- BlockDenseMatrix b(vector<Real>(0), sizes);
+ BlockDiagonalMatrix(const vector<Real> &diagonalPart, const vector<int> &blockSizes):
+ diagonalPart(diagonalPart) {
+ for (unsigned int i = 0; i < blockSizes.size(); i++) {
+ blocks.push_back(Matrix(blockSizes[i], blockSizes[i]));
+ }
+ }
+
+ void setZero() {
+ std::fill(diagonalPart.begin(), diagonalPart.end(), 0);
+
+ for (unsigned int b = 0; b < blocks.size(); b++) {
+ blocks[b].setZero();
+ }
+ }
+
+ void setIdentity() {
+ std::fill(diagonalPart.begin(), diagonalPart.end(), 1);
+
+ for (unsigned int b = 0; b < blocks.size(); b++) {
+ blocks[b].setIdentity();
+ }
+ }
+
+ friend ostream& operator<<(ostream& os, const BlockDiagonalMatrix& a);
+
+};
+
+ostream& operator<<(ostream& os, const BlockDiagonalMatrix& a) {
+
+ os << "{";
+ for (unsigned int i = 0; i < a.diagonalPart.size(); i++) {
+ os << a.diagonalPart[i] << ", ";
+ }
+
+ for (unsigned int b = 0; b < a.blocks.size(); b++) {
+ os << a.blocks[b];
+ if (b < a.blocks.size() - 1)
+ os << ", ";
+ }
+ os << "}";
+ return os;
+}
+
+void inverseCholeskyAndInverse(BlockDiagonalMatrix &a,
+ BlockDiagonalMatrix &work,
+ BlockDiagonalMatrix &invCholesky,
+ BlockDiagonalMatrix &inverse) {
+
+ for (unsigned int i = 0; i < a.diagonalPart.size(); i++) {
+ Real d = a.diagonalPart[i];
+ invCholesky.diagonalPart[i] = 1/sqrt(d);
+ inverse.diagonalPart[i] = 1/d;
+ }
+
+ for (unsigned int b = 0; b < a.blocks.size(); b++) {
+ inverseCholeskyAndInverse(a.blocks[b],
+ work.blocks[b],
+ invCholesky.blocks[b],
+ inverse.blocks[b]);
+ }
+}
+
+class SDPConstraint {
+public:
+ vector<vector<Real> > diagonalConstraints;
+ int row;
+ int col;
+ vector<int> blocks;
+};
+
+class SDP {
+public:
+ vector<Matrix> algebraBasisVectors;
+ vector<SDPConstraint> constraints;
+
+ int numConstraints() const {
+ int result = 0;
+ for (unsigned int i = 0; i < constraints.size(); i++) {
+ result += constraints[i].diagonalConstraints.size();
+ }
+ return result;
+ }
+};
+
+inline Real quadDotProduct(const vector<Real> &v1,
+ const vector<Real> &v2,
+ const vector<Real> &v3,
+ const vector<Real> &v4) {
+ Real result = 0;
+ for (unsigned int i = 0; i < v1.size(); i++)
+ result += v1[i]*v2[i]*v3[i]*v4[i];
+ return result;
+}
+
+void schurComplement(const SDP &sdp,
+ const BlockDiagonalMatrix &y,
+ const BlockDiagonalMatrix &xInv,
+ vector<Matrix> &bilinearPairingsY,
+ vector<Matrix> &bilinearPairingsXInv,
+ vector<Matrix> &work,
+ Matrix &schur) {
+ assert(schur.rows == sdp.numConstraints());
+ assert(schur.rows == schur.cols);
+
+ for (unsigned int b = 0; b < sdp.algebraBasisVectors.size(); b++) {
+ blockMatrixCongruence(y.blocks[b], sdp.algebraBasisVectors[b], work[b], bilinearPairingsY[b]);
+ blockMatrixCongruence(xInv.blocks[b], sdp.algebraBasisVectors[b], work[b], bilinearPairingsXInv[b]);
+ }
+
+ unsigned int r = 0;
+ for (vector<SDPConstraint>::const_iterator f1 = sdp.constraints.begin(); f1 != sdp.constraints.end(); ++f1) {
+
+ int kMax1 = f1->diagonalConstraints.size();
+ for (int k1 = 0; k1 < kMax1; k1++, r++) {
+
+ int row1 = f1->row;
+ int col1 = f1->col;
+
+ unsigned int c = 0;
+ for (vector<SDPConstraint>::const_iterator f2 = sdp.constraints.begin(); f2 != sdp.constraints.end(); ++f2) {
+
+ int kMax2 = f2->diagonalConstraints.size();
+ for (int k2 = 0; k2 < kMax2; k2++, c++) {
+
+ int row2 = f2->row;
+ int col2 = f2->col;
+
+ // TODO: take advantage of the fact that this part of the pairing is symmetric
+ Real tmp = quadDotProduct(f1->diagonalConstraints[k1], y.diagonalPart,
+ f2->diagonalConstraints[k2], xInv.diagonalPart);
+
+ for (unsigned int j1 = 0; j1 < f1->blocks.size(); j1++) {
+ int b1 = f1->blocks[j1];
+
+ for (unsigned int j2 = 0; j2 < f2->blocks.size(); j2++) {
+ int b2 = f2->blocks[j2];
+
+ if (b1 == b2) {
+ int kMax = sdp.algebraBasisVectors[b1].cols;
+ assert(kMax == kMax1);
+ assert(kMax == kMax2);
+
+ // TODO: make the division by 4 happen only once
+ tmp += (bilinearPairingsY[b1].get(row1*kMax+k1,row2*kMax+k2) * bilinearPairingsXInv[b1].get(col2*kMax+k2,col1*kMax+k1) +
+ bilinearPairingsY[b1].get(row1*kMax+k1,col2*kMax+k2) * bilinearPairingsXInv[b1].get(row2*kMax+k2,col1*kMax+k1) +
+ bilinearPairingsY[b1].get(col1*kMax+k1,col2*kMax+k2) * bilinearPairingsXInv[b1].get(row2*kMax+k2,row1*kMax+k1) +
+ bilinearPairingsY[b1].get(col1*kMax+k1,row2*kMax+k2) * bilinearPairingsXInv[b1].get(col2*kMax+k2,row1*kMax+k1))/4;
+ }
+
+ }
+ }
+
+ schur.set(r, c, tmp);
+
+ }
+ }
- cout << x+y << endl;
- for (int i = 0; i < m.dim; i++) {
- for (int j = 0; j < m.dim; j++) {
- cout << b.blocks[0].elem(i,j) << endl;
}
+
}
}
+
+void testBlockCongruence() {
+ Matrix a(4,4);
+ Matrix b(2,3);
+ Matrix result(6,6);
+ Matrix work(4,6);
+ a.setIdentity();
+ b.set(0,0,2);
+ b.set(1,0,3);
+ b.set(0,1,4);
+ b.set(1,1,5);
+ b.set(0,2,6);
+ b.set(1,2,7);
+
+ blockMatrixCongruence(a, b, work, result);
+
+ cout << a << endl;
+ cout << b << endl;
+ cout << work << endl;
+ cout << result << endl;
+
+}
+
+void testBlockDiagonalCholesky() {
+ vector<int> sizes;
+ sizes.push_back(3);
+ sizes.push_back(4);
+ vector<Real> diag(2, 0);
+
+ BlockDiagonalMatrix a(diag, sizes);
+ a.setIdentity();
+ a.diagonalPart[0] = 2;
+ a.diagonalPart[1] = 3;
+ Real aBlock0[] = {14,3,8,3,10,9,8,9,14};
+ a.blocks[0].elements.insert(a.blocks[0].elements.begin(), aBlock0, aBlock0 + 9);
+
+ BlockDiagonalMatrix work(diag, sizes);
+ BlockDiagonalMatrix invCholesky(diag, sizes);
+ BlockDiagonalMatrix inverse(diag, sizes);
+
+ inverseCholeskyAndInverse(a, work, invCholesky, inverse);
+
+ cout << a << endl;
+ cout << invCholesky << endl;
+ cout << inverse << endl;
+}
+
+int main(int argc, char** argv) {
+
+ testBlockCongruence();
+ testBlockDiagonalCholesky();
+
+}
--
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