7 #include "TMatrixTSym.h"
24 int main(
int argc,
char**argv)
31 size_t numberOfInversions;
37 JParser<> zap(
"Example program to test inversion of symmetric matrix.");
47 catch(
const exception &error) {
48 FATAL(error.what() << endl);
52 gRandom->SetSeed(seed);
54 numberOfInversions += (numberOfInversions%2 == 0 ? 1 : 0);
56 const Double_t
xmin = -1.0;
57 const Double_t
xmax = +1.0;
61 for (
int i = 0; i != N; ++i) {
63 A(i,i) = gRandom->Uniform(
xmin,
xmax);
65 for (
int j = 0;
j != i; ++
j) {
70 DEBUG(
"Matrix A" << endl);
77 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
79 for (
size_t i = numberOfInversions; i != 0; --i) {
83 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
85 NOTICE(
"Elapsed time (Jpp) [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl);
86 DEBUG(
"Matrix A^-1" << endl);
93 DEBUG(
"Matrix A x A^-1" << endl);
96 NOTICE(
"A x A^-1 = I ? " <<
C.isIdentity(precision) << endl);
98 if (!
C.isIdentity(precision)) {
99 ERROR(
"Matrix A x A^-1 /= I" << endl);
102 ASSERT(
C.isIdentity(precision));
105 FATAL(error << endl);
110 TMatrixTSym<double> B(A.
size());
112 for (
size_t row = 0; row != A.
size(); ++row) {
113 for (
size_t col = 0; col != A.
size(); ++col) {
114 B(row, col) = A(row, col);
118 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
120 for (
size_t i = numberOfInversions; i != 0; --i) {
124 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
126 NOTICE(
"Elapsed time (ROOT) [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl);
129 FATAL(error << endl);
140 for (
int i = 0; i != N; ++i) {
146 for (
int i = 0; i != N; ++i) {
147 for (
int j = 0;
j != N; ++
j) {
148 Y[i] += B(i,
j) * X[
j];
156 for (
int i = 0; i != N; ++i) {
158 DEBUG(setw(4) << i <<
' ' <<
FIXED(9,5) << X[i] <<
' ' <<
FIXED(9,5) << Y[i] << endl);
160 ASSERT(fabs(X[i] - Y[i]) <= precision);
164 FATAL(error << endl);
173 const double big = 1.0e5;
175 for (
int i = 0; i != N; ++i) {
179 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
183 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
185 NOTICE(
"Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl);
186 NOTICE(
"Pivot " << setw(2) << i <<
'+' << big <<
' ' << (
JMatrixNS(A).mul(B).isIdentity(precision) ?
"okay" :
"error") << endl);
196 FATAL(error << endl);
202 for (
int i = 0; i != N; ++i) {
203 for (
int j = 0;
j != N; ++
j) {
204 ASSERT(A(i,
j) == 0,
"Test reset (" << i <<
"," <<
j <<
")");
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
#define ASSERT(A,...)
Assert macro.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Utility class to parse command line options.
static const double C
Physics constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for floating point format specification.
size_t size() const
Get dimension of matrix.
JMatrixND & reset()
Set matrix to the null matrix.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
void invert()
Invert matrix according LDU decomposition.