\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)
CGAL 4.12 - Linear and Quadratic Programming Solver
QP_solver/invert_matrix.cpp
// example: invert a random matrix (or find out that it's singular)
#include <iostream>
#include <vector>
#include <algorithm>
#include <CGAL/basic.h>
#include <CGAL/Random.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
// choose exact integral type
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpz.h>
typedef CGAL::Gmpz ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif
int n = 4; // dimension of matrix
int s = 10; // coordinates are in {-s, -s+1,...,0,...,s-1,s}
CGAL::Random rd; // random number generator
// program and solution types
int main() {
std::vector<std::vector<CGAL::Quotient<ET> > >
inv_a; // stored columnwise
// we need a free LP (no variable bounds) with Ax = b
Program lp (CGAL::EQUAL, false, 0, false, 0);
// constraint matrix A: the random matrix to be inverted
for (int j=0; j<n; ++j)
for (int i=0; i<n; ++i)
lp.set_a (j, i, rd.get_int (-s, s));
// we need to solve n LP, one for every right-hand side b = e_j to
// get j-th column of the inverse
for (int j=0; j<n; ++j) {
lp.set_b (j, 1);
// solve the lp, using ET as the exact type
Solution s = CGAL::solve_linear_program(lp, ET());
if (s.is_infeasible()) {
std::cout << "matrix is singular" << std::endl;
return 0;
} else {
// store solution
inv_a.push_back(std::vector<CGAL::Quotient<ET> >());
std::copy (s.variable_values_begin(), s.variable_values_end(),
std::back_inserter (inv_a[j]));
}
lp.set_b (j, 0); // reset for next round
}
// output both matrices, and check that they are indeed inverse to each
// other
std::cout << "Random matrix A...:" << std::endl;
Program::A_iterator a = lp.get_a();
for (int i=0; i<n; ++i) {
for (int j=0; j<n; ++j)
std::cout << (*(a+j))[i] << " "; // row i
std::cout << std::endl;
}
std::cout << std::endl << "...and its inverse: " << std::endl;
for (int i=0; i<n; ++i) {
for (int j=0; j<n; ++j)
std::cout << inv_a[j][i] << " "; // row i
std::cout << std::endl;
}
// check inverse property
for (int i=0; i<n; ++i)
for (int j=0; j<n; ++j) {
// i-th row of A times j-th column of inverse
for (int k=0; k<n; ++k) val += (*(a+k))[i] * inv_a[j][k];
assert (val == (i == j ? 1 : 0));
}
return 0;
}