1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
// 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
typedef CGAL::Quadratic_program<int> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;
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
CGAL::Quotient<ET> val;
for (int k=0; k<n; ++k) val += (*(a+k))[i] * inv_a[j][k];
assert (val == (i == j ? 1 : 0));
}
return 0;
}
|
