28 void reportTime( std::string s,
double time);
32 int pr = std::cout.precision(8);
33 std::cout << s <<
"\t" <<
" time = " << time.
RealTime() <<
"\t(sec)\t" 36 std::cout.precision(pr);
45 Timer(
const std::string & s =
"") : fName(s), fTime(0)
49 Timer(
double & t,
const std::string & s =
"") : fName(s), fTime(&t)
59 reportTime(fName, fWatch.RealTime() );
61 if (fTime) *fTime += fWatch.RealTime();
94 std::cout << m << std::endl;
106 for (
int i = 0; i <
N; ++i) {
107 double maxVal = i*10000/(N-1) + 1;
110 for (
int i = 0; i <
N; ++i) {
111 for (
int j = 0; j < i; ++j) {
124 for (
int i = 0; i <
n; ++i) {
125 v[i] = createMatrix<M>();
137 template <
class M,
class Type>
138 struct TestInverter {
139 static bool Inv (
const M & , M & ) {
return false;}
140 static bool Inv2 ( M & ) {
return false;}
144 struct TestInverter<SymMatrix, Choleski> {
145 static bool Inv (
const SymMatrix &
m, SymMatrix &
result ) {
150 static bool Inv2 ( SymMatrix & m ) {
156 struct TestInverter<SymMatrix, BK> {
157 static bool Inv (
const SymMatrix &
m, SymMatrix &
result ) {
162 static bool Inv2 ( SymMatrix & m ) {
168 struct TestInverter<SymMatrix, Cramer> {
169 static bool Inv (
const SymMatrix &
m, SymMatrix &
result ) {
174 static bool Inv2 ( SymMatrix & m ) {
181 struct TestInverter<SymMatrix, QR> {
182 static bool Inv (
const SymMatrix &
m, SymMatrix &
result ) {
183 ROOT::Math::QRDecomposition<double> d;
194 struct TestInverter<TMatrixDSym, Default> {
195 static bool Inv (
const TMatrixDSym &
m, TMatrixDSym &
result ) {
200 static bool Inv2 ( TMatrixDSym & m ) {
207 struct TestInverter<TMatrixDSym, Cramer> {
208 static bool Inv (
const TMatrixDSym &
m, TMatrixDSym &
result ) {
213 static bool Inv2 ( TMatrixDSym & m ) {
220 struct TestInverter<TMatrixDSym, Choleski> {
221 static bool Inv (
const TMatrixDSym &
m, TMatrixDSym &
result ) {
230 struct TestInverter<TMatrixDSym, BK> {
231 static bool Inv (
const TMatrixDSym &
m, TMatrixDSym &
result ) {
240 template<
class M,
class T>
241 double invert(
const std::vector<M* > & matlist,
double & time,std::string s) {
242 M
result = *(matlist.front());
244 int nloop = matlist.size();
246 for (
int l = 0;
l < nloop;
l++)
248 const M &
m = *(matlist[
l]);
249 bool ok = TestInverter<M,T>::Inv(m,result);
251 std::cout <<
"inv failed for matrix " <<
l << std::endl;
262 template<
class M,
class T>
263 double invert2(
const std::vector<M* > & matlist,
double & time,std::string s) {
266 int nloop = matlist.size();
267 std::vector<M *> vmat(nloop);
268 for (
int l = 0;
l < nloop;
l++)
270 vmat[
l] =
new M( *matlist[
l] );
275 for (
int l = 0;
l < nloop;
l++)
278 bool ok = TestInverter<M,T>::Inv2(m);
280 std::cout <<
"inv failed for matrix " <<
l << std::endl;
289 bool equal(
double d1,
double d2,
double stol = 10000) {
290 std::cout.precision(12);
292 if ( std::abs(d1) > 0 && std::abs(d2) > 0 )
293 return ( std::abs( d1-d2) < eps * std::max(std::abs(d1), std::abs(d2) ) );
295 return std::abs(d2) < eps;
297 return std::abs(d1) < eps;
305 std::vector<SymMatrix *> v1(n);
307 std::vector<TMatrixDSym *> v2(n);
313 double s1 = invert<SymMatrix, Choleski> (v1, time,
"SMatrix Chol");
314 double s2 = invert<SymMatrix, BK> (v1, time,
"SMatrix BK");
315 double s3 = invert<SymMatrix, Cramer> (v1, time,
"SMatrix Cram");
318 std::cout <<
"result SMatrix choleski " << s1 <<
" BK " << s2 <<
" cramer " << s3 << std::endl;
319 std::cerr <<
"Error: inversion test for SMatrix FAILED ! " << std::endl;
322 std::cout << std::endl;
324 double m1 = invert<TMatrixDSym, Choleski> (v2, time,
"TMatrix Chol");
325 double m2 = invert<TMatrixDSym, BK> (v2, time,
"TMatrix BK");
326 double m3 = invert<TMatrixDSym, Cramer> (v2, time,
"TMatrix Cram");
327 double m4 = invert<TMatrixDSym, Default> (v2, time,
"TMatrix Def");
331 std::cout <<
"result TMatrix choleski " << m1 <<
" BK " << m2
332 <<
" cramer " << m3 <<
" default " << m4 << std::endl;
333 std::cerr <<
"Error: inversion test for TMatrix FAILED ! " << std::endl;
336 std::cout << std::endl;
341 std::cout <<
"\n - self inversion test \n";
342 double s11 = invert2<SymMatrix, Choleski> (v1, time,
"SMatrix Chol");
343 double s12 = invert2<SymMatrix, BK> (v1, time,
"SMatrix BK");
344 double s13 = invert2<SymMatrix, Cramer> (v1, time,
"SMatrix Cram");
347 std::cout <<
"result SMatrix choleski " << s11 <<
" BK " << s12 <<
" cramer " << s13 << std::endl;
348 std::cerr <<
"Error: self inversion test for SMatrix FAILED ! " << std::endl;
351 std::cout << std::endl;
353 double m13 = invert2<TMatrixDSym, Cramer> (v2, time,
"TMatrix Cram");
354 double m14 = invert2<TMatrixDSym, Default> (v2, time,
"TMatrix Def");
355 ok = (
equal(m13,m14) );
357 std::cout <<
"result TMatrix cramer " << m13 <<
" default " << m14 << std::endl;
358 std::cerr <<
"Error: self inversion test for TMatrix FAILED ! " << std::endl;
361 std::cout << std::endl;
368 std::cout <<
"Test Inversion for matrix with N = " <<
N << std::endl;
370 std::cerr <<
"Test inversion of positive defined matrix ....... ";
371 if (ok) std::cerr <<
"OK \n";
372 else std::cerr <<
"FAILED \n";
373 return (ok) ? 0 : -1;
SMatrix< T, D1, D2, R > InverseChol(int &ifail) const
Invert of a symmetric positive defined Matrix using Choleski decomposition.
static long int sum(long int i)
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
virtual void Print(Option_t *option="") const
Dump this marker with its attributes.
bool equal(double d1, double d2, double stol=10000)
SMatrix< double, N, N, MatRepSym< double, N > > SymMatrix
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
bool InvertFast()
Fast Invertion of a square Matrix ( this method changes the current matrix).
bool InvertChol()
Invertion of a symmetric positive defined Matrix using Choleski decomposition.
SMatrix< T, D1, D2, R > InverseFast(int &ifail) const
Invert a square Matrix and returns a new matrix.
double invert(const std::vector< M * > &matlist, double &time, std::string s)
The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using.
bool stressSymPosInversion(int n, bool selftest)
SMatrix: a generic fixed size D1 x D2 Matrix class.
This is the base class for the ROOT Random number generators.
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
double invert2(const std::vector< M * > &matlist, double &time, std::string s)
Cholesky Decomposition class.
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
RooCmdArg Timer(Bool_t flag=kTRUE)
void printMatrix(const M &m)
R__EXTERN TRandom * gRandom
SMatrix< T, D1, D2, R > Inverse(int &ifail) const
Invert a square Matrix and returns a new matrix.
void printTime(TStopwatch &time, std::string s)
TMatrixTSym< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant.
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
void generate(std::vector< M *> &v)
int testInversion(int n=100000)
virtual Bool_t Decompose()
Matrix A is decomposed in components U and D so that A = U*D*U^T If the decomposition succeeds...