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) {
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));
104 catch (
const JException& error) {
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);
128 catch (
const JException& error) {
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);
163 catch (
const JException& error) {
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);
188 ASSERT(JMatrixNS(A).mul(B).isIdentity(precision));
195 catch (
const JException& error) {
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 <<
")");
Utility class to parse command line options.
int main(int argc, char *argv[])
then usage $script< input file >[option] nPossible options count
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
static const double C
Physics constants.
#define ASSERT(A,...)
Assert macro.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
General purpose messaging.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Utility class to parse command line options.
no fit printf nominal n $STRING awk v X
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
#define DEBUG(A)
Message macros.