00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 #include <stdlib.h>
00092 #include <Riostream.h>
00093 #include <TSystem.h>
00094 #include <TFile.h>
00095 #include <TBenchmark.h>
00096 #include <TArrayD.h>
00097 #include <TF1.h>
00098 #include <TGraph.h>
00099 #include <TROOT.h>
00100 #include "TMath.h"
00101
00102 #include "TMatrixF.h"
00103 #include "TMatrixFSym.h"
00104 #include "TMatrixFSparse.h"
00105 #include "TMatrixFLazy.h"
00106 #include "TVectorF.h"
00107
00108 #include "TMatrixD.h"
00109 #include "TMatrixDSym.h"
00110 #include "TMatrixDSparse.h"
00111 #include "TMatrixDLazy.h"
00112 #include "TMatrixDUtils.h"
00113 #include "TVectorD.h"
00114
00115 #include "TDecompLU.h"
00116 #include "TDecompChol.h"
00117 #include "TDecompQRH.h"
00118 #include "TDecompSVD.h"
00119 #include "TDecompBK.h"
00120 #include "TMatrixDEigen.h"
00121 #include "TMatrixDSymEigen.h"
00122
00123 void stressLinear (Int_t maxSizeReq=100,Int_t verbose=0);
00124 void StatusPrint (Int_t id,const TString &title,Bool_t status);
00125
00126 void mstress_allocation (Int_t msize);
00127 void mstress_matrix_fill (Int_t rsize,Int_t csize);
00128 void mstress_element_op (Int_t rsize,Int_t csize);
00129 void mstress_binary_ebe_op (Int_t rsize, Int_t csize);
00130 void mstress_transposition (Int_t msize);
00131 void mstress_special_creation (Int_t dim);
00132 void mstress_matrix_fill (Int_t rsize,Int_t csize);
00133 void mstress_matrix_promises (Int_t dim);
00134 void mstress_norms (Int_t rsize,Int_t csize);
00135 void mstress_determinant (Int_t msize);
00136 void mstress_mm_multiplications ();
00137 void mstress_sym_mm_multiplications(Int_t msize);
00138 void mstress_vm_multiplications ();
00139 void mstress_inversion ();
00140 void mstress_matrix_io ();
00141
00142 void spstress_allocation (Int_t msize);
00143 void spstress_matrix_fill (Int_t rsize,Int_t csize);
00144 void spstress_element_op (Int_t rsize,Int_t csize);
00145 void spstress_binary_ebe_op (Int_t rsize, Int_t csize);
00146 void spstress_transposition (Int_t msize);
00147 void spstress_norms (Int_t rsize,Int_t csize);
00148 void spstress_mm_multiplications ();
00149 void spstress_vm_multiplications ();
00150 void spstress_matrix_slices (Int_t vsize);
00151 void spstress_matrix_io ();
00152
00153 void vstress_allocation (Int_t msize);
00154 void vstress_element_op (Int_t vsize);
00155 void vstress_binary_op (Int_t vsize);
00156 void vstress_norms (Int_t vsize);
00157 void vstress_matrix_slices (Int_t vsize);
00158 void vstress_vector_io ();
00159
00160 Bool_t test_svd_expansion (const TMatrixD &A);
00161 void astress_decomp ();
00162 void astress_lineqn ();
00163 void astress_pseudo ();
00164 void astress_eigen (Int_t msize);
00165 void astress_decomp_io (Int_t msize);
00166
00167 void stress_backward_io ();
00168
00169 void cleanup ();
00170
00171
00172 #ifndef __CINT__
00173
00174 int main(int argc,const char *argv[])
00175 {
00176 Int_t maxSizeReq = 100;
00177 if (argc > 1) maxSizeReq = atoi(argv[1]);
00178 Int_t verbose = 0;
00179 if (argc > 2) verbose = (atoi(argv[2]) > 0 ? 1 : 0);
00180
00181 gBenchmark = new TBenchmark();
00182 stressLinear(maxSizeReq,verbose);
00183
00184 cleanup();
00185 return 0;
00186 }
00187
00188 #endif
00189
00190 #define EPSILON 1.0e-14
00191
00192 Int_t gVerbose = 0;
00193 Int_t gNrLoop;
00194 const Int_t nrSize = 20;
00195 const Int_t gSizeA[] = {5,6,7,8,9,10,15,20,30,40,50,60,70,80,90,100,300,500,700,1000};
00196
00197
00198
00199 void stressLinear(Int_t maxSizeReq,Int_t verbose)
00200 {
00201 cout << "******************************************************************" <<endl;
00202 cout << "* Starting Linear Algebra - S T R E S S suite *" <<endl;
00203 cout << "******************************************************************" <<endl;
00204 cout << "******************************************************************" <<endl;
00205
00206 gVerbose = verbose;
00207
00208 gBenchmark->Start("stressLinear");
00209
00210 gNrLoop = nrSize-1;
00211 while (gNrLoop > 0 && maxSizeReq < gSizeA[gNrLoop])
00212 gNrLoop--;
00213
00214 const Int_t maxSize = gSizeA[gNrLoop];
00215
00216
00217 {
00218 cout << "* Starting Matrix - S T R E S S *" <<endl;
00219 cout << "******************************************************************" <<endl;
00220 mstress_allocation(maxSize);
00221 mstress_matrix_fill(maxSize,maxSize/2);
00222 mstress_element_op(maxSize,maxSize/2);
00223 mstress_binary_ebe_op(maxSize/2,maxSize);
00224 mstress_transposition(maxSize);
00225 mstress_special_creation(maxSize);
00226 #ifndef __CINT__
00227 mstress_matrix_promises(maxSize);
00228 #endif
00229 mstress_norms(maxSize,maxSize/2);
00230 mstress_determinant(maxSize);
00231 mstress_mm_multiplications();
00232 mstress_sym_mm_multiplications(maxSize);
00233 mstress_vm_multiplications();
00234 mstress_inversion();
00235
00236 mstress_matrix_io();
00237 cout << "******************************************************************" <<endl;
00238 }
00239
00240
00241 {
00242 cout << "* Starting Sparse Matrix - S T R E S S *" <<endl;
00243 cout << "******************************************************************" <<endl;
00244 spstress_allocation(maxSize);
00245 spstress_matrix_fill(maxSize,maxSize/2);
00246 spstress_element_op(maxSize,maxSize/2);
00247 spstress_binary_ebe_op(maxSize/2,maxSize);
00248 spstress_transposition(maxSize);
00249 spstress_norms(maxSize,maxSize/2);
00250 spstress_mm_multiplications();
00251 spstress_vm_multiplications();
00252 spstress_matrix_slices(maxSize);
00253 spstress_matrix_io();
00254 cout << "******************************************************************" <<endl;
00255 }
00256
00257 {
00258 cout << "* Starting Vector - S T R E S S *" <<endl;
00259 cout << "******************************************************************" <<endl;
00260 vstress_allocation(maxSize);
00261 vstress_element_op(maxSize);
00262 vstress_binary_op(maxSize);
00263 vstress_norms(maxSize);
00264 vstress_matrix_slices(maxSize);
00265 vstress_vector_io();
00266 cout << "******************************************************************" <<endl;
00267 }
00268
00269
00270 {
00271 cout << "* Starting Linear Algebra - S T R E S S *" <<endl;
00272 cout << "******************************************************************" <<endl;
00273 astress_decomp();
00274 astress_lineqn();
00275 astress_pseudo();
00276 astress_eigen(5);
00277 astress_decomp_io(10);
00278 cout << "******************************************************************" <<endl;
00279 }
00280
00281
00282 {
00283 cout << "* Starting Backward IO compatibility - S T R E S S *" <<endl;
00284 cout << "******************************************************************" <<endl;
00285 stress_backward_io();
00286 cout << "******************************************************************" <<endl;
00287 }
00288
00289 gBenchmark->Stop("stressLinear");
00290
00291
00292 Bool_t UNIX = strcmp(gSystem->GetName(), "Unix") == 0;
00293 printf("******************************************************************\n");
00294 if (UNIX) {
00295 TString sp = gSystem->GetFromPipe("uname -a");
00296 sp.Resize(60);
00297 printf("* SYS: %s\n",sp.Data());
00298 if (strstr(gSystem->GetBuildNode(),"Linux")) {
00299 sp = gSystem->GetFromPipe("lsb_release -d -s");
00300 printf("* SYS: %s\n",sp.Data());
00301 }
00302 if (strstr(gSystem->GetBuildNode(),"Darwin")) {
00303 sp = gSystem->GetFromPipe("sw_vers -productVersion");
00304 sp += " Mac OS X ";
00305 printf("* SYS: %s\n",sp.Data());
00306 }
00307 } else {
00308 const Char_t *os = gSystem->Getenv("OS");
00309 if (!os) printf("* SYS: Windows 95\n");
00310 else printf("* SYS: %s %s \n",os,gSystem->Getenv("PROCESSOR_IDENTIFIER"));
00311 }
00312
00313 printf("******************************************************************\n");
00314 gBenchmark->Print("stressLinear");
00315 const Int_t nr = 7;
00316 const Double_t x_b12[] = { 10., 30., 50., 100., 300., 500., 700.};
00317 const Double_t y_b12[] = {10.74, 15.72, 20.00, 35.79, 98.77, 415.34, 1390.33};
00318
00319 TGraph gr(nr,x_b12,y_b12);
00320 Double_t ct = gBenchmark->GetCpuTime("stressLinear");
00321 const Double_t rootmarks = 600*gr.Eval(maxSize)/ct;
00322
00323 printf("******************************************************************\n");
00324 printf("* ROOTMARKS =%6.1f * Root%-8s %d/%d\n",rootmarks,gROOT->GetVersion(),
00325 gROOT->GetVersionDate(),gROOT->GetVersionTime());
00326 printf("******************************************************************\n");
00327 }
00328
00329
00330 void StatusPrint(Int_t id,const TString &title,Bool_t status)
00331 {
00332
00333 const Int_t kMAX = 65;
00334 Char_t number[4];
00335 sprintf(number,"%2d",id);
00336 TString header = TString("Test ")+number+" : "+title;
00337 const Int_t nch = header.Length();
00338 for (Int_t i = nch; i < kMAX; i++) header += '.';
00339 cout << header << (status ? "OK" : "FAILED") << endl;
00340 }
00341
00342
00343
00344
00345 void mstress_allocation(Int_t msize)
00346 {
00347 if (gVerbose)
00348 cout << "\n\n---> Test allocation and compatibility check" << endl;
00349
00350 Int_t i,j;
00351 Bool_t ok = kTRUE;
00352
00353 TMatrixD m1(4,msize);
00354 for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
00355 for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
00356 m1(i,j) = TMath::Pi()*i+TMath::E()*j;
00357
00358 TMatrixD m2(0,3,0,msize-1);
00359 TMatrixD m3(1,4,0,msize-1);
00360 TMatrixD m4(m1);
00361
00362 if (gVerbose) {
00363 cout << "\nStatus information reported for matrix m3:" << endl;
00364 cout << " Row lower bound ... " << m3.GetRowLwb() << endl;
00365 cout << " Row upper bound ... " << m3.GetRowUpb() << endl;
00366 cout << " Col lower bound ... " << m3.GetColLwb() << endl;
00367 cout << " Col upper bound ... " << m3.GetColUpb() << endl;
00368 cout << " No. rows ..........." << m3.GetNrows() << endl;
00369 cout << " No. cols ..........." << m3.GetNcols() << endl;
00370 cout << " No. of elements ...." << m3.GetNoElements() << endl;
00371 }
00372
00373 if (gVerbose)
00374 cout << "\nCheck matrices 1 & 2 for compatibility" << endl;
00375 ok &= AreCompatible(m1,m2,gVerbose);
00376
00377 if (gVerbose)
00378 cout << "Check matrices 1 & 4 for compatibility" << endl;
00379 ok &= AreCompatible(m1,m4,gVerbose);
00380
00381 if (gVerbose)
00382 cout << "m2 has to be compatible with m3 after resizing to m3" << endl;
00383 m2.ResizeTo(m3);
00384 ok &= AreCompatible(m2,m3,gVerbose);
00385
00386 TMatrixD m5(m1.GetNrows()+1,m1.GetNcols()+5);
00387 for (i = m5.GetRowLwb(); i <= m5.GetRowUpb(); i++)
00388 for (j = m5.GetColLwb(); j <= m5.GetColUpb(); j++)
00389 m5(i,j) = TMath::Pi()*i+TMath::E()*j;
00390
00391 if (gVerbose)
00392 cout << "m1 has to be compatible with m5 after resizing to m5" << endl;
00393 m1.ResizeTo(m5.GetNrows(),m5.GetNcols());
00394 ok &= AreCompatible(m1,m5,gVerbose);
00395
00396 if (gVerbose)
00397 cout << "m1 has to be equal to m4 after stretching and shrinking" << endl;
00398 m1.ResizeTo(m4.GetNrows(),m4.GetNcols());
00399 ok &= VerifyMatrixIdentity(m1,m4,gVerbose,EPSILON);
00400 if (gVerbose)
00401 cout << "m5 has to be equal to m1 after shrinking" << endl;
00402 m5.ResizeTo(m1.GetNrows(),m1.GetNcols());
00403 ok &= VerifyMatrixIdentity(m1,m5,gVerbose,EPSILON);
00404
00405 if (gVerbose)
00406 cout << "stretching and shrinking for small matrices (stack)" << endl;
00407 if (gVerbose)
00408 cout << "m8 has to be equal to m7 after stretching and shrinking" << endl;
00409 TMatrixD m6(4,4);
00410 for (i = m6.GetRowLwb(); i <= m6.GetRowUpb(); i++)
00411 for (j = m6.GetColLwb(); j <= m6.GetColUpb(); j++)
00412 m6(i,j) = TMath::Pi()*i+TMath::E()*j;
00413 TMatrixD m8(3,3);
00414 for (i = m8.GetRowLwb(); i <= m8.GetRowUpb(); i++)
00415 for (j = m8.GetColLwb(); j <= m8.GetColUpb(); j++)
00416 m8(i,j) = TMath::Pi()*i+TMath::E()*j;
00417 TMatrixD m7(m8);
00418
00419 m8.ResizeTo(4,4);
00420 m8.ResizeTo(3,3);
00421 ok &= VerifyMatrixIdentity(m7,m8,gVerbose,EPSILON);
00422
00423 if (gVerbose)
00424 cout << "m6 has to be equal to m8 after shrinking" << endl;
00425 m6.ResizeTo(3,3);
00426 ok &= VerifyMatrixIdentity(m6,m8,gVerbose,EPSILON);
00427
00428 if (gVerbose)
00429 cout << "\nDone\n" << endl;
00430
00431 StatusPrint(1,"Allocation, Resizing",ok);
00432 }
00433
00434 class FillMatrix : public TElementPosActionD {
00435 Int_t no_elems,no_cols;
00436 void Operation(Double_t &element) const
00437 { element = 4*TMath::Pi()/no_elems * (fI*no_cols+fJ); }
00438 public:
00439 FillMatrix() {}
00440 FillMatrix(const TMatrixD &m) :
00441 no_elems(m.GetNoElements()),no_cols(m.GetNcols()) { }
00442 };
00443
00444
00445
00446
00447
00448 void mstress_matrix_fill(Int_t rsize,Int_t csize)
00449 {
00450 if (gVerbose)
00451 cout << "\n\n---> Test different matrix filling methods\n" << endl;
00452
00453 Bool_t ok = kTRUE;
00454 if (gVerbose)
00455 cout << "Creating m with Apply function..." << endl;
00456 TMatrixD m(-1,rsize-2,1,csize);
00457 #ifndef __CINT__
00458 FillMatrix f(m);
00459 m.Apply(f);
00460 #else
00461 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00462 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
00463 m(i,j) = 4*TMath::Pi()/m.GetNoElements() * (i*m.GetNcols()+j);
00464 #endif
00465
00466 {
00467 if (gVerbose)
00468 cout << "Check identity between m and matrix filled through (i,j)" << endl;
00469
00470 TMatrixD m_overload1(-1,rsize-2,1,csize);
00471 TMatrixD m_overload2(-1,rsize-2,1,csize);
00472
00473 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00474 {
00475 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
00476 {
00477 const Double_t val = 4*TMath::Pi()/rsize/csize*(i*csize+j);
00478 m_overload1(i,j) = val;
00479 m_overload2[i][j] = val;
00480 }
00481 }
00482
00483 ok &= VerifyMatrixIdentity(m,m_overload1,gVerbose,EPSILON);
00484 if (gVerbose)
00485 cout << "Check identity between m and matrix filled through [i][j]" << endl;
00486 ok &= VerifyMatrixIdentity(m,m_overload2,gVerbose,EPSILON);
00487 if (gVerbose)
00488 cout << "Check identity between matrix filled through [i][j] and (i,j)" << endl;
00489 ok &= VerifyMatrixIdentity(m_overload1,m_overload2,gVerbose,EPSILON);
00490 }
00491
00492 {
00493 TArrayD a_fortran(rsize*csize);
00494 TArrayD a_c (rsize*csize);
00495 for (Int_t i = 0; i < rsize; i++)
00496 {
00497 for (Int_t j = 0; j < csize; j++)
00498 {
00499 a_c[i*csize+j] = 4*TMath::Pi()/rsize/csize*((i-1)*csize+j+1);
00500 a_fortran[i+rsize*j] = a_c[i*csize+j];
00501 }
00502 }
00503
00504 if (gVerbose)
00505 cout << "Creating m_fortran by filling with fortran stored matrix" << endl;
00506 TMatrixD m_fortran(-1,rsize-2,1,csize,a_fortran.GetArray(),"F");
00507 if (gVerbose)
00508 cout << "Check identity between m and m_fortran" << endl;
00509 ok &= VerifyMatrixIdentity(m,m_fortran,gVerbose,EPSILON);
00510
00511 if (gVerbose)
00512 cout << "Creating m_c by filling with c stored matrix" << endl;
00513 TMatrixD m_c(-1,rsize-2,1,csize,a_c.GetArray());
00514 if (gVerbose)
00515 cout << "Check identity between m and m_c" << endl;
00516 ok &= VerifyMatrixIdentity(m,m_c,gVerbose,EPSILON);
00517 }
00518
00519 {
00520 if (gVerbose)
00521 cout << "Check insertion/extraction of sub-matrices" << endl;
00522 {
00523 TMatrixD m_sub1 = m;
00524 m_sub1.ResizeTo(0,rsize-2,2,csize);
00525 TMatrixD m_sub2 = m.GetSub(0,rsize-2,2,csize,"");
00526 ok &= VerifyMatrixIdentity(m_sub1,m_sub2,gVerbose,EPSILON);
00527 }
00528
00529 {
00530 TMatrixD m2(-1,rsize-2,1,csize);
00531 TMatrixD m_part1 = m.GetSub(0,rsize-2,2,csize,"");
00532 TMatrixD m_part2 = m.GetSub(0,rsize-2,1,1,"");
00533 TMatrixD m_part3 = m.GetSub(-1,-1,2,csize,"");
00534 TMatrixD m_part4 = m.GetSub(-1,-1,1,1,"");
00535 m2.SetSub(0,2,m_part1);
00536 m2.SetSub(0,1,m_part2);
00537 m2.SetSub(-1,2,m_part3);
00538 m2.SetSub(-1,1,m_part4);
00539 ok &= VerifyMatrixIdentity(m,m2,gVerbose,EPSILON);
00540 }
00541
00542 {
00543 TMatrixD m2(-1,rsize-2,1,csize);
00544 TMatrixD m_part1 = m.GetSub(0,rsize-2,2,csize,"S");
00545 TMatrixD m_part2 = m.GetSub(0,rsize-2,1,1,"S");
00546 TMatrixD m_part3 = m.GetSub(-1,-1,2,csize,"S");
00547 TMatrixD m_part4 = m.GetSub(-1,-1,1,1,"S");
00548 m2.SetSub(0,2,m_part1);
00549 m2.SetSub(0,1,m_part2);
00550 m2.SetSub(-1,2,m_part3);
00551 m2.SetSub(-1,1,m_part4);
00552 ok &= VerifyMatrixIdentity(m,m2,gVerbose,EPSILON);
00553 }
00554 }
00555
00556 {
00557 if (gVerbose)
00558 cout << "Check sub-matrix views" << endl;
00559 {
00560 TMatrixD m3(-1,rsize-2,1,csize);
00561 TMatrixDSub(m3,0,rsize-2,2,csize) = TMatrixDSub(m,0,rsize-2,2,csize);
00562 TMatrixDSub(m3,0,rsize-2,1,1) = TMatrixDSub(m,0,rsize-2,1,1);
00563 TMatrixDSub(m3,-1,-1,2,csize) = TMatrixDSub(m,-1,-1,2,csize);
00564 TMatrixDSub(m3,-1,-1,1,1) = TMatrixDSub(m,-1,-1,1,1);
00565 ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
00566
00567 TMatrixD unit(3,3);
00568 TMatrixDSub(m3,1,3,1,3) = unit.UnitMatrix();
00569 TMatrixDSub(m3,1,3,1,3) *= m.GetSub(1,3,1,3);
00570 ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
00571
00572 TMatrixDSub(m3,0,rsize-2,2,csize) = 1.0;
00573 TMatrixDSub(m3,0,rsize-2,1,1) = 1.0;
00574 TMatrixDSub(m3,-1,-1,2,csize) = 1.0;
00575 TMatrixDSub(m3,-1,-1,1,1) = 1.0;
00576 ok &= (m3 == 1.0);
00577
00578 }
00579 }
00580
00581 {
00582 if (gVerbose)
00583 cout << "Check array Use" << endl;
00584 {
00585 TMatrixD *m1a = new TMatrixD(m);
00586 TMatrixD *m2a = new TMatrixD();
00587 m2a->Use(m1a->GetRowLwb(),m1a->GetRowUpb(),m1a->GetColLwb(),m1a->GetColUpb(),m1a->GetMatrixArray());
00588 ok &= VerifyMatrixIdentity(m,*m2a,gVerbose,EPSILON);
00589 m2a->Sqr();
00590 TMatrixD m4 = m; m4.Sqr();
00591 ok &= VerifyMatrixIdentity(m4,*m1a,gVerbose,EPSILON);
00592 delete m1a;
00593 delete m2a;
00594 }
00595 }
00596
00597 if (gVerbose)
00598 cout << "\nDone\n" << endl;
00599
00600 StatusPrint(2,"Filling, Inserting, Using",ok);
00601 }
00602
00603
00604
00605
00606
00607 typedef Double_t (*dfunc)(Double_t);
00608 class ApplyFunction : public TElementActionD {
00609 dfunc fFunc;
00610 void Operation(Double_t &element) const { element = fFunc(Double_t(element)); }
00611 public:
00612 ApplyFunction(dfunc func) : fFunc(func) { }
00613 };
00614
00615 void mstress_element_op(Int_t rsize,Int_t csize)
00616 {
00617 Bool_t ok = kTRUE;
00618 const Double_t pattern = 8.625;
00619
00620 TMatrixD m(-1,rsize-2,1,csize);
00621
00622 if (gVerbose)
00623 cout << "\nWriting zeros to m..." << endl;
00624 {
00625 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00626 for(Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
00627 m(i,j) = 0;
00628 ok &= VerifyMatrixValue(m,0.,gVerbose,EPSILON);
00629 }
00630
00631 if (gVerbose)
00632 cout << "Creating zero m1 ..." << endl;
00633 TMatrixD m1(TMatrixD::kZero, m);
00634 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
00635
00636 if (gVerbose)
00637 cout << "Comparing m1 with 0 ..." << endl;
00638 R__ASSERT(m1 == 0);
00639 R__ASSERT(!(m1 != 0));
00640
00641 if (gVerbose)
00642 cout << "Writing a pattern " << pattern << " by assigning to m(i,j)..." << endl;
00643 {
00644 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00645 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
00646 m(i,j) = pattern;
00647 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
00648 }
00649
00650 if (gVerbose)
00651 cout << "Writing the pattern by assigning to m1 as a whole ..." << endl;
00652 m1 = pattern;
00653 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00654
00655 if (gVerbose)
00656 cout << "Comparing m and m1 ..." << endl;
00657 R__ASSERT(m == m1);
00658 if (gVerbose)
00659 cout << "Comparing (m=0) and m1 ..." << endl;
00660 R__ASSERT(!(m.Zero() == m1));
00661
00662 if (gVerbose)
00663 cout << "Clearing m1 ..." << endl;
00664 m1.Zero();
00665 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
00666
00667 if (gVerbose)
00668 cout << "\nClear m and add the pattern" << endl;
00669 m.Zero();
00670 m += pattern;
00671 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
00672 if (gVerbose)
00673 cout << " add the doubled pattern with the negative sign" << endl;
00674 m += -2*pattern;
00675 ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
00676 if (gVerbose)
00677 cout << " subtract the trippled pattern with the negative sign" << endl;
00678 m -= -3*pattern;
00679 ok &= VerifyMatrixValue(m,2*pattern,gVerbose,EPSILON);
00680
00681 if (gVerbose)
00682 cout << "\nVerify comparison operations when all elems are the same" << endl;
00683 m = pattern;
00684 R__ASSERT( m == pattern && !(m != pattern) );
00685 R__ASSERT( m > 0 && m >= pattern && m <= pattern );
00686 R__ASSERT( m > -pattern && m >= -pattern );
00687 R__ASSERT( m <= pattern && !(m < pattern) );
00688 m -= 2*pattern;
00689 R__ASSERT( m < -pattern/2 && m <= -pattern/2 );
00690 R__ASSERT( m >= -pattern && !(m > -pattern) );
00691
00692 if (gVerbose)
00693 cout << "\nVerify comparison operations when not all elems are the same" << endl;
00694 m = pattern; m(m.GetRowUpb(),m.GetColUpb()) = pattern-1;
00695 R__ASSERT( !(m == pattern) && !(m != pattern) );
00696 R__ASSERT( m != 0 );
00697 R__ASSERT( !(m >= pattern) && m <= pattern && !(m<pattern) );
00698 R__ASSERT( !(m <= pattern-1) && m >= pattern-1 && !(m>pattern-1) );
00699
00700 if (gVerbose)
00701 cout << "\nAssign 2*pattern to m by repeating additions" << endl;
00702 m = 0; m += pattern; m += pattern;
00703 if (gVerbose)
00704 cout << "Assign 2*pattern to m1 by multiplying by two " << endl;
00705 m1 = pattern; m1 *= 2;
00706 ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
00707 R__ASSERT( m == m1 );
00708 if (gVerbose)
00709 cout << "Multiply m1 by one half returning it to the 1*pattern" << endl;
00710 m1 *= 1/2.;
00711 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00712
00713 if (gVerbose)
00714 cout << "\nAssign -pattern to m and m1" << endl;
00715 m.Zero(); m -= pattern; m1 = -pattern;
00716 ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
00717 R__ASSERT( m == m1 );
00718 if (gVerbose)
00719 cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same" << endl;
00720 m.Sqr();
00721 ok &= VerifyMatrixValue(m,pattern*pattern,gVerbose,EPSILON);
00722 m.Sqrt();
00723 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
00724 m1.Abs();
00725 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00726 ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
00727
00728 if (gVerbose)
00729 cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1" << endl;
00730 {
00731 #ifndef __CINT__
00732 FillMatrix f(m);
00733 m.Apply(f);
00734 #else
00735 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00736 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
00737 m(i,j) = 4*TMath::Pi()/m.GetNoElements() * (i*m.GetNcols()+j);
00738 #endif
00739 }
00740 m1 = m;
00741 {
00742 #ifndef __CINT__
00743 ApplyFunction s(&TMath::Sin);
00744 ApplyFunction c(&TMath::Cos);
00745 m.Apply(s);
00746 m1.Apply(c);
00747 #else
00748 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
00749 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
00750 m(i,j) = TMath::Sin(m(i,j));
00751 m1(i,j) = TMath::Cos(m1(i,j));
00752 }
00753 }
00754 #endif
00755 }
00756 m.Sqr();
00757 m1.Sqr();
00758 m += m1;
00759 ok &= VerifyMatrixValue(m,1.,gVerbose,EPSILON);
00760
00761 if (gVerbose)
00762 cout << "\nDone\n" << endl;
00763
00764 StatusPrint(3,"Uniform matrix operations",ok);
00765 }
00766
00767
00768
00769
00770
00771 void mstress_binary_ebe_op(Int_t rsize, Int_t csize)
00772 {
00773 if (gVerbose)
00774 cout << "\n---> Test Binary Matrix element-by-element operations" << endl;
00775
00776 Bool_t ok = kTRUE;
00777 const Double_t pattern = 4.25;
00778
00779 TMatrixD m(2,rsize+1,0,csize-1);
00780 TMatrixD m1(TMatrixD::kZero,m);
00781 TMatrixD mp(TMatrixD::kZero,m);
00782
00783 {
00784 for (Int_t i = mp.GetRowLwb(); i <= mp.GetRowUpb(); i++)
00785 for (Int_t j = mp.GetColLwb(); j <= mp.GetColUpb(); j++)
00786 mp(i,j) = (i-m.GetNrows()/2.)*j*pattern;
00787 }
00788
00789 if (gVerbose)
00790 cout << "\nVerify assignment of a matrix to the matrix" << endl;
00791 m = pattern;
00792 m1.Zero();
00793 m1 = m;
00794 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00795 R__ASSERT( m1 == m );
00796
00797 if (gVerbose)
00798 cout << "\nAdding the matrix to itself, uniform pattern " << pattern << endl;
00799 m.Zero(); m = pattern;
00800 m1 = m; m1 += m1;
00801 ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
00802 if (gVerbose)
00803 cout << " subtracting two matrices ..." << endl;
00804 m1 -= m;
00805 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00806 if (gVerbose)
00807 cout << " subtracting the matrix from itself" << endl;
00808 m1 -= m1;
00809 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
00810 if (gVerbose)
00811 cout << " adding two matrices together" << endl;
00812 m1 += m;
00813 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
00814
00815 if (gVerbose) {
00816 cout << "\nArithmetic operations on matrices with not the same elements" << endl;
00817 cout << " adding mp to the zero matrix..." << endl;
00818 }
00819 m.Zero(); m += mp;
00820 ok &= VerifyMatrixIdentity(m,mp,gVerbose,EPSILON);
00821 m1 = m;
00822 if (gVerbose)
00823 cout << " making m = 3*mp and m1 = 3*mp, via add() and succesive mult" << endl;
00824 Add(m,2.,mp);
00825 m1 += m1; m1 += mp;
00826 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
00827 if (gVerbose)
00828 cout << " clear both m and m1, by subtracting from itself and via add()" << endl;
00829 m1 -= m1;
00830 Add(m,-3.,mp);
00831 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
00832
00833 if (gVerbose) {
00834 cout << "\nTesting element-by-element multiplications and divisions" << endl;
00835 cout << " squaring each element with sqr() and via multiplication" << endl;
00836 }
00837 m = mp; m1 = mp;
00838 m.Sqr();
00839 ElementMult(m1,m1);
00840 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
00841 if (gVerbose)
00842 cout << " compare (m = pattern^2)/pattern with pattern" << endl;
00843 m = pattern; m1 = pattern;
00844 m.Sqr();
00845 ElementDiv(m,m1);
00846 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
00847 if (gVerbose)
00848 Compare(m1,m);
00849
00850 if (gVerbose)
00851 cout << "\nDone\n" << endl;
00852
00853 StatusPrint(4,"Binary Matrix element-by-element operations",ok);
00854 }
00855
00856
00857
00858
00859
00860 void mstress_transposition(Int_t msize)
00861 {
00862 if (gVerbose) {
00863 cout << "\n---> Verify matrix transpose "
00864 "for matrices of a characteristic size " << msize << endl;
00865 }
00866
00867 Bool_t ok = kTRUE;
00868 {
00869 if (gVerbose)
00870 cout << "\nCheck to see that a square UnitMatrix stays the same";
00871 TMatrixD m(msize,msize);
00872 m.UnitMatrix();
00873 TMatrixD mt(TMatrixD::kTransposed,m);
00874 ok &= ( m == mt ) ? kTRUE : kFALSE;
00875 }
00876
00877 {
00878 if (gVerbose)
00879 cout << "\nTest a non-square UnitMatrix";
00880 TMatrixD m(msize,msize+1);
00881 m.UnitMatrix();
00882 TMatrixD mt(TMatrixD::kTransposed,m);
00883 R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows() );
00884 for (Int_t i = m.GetRowLwb(); i <= TMath::Min(m.GetRowUpb(),m.GetColUpb()); i++)
00885 for (Int_t j = m.GetColLwb(); j <= TMath::Min(m.GetRowUpb(),m.GetColUpb()); j++)
00886 ok &= ( m(i,j) == mt(i,j) ) ? kTRUE : kFALSE;
00887 }
00888
00889 {
00890 if (gVerbose)
00891 cout << "\nCheck to see that a symmetric (Hilbert)Matrix stays the same";
00892 TMatrixD m = THilbertMatrixD(msize,msize);
00893 TMatrixD mt(TMatrixD::kTransposed,m);
00894 ok &= ( m == mt ) ? kTRUE : kFALSE;
00895 }
00896
00897 {
00898 if (gVerbose)
00899 cout << "\nCheck transposing a non-symmetric matrix";
00900 TMatrixD m = THilbertMatrixD(msize+1,msize);
00901 m(1,2) = TMath::Pi();
00902 TMatrixD mt(TMatrixD::kTransposed,m);
00903 R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows());
00904 R__ASSERT(mt(2,1) == (Double_t)TMath::Pi() && mt(1,2) != (Double_t)TMath::Pi());
00905 R__ASSERT(mt[2][1] == (Double_t)TMath::Pi() && mt[1][2] != (Double_t)TMath::Pi());
00906
00907 if (gVerbose)
00908 cout << "\nCheck double transposing a non-symmetric matrix" << endl;
00909 TMatrixD mtt(TMatrixD::kTransposed,mt);
00910 ok &= ( m == mtt ) ? kTRUE : kFALSE;
00911 }
00912
00913 if (gVerbose)
00914 cout << "\nDone\n" << endl;
00915
00916 StatusPrint(5,"Matrix transposition",ok);
00917 }
00918
00919
00920
00921
00922
00923 class MakeHilbert : public TElementPosActionD {
00924 void Operation(Double_t &element) const { element = 1./(fI+fJ+1); }
00925 public:
00926 MakeHilbert() { }
00927 };
00928
00929 #if !defined (__CINT__) || defined (__MAKECINT__)
00930 class TestUnit : public TElementPosActionD {
00931 mutable Int_t fIsUnit;
00932 void Operation(Double_t &element) const
00933 { if (fIsUnit)
00934 fIsUnit = ((fI==fJ) ? (element == 1.0) : (element == 0)); }
00935 public:
00936 TestUnit() : fIsUnit(0==0) { }
00937 Int_t is_indeed_unit() const { return fIsUnit; }
00938 };
00939 #else
00940 Bool_t is_indeed_unit(TMatrixD &m) {
00941 Bool_t isUnit = kTRUE;
00942 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
00943 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
00944 if (isUnit)
00945 isUnit = ((i==j) ? (m(i,j) == 1.0) : (m(i,j) == 0));
00946 }
00947 return isUnit;
00948 }
00949 #endif
00950
00951 void mstress_special_creation(Int_t dim)
00952 {
00953 if (gVerbose)
00954 cout << "\n---> Check creating some special matrices of dimension " << dim << endl;
00955
00956 Bool_t ok = kTRUE;
00957 {
00958 if (gVerbose)
00959 cout << "\ntest creating Hilbert matrices" << endl;
00960 TMatrixD m = THilbertMatrixD(dim+1,dim);
00961 TMatrixD m1(TMatrixD::kZero,m);
00962 ok &= ( !(m == m1) ) ? kTRUE : kFALSE;
00963 ok &= ( m != 0 ) ? kTRUE : kFALSE;
00964 #ifndef __CINT__
00965 MakeHilbert mh;
00966 m1.Apply(mh);
00967 #else
00968 for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
00969 for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
00970 m1(i,j) = 1./(i+j+1);
00971 #endif
00972 ok &= ( m1 != 0 ) ? kTRUE : kFALSE;
00973 ok &= ( m == m1 ) ? kTRUE : kFALSE;
00974 }
00975
00976 {
00977 if (gVerbose)
00978 cout << "test creating zero matrix and copy constructor" << endl;
00979 TMatrixD m = THilbertMatrixD(dim,dim+1);
00980 ok &= ( m != 0 ) ? kTRUE : kFALSE;
00981 TMatrixD m1(m);
00982 ok &= ( m1 == m ) ? kTRUE : kFALSE;
00983 TMatrixD m2(TMatrixD::kZero,m);
00984 ok &= ( m2 == 0 ) ? kTRUE : kFALSE;
00985 ok &= ( m != 0 ) ? kTRUE : kFALSE;
00986 }
00987
00988 {
00989 if (gVerbose)
00990 cout << "test creating unit matrices" << endl;
00991 TMatrixD m(dim,dim);
00992 #ifndef __CINT__
00993 {
00994 TestUnit test_unit;
00995 m.Apply(test_unit);
00996 ok &= ( !test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
00997 }
00998 #else
00999 ok &= ( !is_indeed_unit(m) ) ? kTRUE : kFALSE;
01000 #endif
01001 m.UnitMatrix();
01002 #ifndef __CINT__
01003 {
01004 TestUnit test_unit;
01005 m.Apply(test_unit);
01006 ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
01007 }
01008 #else
01009 ok &= ( is_indeed_unit(m) ) ? kTRUE : kFALSE;
01010 #endif
01011 m.ResizeTo(dim-1,dim);
01012 TMatrixD m2(TMatrixD::kUnit,m);
01013 #ifndef __CINT__
01014 {
01015 TestUnit test_unit;
01016 m2.Apply(test_unit);
01017 ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
01018 }
01019 #else
01020 ok &= ( is_indeed_unit(m2) ) ? kTRUE : kFALSE;
01021 #endif
01022 m.ResizeTo(dim,dim-2);
01023 m.UnitMatrix();
01024 #ifndef __CINT__
01025 {
01026 TestUnit test_unit;
01027 m.Apply(test_unit);
01028 ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
01029 }
01030 #else
01031 ok &= ( is_indeed_unit(m) ) ? kTRUE : kFALSE;
01032 #endif
01033 }
01034
01035 {
01036 if (gVerbose)
01037 cout << "check to see that Haar matrix has *exactly* orthogonal columns" << endl;
01038 Int_t j;
01039 const Int_t order = 5;
01040 const TMatrixD haar = THaarMatrixD(order);
01041 ok &= ( haar.GetNrows() == (1<<order) &&
01042 haar.GetNrows() == haar.GetNcols() ) ? kTRUE : kFALSE;
01043 TVectorD colj(1<<order);
01044 TVectorD coll(1<<order);
01045 for (j = haar.GetColLwb(); j <= haar.GetColUpb(); j++) {
01046 colj = TMatrixDColumn_const(haar,j);
01047 ok &= (TMath::Abs(colj*colj-1.0) <= 1.0e-15 ) ? kTRUE : kFALSE;
01048 for (Int_t l = j+1; l <= haar.GetColUpb(); l++) {
01049 coll = TMatrixDColumn_const(haar,l);
01050 const Double_t val = colj*coll;
01051 ok &= ( TMath::Abs(val) <= 1.0e-15 ) ? kTRUE : kFALSE;
01052 }
01053 }
01054
01055 if (gVerbose)
01056 cout << "make Haar (sub)matrix and test it *is* a submatrix" << endl;
01057 const Int_t no_sub_cols = (1<<order) - 3;
01058 const TMatrixD haar_sub = THaarMatrixD(order,no_sub_cols);
01059 ok &= ( haar_sub.GetNrows() == (1<<order) &&
01060 haar_sub.GetNcols() == no_sub_cols ) ? kTRUE : kFALSE;
01061 for (j = haar_sub.GetColLwb(); j <= haar_sub.GetColUpb(); j++) {
01062 colj = TMatrixDColumn_const(haar,j);
01063 coll = TMatrixDColumn_const(haar_sub,j);
01064 ok &= VerifyVectorIdentity(colj,coll,gVerbose);
01065 }
01066 }
01067
01068 if (gVerbose)
01069 cout << "\nDone\n" << endl;
01070
01071 StatusPrint(6,"Haar/Hilbert Matrix",ok);
01072 }
01073
01074 #ifndef __CINT__
01075
01076
01077
01078
01079 class hilbert_matrix_promise : public TMatrixDLazy {
01080 void FillIn(TMatrixD &m) const { m = THilbertMatrixD(m.GetRowLwb(),m.GetRowUpb(),
01081 m.GetColLwb(),m.GetColUpb()); }
01082
01083 public:
01084 hilbert_matrix_promise(Int_t nrows,Int_t ncols)
01085 : TMatrixDLazy(nrows,ncols) {}
01086 hilbert_matrix_promise(Int_t row_lwb,Int_t row_upb,
01087 Int_t col_lwb,Int_t col_upb)
01088 : TMatrixDLazy(row_lwb,row_upb,col_lwb,col_upb) { }
01089 };
01090
01091 void mstress_matrix_promises(Int_t dim)
01092 {
01093 if (gVerbose)
01094 cout << "\n---> Check making/forcing promises, (lazy)matrices of dimension " << dim << endl;
01095
01096 Bool_t ok = kTRUE;
01097 {
01098 if (gVerbose)
01099 cout << "\nmake a promise and force it by a constructor" << endl;
01100 TMatrixD m = hilbert_matrix_promise(dim,dim+1);
01101 TMatrixD m1 = THilbertMatrixD(dim,dim+1);
01102 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
01103 }
01104
01105 {
01106 if (gVerbose)
01107 cout << "make a promise and force it by an assignment" << endl;
01108 TMatrixD m(-1,dim,0,dim);
01109 m = hilbert_matrix_promise(-1,dim,0,dim);
01110 TMatrixD m1 = THilbertMatrixD(-1,dim,0,dim);
01111 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
01112 }
01113
01114 if (gVerbose)
01115 cout << "\nDone\n" << endl;
01116
01117 StatusPrint(7,"Matrix promises",ok);
01118 }
01119 #endif
01120
01121
01122
01123
01124
01125 void mstress_norms(Int_t rsize_req,Int_t csize_req)
01126 {
01127 if (gVerbose)
01128 cout << "\n---> Verify norm calculations" << endl;
01129
01130 Bool_t ok = kTRUE;
01131 const Double_t pattern = 10.25;
01132
01133 Int_t rsize = rsize_req;
01134 if (rsize%2 != 0) rsize--;
01135 Int_t csize = csize_req;
01136 if (csize%2 != 0) csize--;
01137 if (rsize%2 == 1 || csize%2 == 1) {
01138 cout << "rsize: " << rsize <<endl;
01139 cout << "csize: " << csize <<endl;
01140 Fatal("mstress_norms","Sorry, size of the matrix to test must be even for this test\n");
01141 }
01142
01143 TMatrixD m(rsize,csize);
01144
01145 if (gVerbose)
01146 cout << "\nAssign " << pattern << " to all the elements and check norms" << endl;
01147 m = pattern;
01148 if (gVerbose)
01149 cout << " 1. (col) norm should be pattern*nrows" << endl;
01150 ok &= ( m.Norm1() == pattern*m.GetNrows() ) ? kTRUE : kFALSE;
01151 ok &= ( m.Norm1() == m.ColNorm() ) ? kTRUE : kFALSE;
01152 if (gVerbose)
01153 cout << " Inf (row) norm should be pattern*ncols" << endl;
01154 ok &= ( m.NormInf() == pattern*m.GetNcols() ) ? kTRUE : kFALSE;
01155 ok &= ( m.NormInf() == m.RowNorm() ) ? kTRUE : kFALSE;
01156 if (gVerbose)
01157 cout << " Square of the Eucl norm has got to be pattern^2 * no_elems" << endl;
01158 ok &= ( m.E2Norm() == (pattern*pattern)*m.GetNoElements() ) ? kTRUE : kFALSE;
01159 TMatrixD m1(TMatrixD::kZero,m);
01160 ok &= ( m.E2Norm() == E2Norm(m,m1) ) ? kTRUE : kFALSE;
01161
01162 if (gVerbose)
01163 cout << "\nDone\n" << endl;
01164
01165 StatusPrint(8,"Matrix Norms",ok);
01166 }
01167
01168
01169
01170
01171
01172 void mstress_determinant(Int_t msize)
01173 {
01174 if (gVerbose)
01175 cout << "\n---> Verify determinant evaluation for a square matrix of size " << msize << endl;
01176
01177 Bool_t ok = kTRUE;
01178 TMatrixD m(msize,msize);
01179 const Double_t pattern = 2.5;
01180
01181 {
01182 if (gVerbose)
01183 cout << "\nCheck to see that the determinant of the unit matrix is one" <<endl;
01184 m.UnitMatrix();
01185 Double_t d1,d2;
01186 m.Determinant(d1,d2);
01187 const Double_t det = d1*TMath::Power(2.,d2);
01188 if (gVerbose) {
01189 cout << "det = " << det << " deviation= " << TMath::Abs(det-1);
01190 cout << ( (TMath::Abs(det-1) < EPSILON) ? " OK" : " too large") <<endl;
01191 }
01192 ok &= (TMath::Abs(det-1) < EPSILON) ? kTRUE : kFALSE;
01193 }
01194
01195 if (ok)
01196 {
01197 if (gVerbose)
01198 cout << "\nCheck the determinant for the matrix with " << pattern << " at the diagonal" << endl;
01199 {
01200 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
01201 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
01202 m(i,j) = ( i==j ? pattern : 0 );
01203 }
01204 Double_t d1,d2;
01205 m.Determinant(d1,d2);
01206 const Double_t d1_abs = TMath::Abs(d1);
01207 const Double_t log_det1 = d2*TMath::Log(2.)+TMath::Log(d1_abs);
01208 const Double_t log_det2 = ((Double_t)m.GetNrows())*TMath::Log(pattern);
01209 if (gVerbose) {
01210 cout << "log_det1 = " << log_det1 << " deviation= " << TMath::Abs(log_det1-log_det2);
01211 cout << ( (TMath::Abs(log_det1-log_det2) < msize*EPSILON) ? " OK" : " too large") <<endl;
01212 }
01213 ok &= ( TMath::Abs(log_det1-log_det2) < msize*EPSILON ) ? kTRUE : kFALSE;
01214 }
01215
01216 if (ok)
01217 {
01218 if (gVerbose)
01219 cout << "\nCheck the determinant of the transposed matrix" << endl;
01220 m.UnitMatrix();
01221 m(1,2) = pattern;
01222 TMatrixD m_tran(TMatrixD::kTransposed,m);
01223 ok &= ( !(m == m_tran) ) ? kTRUE : kFALSE;
01224 const Double_t det1 = m.Determinant();
01225 const Double_t det2 = m_tran.Determinant();
01226 if (gVerbose) {
01227 cout << "det1 = " << det1 << " deviation= " << TMath::Abs(det1-det2);
01228 cout << ( (TMath::Abs(det1-det2) < msize*EPSILON) ? " OK" : " too large") <<endl;
01229 }
01230 ok &= ( TMath::Abs(det1-det2) < msize*EPSILON ) ? kTRUE : kFALSE;
01231 }
01232
01233 if (ok)
01234 {
01235 if (gVerbose)
01236 cout << "\nswap two rows/cols of a matrix through method 1 and watch det's sign" << endl;
01237 m.UnitMatrix();
01238 TMatrixDRow(m,3) = pattern;
01239 Double_t d1,d2;
01240 m.Determinant(d1,d2);
01241 TMatrixDRow row1(m,1);
01242 TVectorD vrow1(m.GetRowLwb(),m.GetRowUpb()); vrow1 = row1;
01243 TVectorD vrow3(m.GetRowLwb(),m.GetRowUpb()); vrow3 = TMatrixDRow(m,3);
01244 row1 = vrow3; TMatrixDRow(m,3) = vrow1;
01245 Double_t d1_p,d2_p;
01246 m.Determinant(d1_p,d2_p);
01247 if (gVerbose) {
01248 cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1+d1_p);
01249 cout << ( ( d1 == -d1_p ) ? " OK" : " too large") <<endl;
01250 }
01251 ok &= ( d1 == -d1_p ) ? kTRUE : kFALSE;
01252 TMatrixDColumn col2(m,2);
01253 TVectorD vcol2(m.GetRowLwb(),m.GetRowUpb()); vcol2 = col2;
01254 TVectorD vcol4(m.GetRowLwb(),m.GetRowUpb()); vcol4 = TMatrixDColumn(m,4);
01255 col2 = vcol4; TMatrixDColumn(m,4) = vcol2;
01256 m.Determinant(d1_p,d2_p);
01257 if (gVerbose) {
01258 cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1-d1_p);
01259 cout << ( ( d1 == d1_p ) ? " OK" : " too large") <<endl;
01260 }
01261 ok &= ( d1 == d1_p ) ? kTRUE : kFALSE;
01262 }
01263
01264 if (ok)
01265 {
01266 if (gVerbose)
01267 cout << "\nswap two rows/cols of a matrix through method 2 and watch det's sign" << endl;
01268 m.UnitMatrix();
01269 TMatrixDRow(m,3) = pattern;
01270 Double_t d1,d2;
01271 m.Determinant(d1,d2);
01272
01273 TMatrixD m_save( m);
01274 TMatrixDRow(m,1) = TMatrixDRow(m_save,3);
01275 TMatrixDRow(m,3) = TMatrixDRow(m_save,1);
01276 Double_t d1_p,d2_p;
01277 m.Determinant(d1_p,d2_p);
01278 if (gVerbose) {
01279 cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1+d1_p);
01280 cout << ( ( d1 == -d1_p ) ? " OK" : " too large") <<endl;
01281 }
01282 ok &= ( d1 == -d1_p ) ? kTRUE : kFALSE;
01283
01284 m_save = m;
01285 TMatrixDColumn(m,2) = TMatrixDColumn(m_save,4);
01286 TMatrixDColumn(m,4) = TMatrixDColumn(m_save,2);
01287 m.Determinant(d1_p,d2_p);
01288 if (gVerbose) {
01289 cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1-d1_p);
01290 cout << ( ( d1 == d1_p ) ? " OK" : " too large") <<endl;
01291 }
01292 ok &= ( d1 == d1_p ) ? kTRUE : kFALSE;
01293 }
01294
01295 if (ok)
01296 {
01297 if (gVerbose)
01298 cout << "\nCheck the determinant for the matrix with " << pattern << " at the anti-diagonal" << endl;
01299 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
01300 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
01301 m(i,j) = ( i==(m.GetColUpb()+m.GetColLwb()-j) ? pattern : 0 );
01302
01303 Double_t d1,d2;
01304 m.Determinant(d1,d2);
01305 const Double_t d1_abs = TMath::Abs(d1);
01306 const Double_t log_det1 = d2*TMath::Log(2.)+TMath::Log(d1_abs);
01307 const Double_t log_det2 = ((Double_t)m.GetNrows())*TMath::Log(pattern);
01308 const Double_t sign = ( m.GetNrows()*(m.GetNrows()-1)/2 & 1 ? -1 : 1 );
01309 if (gVerbose) {
01310 cout << "log_det1 = " << log_det1 << " deviation= " << TMath::Abs(log_det1-log_det2);
01311 cout << ( (TMath::Abs(log_det1-log_det2) < msize*EPSILON) ? " OK" : " too large") <<endl;
01312 cout << ( sign * d1 > 0. ? " OK sign" : " wrong sign") <<endl;
01313 }
01314 ok &= ( TMath::Abs(log_det1-log_det2) < msize*EPSILON ) ? kTRUE : kFALSE;
01315 ok &= ( sign * d1 > 0. ) ? kTRUE : kFALSE;
01316 }
01317
01318
01319 if (0) {
01320 if (gVerbose)
01321 cout << "\nCheck the determinant for the singular matrix"
01322 "\n defined as above with zero first row" << endl;
01323 m.Zero();
01324 {
01325 for (Int_t i = m.GetRowLwb()+1; i <= m.GetRowUpb(); i++)
01326 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
01327 m(i,j) = ( i==(m.GetColUpb()+m.GetColLwb()-j) ? pattern : 0 );
01328 }
01329 cout << "\n determinant is " << m.Determinant();
01330 ok &= ( m.Determinant() == 0 ) ? kTRUE : kFALSE;
01331 }
01332
01333 if (gVerbose)
01334 cout << "\nCheck out the determinant of the Hilbert matrix";
01335 TMatrixDSym H = THilbertMatrixDSym(3);
01336 H.SetTol(1.0e-20);
01337 if (gVerbose) {
01338 cout << "\n 3x3 Hilbert matrix: exact determinant 1/2160 ";
01339 cout << "\n computed 1/"<< 1/H.Determinant();
01340 }
01341
01342 H.ResizeTo(4,4);
01343 H = THilbertMatrixDSym(4);
01344 H.SetTol(1.0e-20);
01345 if (gVerbose) {
01346 cout << "\n 4x4 Hilbert matrix: exact determinant 1/6048000 ";
01347 cout << "\n computed 1/"<< 1/H.Determinant();
01348 }
01349
01350 H.ResizeTo(5,5);
01351 H = THilbertMatrixDSym(5);
01352 H.SetTol(1.0e-20);
01353 if (gVerbose) {
01354 cout << "\n 5x5 Hilbert matrix: exact determinant 3.749295e-12";
01355 cout << "\n computed "<< H.Determinant();
01356 }
01357
01358 if (gVerbose) {
01359 TDecompQRH qrh(H);
01360 Double_t d1,d2;
01361 qrh.Det(d1,d2);
01362 cout << "\n qrh det = " << d1*TMath::Power(2.0,d2) <<endl;
01363 }
01364
01365 if (gVerbose) {
01366 TDecompChol chol(H);
01367 Double_t d1,d2;
01368 chol.Det(d1,d2);
01369 cout << "\n chol det = " << d1*TMath::Power(2.0,d2) <<endl;
01370 }
01371
01372 if (gVerbose) {
01373 TDecompSVD svd(H);
01374 Double_t d1,d2;
01375 svd.Det(d1,d2);
01376 cout << "\n svd det = " << d1*TMath::Power(2.0,d2) <<endl;
01377 }
01378
01379 H.ResizeTo(7,7);
01380 H = THilbertMatrixDSym(7);
01381 H.SetTol(1.0e-20);
01382 if (gVerbose) {
01383 cout << "\n 7x7 Hilbert matrix: exact determinant 4.8358e-25";
01384 cout << "\n computed "<< H.Determinant();
01385 }
01386
01387 H.ResizeTo(9,9);
01388 H = THilbertMatrixDSym(9);
01389 H.SetTol(1.0e-20);
01390 if (gVerbose) {
01391 cout << "\n 9x9 Hilbert matrix: exact determinant 9.72023e-43";
01392 cout << "\n computed "<< H.Determinant();
01393 }
01394
01395 H.ResizeTo(10,10);
01396 H = THilbertMatrixDSym(10);
01397 H.SetTol(1.0e-20);
01398 if (gVerbose) {
01399 cout << "\n 10x10 Hilbert matrix: exact determinant 2.16418e-53";
01400 cout << "\n computed "<< H.Determinant();
01401 }
01402
01403 if (gVerbose)
01404 cout << "\nDone\n" << endl;
01405
01406 StatusPrint(9,"Matrix Determinant",ok);
01407 }
01408
01409
01410
01411
01412
01413 void mstress_mm_multiplications()
01414 {
01415 Bool_t ok = kTRUE;
01416
01417 Int_t iloop = 0;
01418 Int_t nr = 0;
01419 while (iloop <= gNrLoop) {
01420 const Int_t msize = gSizeA[iloop];
01421 const Double_t epsilon = EPSILON*msize/100;
01422
01423 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
01424
01425 Int_t i,j;
01426 if (msize <= 5) {
01427 iloop++;
01428 continue;
01429 }
01430
01431 if (verbose)
01432 cout << "\n---> Verify matrix multiplications "
01433 "for matrices of the characteristic size " << msize << endl;
01434
01435 {
01436 if (verbose)
01437 cout << "\nTest inline multiplications of the UnitMatrix" << endl;
01438 TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
01439 TMatrixD u(TMatrixD::kUnit,m);
01440 m(3,1) = TMath::Pi();
01441 u *= m;
01442 ok &= VerifyMatrixIdentity(u,m,verbose,epsilon);
01443 }
01444
01445 if (ok)
01446 {
01447 if (verbose)
01448 cout << "Test inline multiplications by a DiagMat" << endl;
01449 TMatrixD m = THilbertMatrixD(msize+3,msize);
01450 m(1,3) = TMath::Pi();
01451 TVectorD v(msize);
01452 for (i = v.GetLwb(); i <= v.GetUpb(); i++)
01453 v(i) = 1+i;
01454 TMatrixD diag(msize,msize);
01455 TMatrixDDiag d(diag);
01456 d = v;
01457 TMatrixD eth = m;
01458 for (i = eth.GetRowLwb(); i <= eth.GetRowUpb(); i++)
01459 for (j = eth.GetColLwb(); j <= eth.GetColUpb(); j++)
01460 eth(i,j) *= v(j);
01461 m *= diag;
01462 ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
01463 }
01464
01465 if (ok)
01466 {
01467 if (verbose)
01468 cout << "Test XPP = X where P is a permutation matrix" << endl;
01469 TMatrixD m = THilbertMatrixD(msize-1,msize);
01470 m(2,3) = TMath::Pi();
01471 TMatrixD eth = m;
01472 TMatrixD p(msize,msize);
01473 for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++)
01474 p(p.GetRowUpb()+p.GetRowLwb()-i,i) = 1;
01475 m *= p;
01476 m *= p;
01477 ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
01478 }
01479
01480 if (ok)
01481 {
01482 if (verbose)
01483 cout << "Test general matrix multiplication through inline mult" << endl;
01484 TMatrixD m = THilbertMatrixD(msize-2,msize);
01485 m(3,3) = TMath::Pi();
01486 TMatrixD mt(TMatrixD::kTransposed,m);
01487 TMatrixD p = THilbertMatrixD(msize,msize);
01488 TMatrixDDiag(p) += 1;
01489 TMatrixD mp(m,TMatrixD::kMult,p);
01490 TMatrixD m1 = m;
01491 m *= p;
01492 ok &= VerifyMatrixIdentity(m,mp,verbose,epsilon);
01493 TMatrixD mp1(mt,TMatrixD::kTransposeMult,p);
01494 VerifyMatrixIdentity(m,mp1,verbose,epsilon);
01495 ok &= ( !(m1 == m) );
01496 TMatrixD mp2(TMatrixD::kZero,m1);
01497 ok &= ( mp2 == 0 );
01498 mp2.Mult(m1,p);
01499 ok &= VerifyMatrixIdentity(m,mp2,verbose,epsilon);
01500
01501 if (verbose)
01502 cout << "Test XP=X*P vs XP=X;XP*=P" << endl;
01503 TMatrixD mp3 = m1*p;
01504 ok &= VerifyMatrixIdentity(m,mp3,verbose,epsilon);
01505 }
01506
01507 if (ok)
01508 {
01509 if (verbose)
01510 cout << "Check n * m == n * m_sym; m_sym * n == m * n; m_sym * m_sym == m * m" <<endl;
01511
01512 const TMatrixD n = THilbertMatrixD(0,msize-1,0,msize-1);
01513 const TMatrixD m = n;
01514 const TMatrixDSym m_sym = THilbertMatrixDSym(0,msize-1);
01515
01516 const TMatrixD nm1 = n * m_sym;
01517 const TMatrixD nm2 = n * m;
01518 const TMatrixD mn1 = m_sym * n;
01519 const TMatrixD mn2 = m * n;
01520 const TMatrixD mm1 = m_sym * m_sym;
01521 const TMatrixD mm2 = m * m;
01522
01523 ok &= VerifyMatrixIdentity(nm1,nm2,verbose,epsilon);
01524 ok &= VerifyMatrixIdentity(mn1,mn2,verbose,epsilon);
01525 ok &= VerifyMatrixIdentity(mm1,mm2,verbose,epsilon);
01526 }
01527
01528 #ifndef __CINT__
01529 if (nr >= Int_t(1.0e+5/msize/msize)) {
01530 #else
01531 if (nr >= Int_t(1.0e+3/msize/msize)) {
01532 #endif
01533 nr = 0;
01534 iloop++;
01535 } else
01536 nr++;
01537
01538 if (!ok) {
01539 if (gVerbose)
01540 cout << "General Matrix Multiplications failed for size " << msize << endl;
01541 break;
01542 }
01543 }
01544
01545 if (ok)
01546 {
01547 if (gVerbose)
01548 cout << "Check to see UU' = U'U = E when U is the Haar matrix" << endl;
01549 const Int_t order = 5;
01550 const Int_t no_sub_cols = (1<<order)-5;
01551 TMatrixD haar_sub = THaarMatrixD(5,no_sub_cols);
01552 TMatrixD haar_sub_t(TMatrixD::kTransposed,haar_sub);
01553 TMatrixD hsths(haar_sub_t,TMatrixD::kMult,haar_sub);
01554 TMatrixD hsths1(TMatrixD::kZero,hsths); hsths1.Mult(haar_sub_t,haar_sub);
01555 TMatrixD hsths_eth(TMatrixD::kUnit,hsths);
01556 ok &= ( hsths.GetNrows() == no_sub_cols && hsths.GetNcols() == no_sub_cols );
01557 ok &= VerifyMatrixIdentity(hsths,hsths_eth,gVerbose,EPSILON);
01558 ok &= VerifyMatrixIdentity(hsths1,hsths_eth,gVerbose,EPSILON);
01559 TMatrixD haar = THaarMatrixD(order);
01560 TMatrixD unit(TMatrixD::kUnit,haar);
01561 TMatrixD haar_t(TMatrixD::kTransposed,haar);
01562 TMatrixD hth(haar,TMatrixD::kTransposeMult,haar);
01563 TMatrixD hht(haar,TMatrixD::kMult,haar_t);
01564 TMatrixD hht1 = haar; hht1 *= haar_t;
01565 TMatrixD hht2(TMatrixD::kZero,haar); hht2.Mult(haar,haar_t);
01566 ok &= VerifyMatrixIdentity(unit,hth,gVerbose,EPSILON);
01567 ok &= VerifyMatrixIdentity(unit,hht,gVerbose,EPSILON);
01568 ok &= VerifyMatrixIdentity(unit,hht1,gVerbose,EPSILON);
01569 ok &= VerifyMatrixIdentity(unit,hht2,gVerbose,EPSILON);
01570 }
01571
01572 if (gVerbose)
01573 cout << "\nDone\n" << endl;
01574
01575 StatusPrint(10,"General Matrix Multiplications",ok);
01576 }
01577
01578
01579
01580
01581
01582 void mstress_sym_mm_multiplications(Int_t msize)
01583 {
01584 if (gVerbose)
01585 cout << "\n---> Verify symmetric matrix multiplications "
01586 "for matrices of the characteristic size " << msize << endl;
01587
01588 Int_t i,j;
01589 Bool_t ok = kTRUE;
01590
01591 const Double_t epsilon = EPSILON*msize/100;
01592
01593 {
01594 if (gVerbose)
01595 cout << "\nTest inline multiplications of the UnitMatrix" << endl;
01596 TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
01597 TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
01598 TMatrixDSym u(TMatrixDSym::kUnit,m_sym);
01599 TMatrixD u2 = u * m_sym;
01600 ok &= VerifyMatrixIdentity(u2,m_sym,gVerbose,epsilon);
01601 }
01602
01603 if (ok)
01604 {
01605 if (gVerbose)
01606 cout << "\nTest symmetric multiplications" << endl;
01607 {
01608 if (gVerbose)
01609 cout << "\n Test m * m_sym == m_sym * m == m_sym * m_sym multiplications" << endl;
01610 TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
01611 TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
01612 TMatrixD mm = m * m;
01613 TMatrixD mm_sym1 = m_sym * m_sym;
01614 TMatrixD mm_sym2 = m * m_sym;
01615 TMatrixD mm_sym3 = m_sym * m;
01616 ok &= VerifyMatrixIdentity(mm,mm_sym1,gVerbose,epsilon);
01617 ok &= VerifyMatrixIdentity(mm,mm_sym2,gVerbose,epsilon);
01618 ok &= VerifyMatrixIdentity(mm,mm_sym3,gVerbose,epsilon);
01619 }
01620
01621 {
01622 if (gVerbose)
01623 cout << "\n Test m^T * m_sym == m_sym^T * m == m_sym^T * m_sym multiplications" << endl;
01624 TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
01625 TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
01626 TMatrixD mtm = TMatrixD(m ,TMatrixD::kTransposeMult,m);
01627 TMatrixD mtm_sym1 = TMatrixD(m_sym,TMatrixD::kTransposeMult,m_sym);
01628 TMatrixD mtm_sym2 = TMatrixD(m ,TMatrixD::kTransposeMult,m_sym);
01629 TMatrixD mtm_sym3 = TMatrixD(m_sym,TMatrixD::kTransposeMult,m);
01630 ok &= VerifyMatrixIdentity(mtm,mtm_sym1,gVerbose,epsilon);
01631 ok &= VerifyMatrixIdentity(mtm,mtm_sym2,gVerbose,epsilon);
01632 ok &= VerifyMatrixIdentity(mtm,mtm_sym3,gVerbose,epsilon);
01633 }
01634
01635 {
01636 if (gVerbose)
01637 cout << "\n Test m * m_sym^T == m_sym * m^T == m_sym * m_sym^T multiplications" << endl;
01638 TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
01639 TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
01640 TMatrixD mmt = TMatrixD(m ,TMatrixD::kMultTranspose,m);
01641 TMatrixD mmt_sym1 = TMatrixD(m_sym,TMatrixD::kMultTranspose,m_sym);
01642 TMatrixD mmt_sym2 = TMatrixD(m ,TMatrixD::kMultTranspose,m_sym);
01643 TMatrixD mmt_sym3 = TMatrixD(m_sym,TMatrixD::kMultTranspose,m);
01644 ok &= VerifyMatrixIdentity(mmt,mmt_sym1,gVerbose,epsilon);
01645 ok &= VerifyMatrixIdentity(mmt,mmt_sym2,gVerbose,epsilon);
01646 ok &= VerifyMatrixIdentity(mmt,mmt_sym3,gVerbose,epsilon);
01647 }
01648
01649 {
01650 if (gVerbose)
01651 cout << "\n Test n * m_sym == n * m multiplications" << endl;
01652 TMatrixD n = THilbertMatrixD(-1,msize,-1,msize);
01653 TMatrixD m = n;
01654 n(1,3) = TMath::Pi();
01655 n(3,1) = TMath::Pi();
01656 TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
01657 TMatrixD nm1 = n * m_sym;
01658 TMatrixD nm2 = n * m;
01659 ok &= VerifyMatrixIdentity(nm1,nm2,gVerbose,epsilon);
01660 }
01661 }
01662
01663 if (ok)
01664 {
01665 if (gVerbose)
01666 cout << "Test inline multiplications by a DiagMatrix" << endl;
01667 TMatrixDSym ms = THilbertMatrixDSym(msize);
01668 ms(1,3) = TMath::Pi();
01669 ms(3,1) = TMath::Pi();
01670 TVectorD v(msize);
01671 for (i = v.GetLwb(); i <= v.GetUpb(); i++)
01672 v(i) = 1+i;
01673 TMatrixDSym diag(msize);
01674 TMatrixDDiag d(diag); d = v;
01675 TMatrixDSym eth = ms;
01676 for (i = eth.GetRowLwb(); i <= eth.GetRowUpb(); i++)
01677 for (j = eth.GetColLwb(); j <= eth.GetColUpb(); j++)
01678 eth(i,j) *= v(j);
01679 TMatrixD m2 = ms * diag;
01680 ok &= VerifyMatrixIdentity(m2,eth,gVerbose,epsilon);
01681 }
01682
01683 if (ok)
01684 {
01685 if (gVerbose)
01686 cout << "Test XPP = X where P is a permutation matrix" << endl;
01687 TMatrixDSym ms = THilbertMatrixDSym(msize);
01688 ms(2,3) = TMath::Pi();
01689 ms(3,2) = TMath::Pi();
01690 TMatrixDSym eth = ms;
01691 TMatrixDSym p(msize);
01692 for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++)
01693 p(p.GetRowUpb()+p.GetRowLwb()-i,i) = 1;
01694 TMatrixD m2 = ms * p;
01695 m2 *= p;
01696 ok &= VerifyMatrixIdentity(m2,eth,gVerbose,epsilon);
01697 }
01698
01699 if (ok)
01700 {
01701 if (gVerbose)
01702 cout << "Test general matrix multiplication through inline mult" << endl;
01703 TMatrixDSym ms = THilbertMatrixDSym(msize);
01704 ms(2,3) = TMath::Pi();
01705 ms(3,2) = TMath::Pi();
01706 TMatrixDSym mt(TMatrixDSym::kTransposed,ms);
01707 TMatrixDSym p = THilbertMatrixDSym(msize);
01708 TMatrixDDiag(p) += 1;
01709 TMatrixD mp(ms,TMatrixD::kMult,p);
01710 TMatrixDSym m1 = ms;
01711 TMatrixD m3(ms,TMatrixD::kMult,p);
01712 memcpy(ms.GetMatrixArray(),m3.GetMatrixArray(),msize*msize*sizeof(Double_t));
01713 ok &= VerifyMatrixIdentity(ms,mp,gVerbose,epsilon);
01714 TMatrixD mp1(mt,TMatrixD::kTransposeMult,p);
01715 ok &= VerifyMatrixIdentity(ms,mp1,gVerbose,epsilon);
01716 ok &= ( !(m1 == ms) ) ? kTRUE : kFALSE;
01717 TMatrixDSym mp2(TMatrixDSym::kZero,ms);
01718 ok &= ( mp2 == 0 ) ? kTRUE : kFALSE;
01719
01720 if (gVerbose)
01721 cout << "Test XP=X*P vs XP=X;XP*=P" << endl;
01722 TMatrixD mp3 = m1*p;
01723 ok &= VerifyMatrixIdentity(ms,mp3,gVerbose,epsilon);
01724 }
01725
01726 if (ok)
01727 {
01728 if (gVerbose)
01729 cout << "Check to see UU' = U'U = E when U is the Haar matrix" << endl;
01730 {
01731 const Int_t order = 5;
01732 const Int_t no_sub_cols = (1<<order)-5;
01733 TMatrixD haarb = THaarMatrixD(5,no_sub_cols);
01734 TMatrixD haarb_t(TMatrixD::kTransposed,haarb);
01735 TMatrixD hth(haarb_t,TMatrixD::kMult,haarb);
01736 TMatrixDSym hth1(TMatrixDSym::kAtA,haarb);
01737 ok &= VerifyMatrixIdentity(hth,hth1,gVerbose,epsilon);
01738 }
01739
01740 {
01741 TMatrixD haar = THaarMatrixD(5);
01742 TMatrixD unit(TMatrixD::kUnit,haar);
01743 TMatrixD haar_t(TMatrixD::kTransposed,haar);
01744 TMatrixDSym hths(TMatrixDSym::kAtA,haar);
01745 TMatrixD hht(haar,TMatrixD::kMult,haar_t);
01746 TMatrixD hht1 = haar; hht1 *= haar_t;
01747 ok &= VerifyMatrixIdentity(unit,hths,gVerbose,epsilon);
01748 ok &= VerifyMatrixIdentity(unit,hht,gVerbose,epsilon);
01749 ok &= VerifyMatrixIdentity(unit,hht1,gVerbose,epsilon);
01750 }
01751 }
01752
01753 if (ok)
01754 {
01755 if (gVerbose)
01756 cout << "Check to see A.Similarity(Haar) = Haar * A * Haar^T" <<
01757 " and A.SimilarityT(Haar) = Haar^T * A * Haar" << endl;
01758 {
01759 TMatrixD h = THaarMatrixD(5);
01760 TMatrixDSym ms = THilbertMatrixDSym(1<<5);
01761 TMatrixD hmht = h*TMatrixD(ms,TMatrixD::kMultTranspose,h);
01762 ok &= VerifyMatrixIdentity(ms.Similarity(h),hmht,gVerbose,epsilon);
01763 TMatrixD htmh = TMatrixD(h,TMatrixD::kTransposeMult,ms)*h;
01764 ok &= VerifyMatrixIdentity(ms.SimilarityT(h),htmh,gVerbose,epsilon);
01765 }
01766 if (gVerbose)
01767 cout << "Check to see A.Similarity(B_sym) = A.Similarity(B)" << endl;
01768 {
01769 TMatrixDSym nsym = THilbertMatrixDSym(5);
01770 TMatrixD n = THilbertMatrixD(5,5);
01771 TMatrixDSym ms = THilbertMatrixDSym(5);
01772 ok &= VerifyMatrixIdentity(ms.Similarity(nsym),ms.Similarity(n),gVerbose,epsilon);
01773 }
01774 }
01775
01776 if (gVerbose)
01777 cout << "\nDone\n" << endl;
01778
01779 StatusPrint(11,"Symmetric Matrix Multiplications",ok);
01780 }
01781
01782
01783
01784
01785
01786 void mstress_vm_multiplications()
01787 {
01788 Bool_t ok = kTRUE;
01789
01790 Int_t iloop = gNrLoop;
01791 Int_t nr = 0;
01792 while (iloop >= 0) {
01793 const Int_t msize = gSizeA[iloop];
01794
01795 const Double_t epsilon = EPSILON*msize;
01796
01797 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
01798
01799 if (verbose)
01800 cout << "\n---> Verify vector-matrix multiplications "
01801 "for matrices of the characteristic size " << msize << endl;
01802
01803 {
01804 if (verbose)
01805 cout << "\nCheck shrinking a vector by multiplying by a non-sq unit matrix" << endl;
01806 TVectorD vb(-2,msize);
01807 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
01808 vb(i) = TMath::Pi()-i;
01809 ok &= ( vb != 0 ) ? kTRUE : kFALSE;
01810 TMatrixD mc(1,msize-2,-2,msize);
01811 mc.UnitMatrix();
01812 TVectorD v1 = vb;
01813 TVectorD v2 = vb;
01814 v1 *= mc;
01815 v2.ResizeTo(1,msize-2);
01816 ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
01817 }
01818
01819 if (ok)
01820 {
01821 if (verbose)
01822 cout << "Check expanding a vector by multiplying by a non-sq unit matrix" << endl;
01823 TVectorD vb(msize);
01824 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
01825 vb(i) = TMath::Pi()+i;
01826 ok &= ( vb != 0 ) ? kTRUE : kFALSE;
01827 TMatrixD me(2,msize+5,0,msize-1);
01828 me.UnitMatrix();
01829 TVectorD v1 = vb;
01830 TVectorD v2 = vb;
01831 v1 *= me;
01832 v2.ResizeTo(v1);
01833 ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
01834 }
01835
01836 if (ok)
01837 {
01838 if (verbose)
01839 cout << "Check general matrix-vector multiplication" << endl;
01840 TVectorD vb(msize);
01841 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
01842 vb(i) = TMath::Pi()+i;
01843 TMatrixD vm(msize,1);
01844 TMatrixDColumn(vm,0) = vb;
01845 TMatrixD m = THilbertMatrixD(0,msize,0,msize-1);
01846 vb *= m;
01847 ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
01848 TMatrixD mvm(m,TMatrixD::kMult,vm);
01849 TMatrixD mvb(msize+1,1);
01850 TMatrixDColumn(mvb,0) = vb;
01851 ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
01852 }
01853
01854 if (ok)
01855 {
01856 if (verbose)
01857 cout << "Check symmetric matrix-vector multiplication" << endl;
01858 TVectorD vb(msize);
01859 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
01860 vb(i) = TMath::Pi()+i;
01861 TMatrixD vm(msize,1);
01862 TMatrixDColumn(vm,0) = vb;
01863 TMatrixDSym ms = THilbertMatrixDSym(0,msize-1);
01864 vb *= ms;
01865 ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
01866 TMatrixD mvm(ms,TMatrixD::kMult,vm);
01867 TMatrixD mvb(msize,1);
01868 TMatrixDColumn(mvb,0) = vb;
01869 ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
01870 }
01871
01872 #ifndef __CINT__
01873 if (nr >= Int_t(1.0e+5/msize/msize)) {
01874 #else
01875 if (nr >= Int_t(1.0e+3/msize/msize)) {
01876 #endif
01877 nr = 0;
01878 iloop--;
01879 } else
01880 nr++;
01881
01882 if (!ok) {
01883 if (gVerbose)
01884 cout << "Vector Matrix Multiplications failed for size " << msize << endl;
01885 break;
01886 }
01887 }
01888
01889 if (gVerbose)
01890 cout << "\nDone\n" << endl;
01891
01892 StatusPrint(12,"Matrix Vector Multiplications",ok);
01893 }
01894
01895
01896
01897
01898
01899 void mstress_inversion()
01900 {
01901 Bool_t ok = kTRUE;
01902
01903 Int_t iloop = gNrLoop;
01904 Int_t nr = 0;
01905 while (iloop >= 0) {
01906 const Int_t msize = gSizeA[iloop];
01907
01908 const Double_t epsilon = EPSILON*msize/10;
01909
01910 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
01911
01912 if (verbose)
01913 cout << "\n---> Verify matrix inversion for square matrices of size " << msize << endl;
01914 {
01915 if (verbose)
01916 cout << "\nTest inversion of a diagonal matrix" << endl;
01917 TMatrixD m(-1,msize,-1,msize);
01918 TMatrixD mi(TMatrixD::kZero,m);
01919 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
01920 m(i,i)=i-m.GetRowLwb()+1;
01921 mi(i,i) = 1/m(i,i);
01922 }
01923 TMatrixD mi1(TMatrixD::kInverted,m);
01924 m.Invert();
01925 ok &= VerifyMatrixIdentity(m,mi,verbose,epsilon);
01926 ok &= VerifyMatrixIdentity(mi1,mi,verbose,epsilon);
01927 }
01928
01929 if (ok)
01930 {
01931 if (verbose)
01932 cout << "Test inversion of an orthonormal (Haar) matrix" << endl;
01933 const Int_t order = Int_t(TMath::Log(msize)/TMath::Log(2));
01934 TMatrixD mr = THaarMatrixD(order);
01935 TMatrixD morig = mr;
01936 TMatrixD mt(TMatrixD::kTransposed,mr);
01937 Double_t det = -1;
01938 mr.Invert(&det);
01939 ok &= VerifyMatrixIdentity(mr,mt,verbose,epsilon);
01940 ok &= ( TMath::Abs(det-1) <= msize*epsilon ) ? kTRUE : kFALSE;
01941 if (verbose) {
01942 cout << "det = " << det << " deviation= " << TMath::Abs(det-1);
01943 cout << ( (TMath::Abs(det-1) < msize*epsilon) ? " OK" : " too large") <<endl;
01944 }
01945 TMatrixD mti(TMatrixD::kInverted,mt);
01946 ok &= VerifyMatrixIdentity(mti,morig,verbose,msize*epsilon);
01947 }
01948
01949 if (ok)
01950 {
01951 if (verbose)
01952 cout << "Test inversion of a good matrix with diagonal dominance" << endl;
01953 TMatrixD mh = THilbertMatrixD(msize,msize);
01954 TMatrixDDiag(mh) += 1;
01955 TMatrixD morig = mh;
01956 Double_t det_inv = 0;
01957 const Double_t det_comp = mh.Determinant();
01958 mh.Invert(&det_inv);
01959 if (verbose) {
01960 cout << "\tcomputed determinant " << det_comp << endl;
01961 cout << "\tdeterminant returned by Invert() " << det_inv << endl;
01962 }
01963
01964 if (verbose)
01965 cout << "\tcheck to see M^(-1) * M is E" << endl;
01966 TMatrixD mim(mh,TMatrixD::kMult,morig);
01967 TMatrixD unit(TMatrixD::kUnit,mh);
01968 ok &= VerifyMatrixIdentity(mim,unit,verbose,epsilon);
01969
01970 if (verbose)
01971 cout << "Test inversion through the matrix decompositions" << endl;
01972
01973 TMatrixDSym ms = THilbertMatrixDSym(msize);
01974 TMatrixDDiag(ms) += 1;
01975 if (verbose)
01976 cout << "Test inversion through SVD" << endl;
01977 TMatrixD inv_svd (msize,msize); TDecompSVD svd (ms); svd.Invert(inv_svd);
01978 ok &= VerifyMatrixIdentity(inv_svd,mh,verbose,epsilon);
01979 if (verbose)
01980 cout << "Test inversion through LU" << endl;
01981 TMatrixD inv_lu (msize,msize); TDecompLU lu (ms); lu.Invert(inv_lu);
01982 ok &= VerifyMatrixIdentity(inv_lu,mh,verbose,epsilon);
01983 if (verbose)
01984 cout << "Test inversion through Cholesky" << endl;
01985 TMatrixDSym inv_chol(msize); TDecompChol chol(ms); chol.Invert(inv_chol);
01986 ok &= VerifyMatrixIdentity(inv_chol,mh,verbose,epsilon);
01987 if (verbose)
01988 cout << "Test inversion through QRH" << endl;
01989 TMatrixD inv_qrh (msize,msize); TDecompQRH qrh (ms); qrh.Invert(inv_qrh);
01990 ok &= VerifyMatrixIdentity(inv_qrh,mh,verbose,epsilon);
01991 if (verbose)
01992 cout << "Test inversion through Bunch-Kaufman" << endl;
01993 TMatrixDSym inv_bk(msize); TDecompBK bk(ms); chol.Invert(inv_bk);
01994 ok &= VerifyMatrixIdentity(inv_bk,mh,verbose,epsilon);
01995
01996 if (verbose)
01997 cout << "\tcheck to see M * M^(-1) is E" << endl;
01998 TMatrixD mmi = morig; mmi *= mh;
01999 ok &= VerifyMatrixIdentity(mmi,unit,verbose,epsilon);
02000 }
02001
02002 #ifndef __CINT__
02003 if (nr >= Int_t(1.0e+5/msize/msize)) {
02004 #else
02005 if (nr >= Int_t(1.0e+3/msize/msize)) {
02006 #endif
02007 nr = 0;
02008 iloop--;
02009 } else
02010 nr++;
02011
02012 if (!ok) {
02013 if (gVerbose)
02014 cout << "Matrix Inversion failed for size " << msize << endl;
02015 break;
02016 }
02017 }
02018
02019 {
02020 if (gVerbose) {
02021 cout << "Check to see that Invert() and InvertFast() give identical results" <<endl;
02022 cout << " for size < (7x7)" <<endl;
02023 }
02024 Int_t size;
02025 for (size = 2; size < 7; size++) {
02026 TMatrixD m1 = THilbertMatrixD(size,size);
02027 m1(0,1) = TMath::Pi();
02028 TMatrixDDiag(m1) += 1;
02029 TMatrixD m2 = m1;
02030 Double_t det1 = 0.0;
02031 Double_t det2 = 0.0;
02032 m1.Invert(&det1);
02033 m2.InvertFast(&det2);
02034 ok &= VerifyMatrixIdentity(m1,m2,gVerbose,EPSILON);
02035 ok &= (TMath::Abs(det1-det2) < EPSILON);
02036 if (gVerbose) {
02037 cout << "det(Invert)= " << det1 << " det(InvertFast)= " << det2 <<endl;
02038 cout << " deviation= " << TMath::Abs(det1-det2);
02039 cout << ( (TMath::Abs(det1-det2) < EPSILON) ? " OK" : " too large") <<endl;
02040 }
02041 }
02042 for (size = 2; size < 7; size++) {
02043 TMatrixDSym ms1 = THilbertMatrixDSym(size);
02044 TMatrixDDiag(ms1) += 1;
02045 TMatrixDSym ms2 = ms1;
02046 Double_t det1 = 0.0;
02047 Double_t det2 = 0.0;
02048 ms1.Invert(&det1);
02049 ms2.InvertFast(&det2);
02050 ok &= VerifyMatrixIdentity(ms1,ms2,gVerbose,EPSILON);
02051 ok &= (TMath::Abs(det1-det2) < EPSILON);
02052 if (gVerbose) {
02053 cout << "det(Invert)= " << det1 << " det(InvertFast)= " << det2 <<endl;
02054 cout << " deviation= " << TMath::Abs(det1-det2);
02055 cout << ( (TMath::Abs(det1-det2) < EPSILON) ? " OK" : " too large") <<endl;
02056 }
02057 }
02058 }
02059
02060 if (gVerbose)
02061 cout << "\nDone\n" << endl;
02062
02063 StatusPrint(13,"Matrix Inversion",ok);
02064 }
02065
02066
02067
02068
02069
02070 void mstress_matrix_io()
02071 {
02072 if (gVerbose)
02073 cout << "\n---> Test matrix I/O" << endl;
02074
02075 Bool_t ok = kTRUE;
02076 const Double_t pattern = TMath::Pi();
02077
02078 TFile *f = new TFile("vmatrix.root", "RECREATE");
02079
02080 Char_t name[80];
02081 Int_t iloop = gNrLoop;
02082 while (iloop >= 0) {
02083 const Int_t msize = gSizeA[iloop];
02084
02085 const Int_t verbose = (gVerbose && iloop==gNrLoop);
02086
02087 TMatrixD m(msize,msize);
02088 m = pattern;
02089
02090 Double_t *pattern_array = new Double_t[msize*msize];
02091 for (Int_t i = 0; i < msize*msize; i++)
02092 pattern_array[i] = pattern;
02093 TMatrixD ma;
02094 ma.Use(msize,msize,pattern_array);
02095
02096 TMatrixDSym ms(msize);
02097 ms = pattern;
02098
02099 if (verbose)
02100 cout << "\nWrite matrix m to database" << endl;
02101 sprintf(name,"m_%d",msize);
02102 m.Write(name);
02103
02104 if (verbose)
02105 cout << "\nWrite matrix ma which adopts to database" << endl;
02106 sprintf(name,"ma_%d",msize);
02107 ma.Write(name);
02108
02109 if (verbose)
02110 cout << "\nWrite symmetric matrix ms to database" << endl;
02111 sprintf(name,"ms_%d",msize);
02112 ms.Write(name);
02113
02114 delete [] pattern_array;
02115
02116 iloop--;
02117 }
02118
02119 if (gVerbose)
02120 cout << "\nClose database" << endl;
02121 delete f;
02122
02123 if (gVerbose)
02124 cout << "\nOpen database in read-only mode and read matrix" << endl;
02125 TFile *f1 = new TFile("vmatrix.root");
02126
02127 iloop = gNrLoop;
02128 while (iloop >= 0) {
02129 const Int_t msize = gSizeA[iloop];
02130
02131 const Int_t verbose = (gVerbose && iloop==gNrLoop);
02132
02133 TMatrixD m(msize,msize);
02134 m = pattern;
02135 sprintf(name,"m_%d",msize);
02136 TMatrixD *mr = (TMatrixD*) f1->Get(name);
02137 sprintf(name,"ma_%d",msize);
02138 TMatrixD *mar = (TMatrixD*) f1->Get(name);
02139 sprintf(name,"ms_%d",msize);
02140 TMatrixDSym *msr = (TMatrixDSym*) f1->Get(name);
02141
02142 if (verbose)
02143 cout << "\nRead matrix should be same as original" << endl;
02144 ok &= ((*mr) == m) ? kTRUE : kFALSE;
02145 ok &= ((*mar) == m) ? kTRUE : kFALSE;
02146 ok &= ((*msr) == m) ? kTRUE : kFALSE;
02147
02148 iloop--;
02149 }
02150
02151 delete f1;
02152
02153 if (gVerbose)
02154 cout << "\nDone\n" << endl;
02155
02156 StatusPrint(14,"Matrix Persistence",ok);
02157 }
02158
02159
02160
02161
02162 void spstress_allocation(Int_t msize)
02163 {
02164 if (gVerbose)
02165 cout << "\n\n---> Test allocation and compatibility check" << endl;
02166
02167 Int_t i,j;
02168 Bool_t ok = kTRUE;
02169
02170 TMatrixDSparse m1(0,3,0,msize-1);
02171 {
02172 Int_t nr = 4*msize;
02173 Int_t *irow = new Int_t[nr];
02174 Int_t *icol = new Int_t[nr];
02175 Double_t *val = new Double_t[nr];
02176
02177 Int_t n = 0;
02178 for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
02179 for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
02180 irow[n] = i;
02181 icol[n] = j;
02182 val[n] = TMath::Pi()*i+TMath::E()*j;
02183 n++;
02184 }
02185 }
02186 m1.SetMatrixArray(nr,irow,icol,val);
02187 delete [] irow;
02188 delete [] icol;
02189 delete [] val;
02190 }
02191
02192 TMatrixDSparse m2(0,3,0,msize-1);
02193 TMatrixDSparse m3(1,4,0,msize-1);
02194 TMatrixDSparse m4(m1);
02195
02196 if (gVerbose) {
02197 cout << "\nStatus information reported for matrix m3:" << endl;
02198 cout << " Row lower bound ... " << m3.GetRowLwb() << endl;
02199 cout << " Row upper bound ... " << m3.GetRowUpb() << endl;
02200 cout << " Col lower bound ... " << m3.GetColLwb() << endl;
02201 cout << " Col upper bound ... " << m3.GetColUpb() << endl;
02202 cout << " No. rows ..........." << m3.GetNrows() << endl;
02203 cout << " No. cols ..........." << m3.GetNcols() << endl;
02204 cout << " No. of elements ...." << m3.GetNoElements() << endl;
02205 }
02206
02207 if (gVerbose)
02208 cout << "Check matrices 1 & 4 for compatibility" << endl;
02209 ok &= AreCompatible(m1,m4,gVerbose);
02210
02211 if (gVerbose)
02212 cout << "m2 has to be compatible with m3 after resizing to m3" << endl;
02213 m2.ResizeTo(m3);
02214 ok &= AreCompatible(m2,m3,gVerbose);
02215
02216 TMatrixD m5_d(m1.GetNrows()+1,m1.GetNcols()+5);
02217 for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
02218 for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
02219 m5_d(i,j) = TMath::Pi()*i+TMath::E()*j;
02220 TMatrixDSparse m5(m5_d);
02221
02222 if (gVerbose)
02223 cout << "m1 has to be compatible with m5 after resizing to m5" << endl;
02224 m1.ResizeTo(m5.GetNrows(),m5.GetNcols());
02225 ok &= AreCompatible(m1,m5,gVerbose);
02226
02227 if (gVerbose)
02228 cout << "m1 has to be equal to m4 after stretching and shrinking" << endl;
02229 m1.ResizeTo(m4.GetNrows(),m4.GetNcols());
02230 ok &= VerifyMatrixIdentity(m1,m4,gVerbose,EPSILON);
02231 if (gVerbose)
02232 cout << "m5 has to be equal to m1 after shrinking" << endl;
02233 m5.ResizeTo(m1.GetNrows(),m1.GetNcols());
02234 ok &= VerifyMatrixIdentity(m1,m5,gVerbose,msize*EPSILON);
02235
02236 if (gVerbose)
02237 cout << "\nDone\n" << endl;
02238
02239 StatusPrint(1,"Allocation, Resizing",ok);
02240 }
02241
02242
02243
02244
02245
02246 void spstress_matrix_fill(Int_t rsize,Int_t csize)
02247 {
02248 Bool_t ok = kTRUE;
02249
02250 if (csize < 4) {
02251 Error("spstress_matrix_fill","rsize should be >= 4");
02252 ok = kFALSE;
02253 StatusPrint(2,"Filling, Inserting, Using",ok);
02254 return;
02255 }
02256
02257 if (csize < 4) {
02258 Error("spstress_matrix_fill","csize should be >= 4");
02259 ok = kFALSE;
02260 StatusPrint(2,"Filling, Inserting, Using",ok);
02261 return;
02262 }
02263
02264 TMatrixD m_d(-1,rsize-2,1,csize);
02265 for (Int_t i = m_d.GetRowLwb(); i <= m_d.GetRowUpb(); i++)
02266 for (Int_t j = m_d.GetColLwb(); j <= m_d.GetColUpb(); j++)
02267
02268 if (j != 2)
02269 m_d(i,j) = TMath::Pi()*i+TMath::E()*j;
02270 TMatrixDSparse m(m_d);
02271
02272 {
02273 if (gVerbose)
02274 cout << "Check filling through operator(i,j) without setting sparse index" << endl;
02275 TMatrixDSparse m1(-1,rsize-2,1,csize);
02276
02277 for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
02278 for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
02279 if (j != 2)
02280 m1(i,j) = TMath::Pi()*i+TMath::E()*j;
02281 ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
02282 }
02283
02284 {
02285 if (gVerbose)
02286 cout << "Check filling through operator(i,j)" << endl;
02287 TMatrixDSparse m2(-1,rsize-2,1,csize);
02288 m2.SetSparseIndex(m);
02289
02290 for (Int_t i = m2.GetRowLwb(); i <= m2.GetRowUpb(); i++)
02291 for (Int_t j = m2.GetColLwb(); j <= m2.GetColUpb(); j++)
02292 if (j != 2)
02293 m2(i,j) = TMath::Pi()*i+TMath::E()*j;
02294 ok &= VerifyMatrixIdentity(m2,m,gVerbose,EPSILON);
02295 }
02296
02297 {
02298 if (gVerbose)
02299 cout << "Check insertion/extraction of sub-matrices" << endl;
02300 {
02301 TMatrixDSparse m_sub1 = m;
02302 m_sub1.ResizeTo(0,rsize-2,2,csize);
02303 TMatrixDSparse m_sub2 = m.GetSub(0,rsize-2,2,csize,"");
02304 ok &= VerifyMatrixIdentity(m_sub1,m_sub2,gVerbose,EPSILON);
02305 }
02306
02307 {
02308 TMatrixDSparse m3(-1,rsize-2,1,csize);
02309 TMatrixDSparse m_part1 = m.GetSub(-1,rsize-2,1,csize,"");
02310 m3.SetSub(-1,1,m_part1);
02311 ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
02312 }
02313
02314 {
02315 TMatrixDSparse m4(-1,rsize-2,1,csize);
02316 TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"");
02317 TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"");
02318 TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"");
02319 TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"");
02320 m4.SetSub(0,2,m_part1);
02321 m4.SetSub(0,1,m_part2);
02322 m4.SetSub(-1,2,m_part3);
02323 m4.SetSub(-1,1,m_part4);
02324 ok &= VerifyMatrixIdentity(m,m4,gVerbose,EPSILON);
02325 }
02326
02327 {
02328
02329 TMatrixDSparse m5(-1,rsize-2,1,csize);
02330 TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"");
02331 TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"");
02332 TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"");
02333 TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"");
02334 m5.SetSub(-1,1,m_part4);
02335 m5.SetSub(-1,2,m_part3);
02336 m5.SetSub(0,1,m_part2);
02337 m5.SetSub(0,2,m_part1);
02338 ok &= VerifyMatrixIdentity(m,m5,gVerbose,EPSILON);
02339 }
02340
02341 {
02342 TMatrixDSparse m6(-1,rsize-2,1,csize);
02343 TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"S");
02344 TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"S");
02345 TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"S");
02346 TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"S");
02347 m6.SetSub(0,2,m_part1);
02348 m6.SetSub(0,1,m_part2);
02349 m6.SetSub(-1,2,m_part3);
02350 m6.SetSub(-1,1,m_part4);
02351 ok &= VerifyMatrixIdentity(m,m6,gVerbose,EPSILON);
02352 }
02353 }
02354
02355 {
02356 if (gVerbose)
02357 cout << "Check insertion/extraction of rows" << endl;
02358
02359 TMatrixDSparse m1 = m;
02360 TVectorD v1(1,csize);
02361 TVectorD v2(1,csize);
02362 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
02363 v1 = TMatrixDSparseRow(m,i);
02364 m1.InsertRow(i,1,v1.GetMatrixArray());
02365 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
02366 m1.InsertRow(i,3,v1.GetMatrixArray()+2,v1.GetNrows()-2);
02367 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
02368
02369 m1.ExtractRow(i,1,v2.GetMatrixArray());
02370 ok &= VerifyVectorIdentity(v1,v2,gVerbose,EPSILON);
02371 m1.ExtractRow(i,3,v2.GetMatrixArray()+2,v1.GetNrows()-2);
02372 ok &= VerifyVectorIdentity(v1,v2,gVerbose,EPSILON);
02373 }
02374 }
02375
02376 {
02377 if (gVerbose)
02378 cout << "Check array Use" << endl;
02379 {
02380 TMatrixDSparse *m1a = new TMatrixDSparse(m);
02381 TMatrixDSparse *m2a = new TMatrixDSparse();
02382 m2a->Use(*m1a);
02383 m2a->Sqr();
02384 TMatrixDSparse m7 = m; m7.Sqr();
02385 ok &= VerifyMatrixIdentity(m7,*m1a,gVerbose,EPSILON);
02386 delete m1a;
02387 delete m2a;
02388 }
02389 }
02390
02391 if (gVerbose)
02392 cout << "\nDone\n" << endl;
02393
02394 StatusPrint(2,"Filling, Inserting, Using",ok);
02395 }
02396
02397
02398
02399
02400
02401 void spstress_element_op(Int_t rsize,Int_t csize)
02402 {
02403 Bool_t ok = kTRUE;
02404 const Double_t pattern = 8.625;
02405
02406 TMatrixDSparse m(-1,rsize-2,1,csize);
02407
02408 if (gVerbose)
02409 cout << "Creating zero m1 ..." << endl;
02410 TMatrixDSparse m1(TMatrixDSparse::kZero, m);
02411 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
02412
02413 if (gVerbose)
02414 cout << "Comparing m1 with 0 ..." << endl;
02415 R__ASSERT(m1 == 0);
02416 R__ASSERT(!(m1 != 0));
02417
02418 if (gVerbose)
02419 cout << "Writing a pattern " << pattern << " by assigning through SetMatrixArray..." << endl;
02420 {
02421 const Int_t nr = rsize*csize;
02422 Int_t *irow = new Int_t[nr];
02423 Int_t *icol = new Int_t[nr];
02424 Double_t *val = new Double_t[nr];
02425
02426 Int_t n = 0;
02427 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
02428 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
02429 irow[n] = i;
02430 icol[n] = j;
02431 val[n] = pattern;
02432 n++;
02433 }
02434 }
02435 m.SetMatrixArray(nr,irow,icol,val);
02436 delete [] irow;
02437 delete [] icol;
02438 delete [] val;
02439 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
02440 }
02441
02442 if (gVerbose)
02443 cout << "Writing the pattern by assigning to m1 as a whole ..." << endl;
02444 m1.SetSparseIndex(m);
02445 m1 = pattern;
02446 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
02447
02448 if (gVerbose)
02449 cout << "Comparing m and m1 ..." << endl;
02450 R__ASSERT(m == m1);
02451 if (gVerbose)
02452 cout << "Comparing (m=0) and m1 ..." << endl;
02453 R__ASSERT(!((m=0) == m1));
02454
02455 if (gVerbose)
02456 cout << "Clearing m1 ..." << endl;
02457 m1.Zero();
02458 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
02459
02460 if (gVerbose)
02461 cout << "\nSet m = pattern" << endl;
02462 m = pattern;
02463 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
02464 if (gVerbose)
02465 cout << " add the doubled pattern with the negative sign" << endl;
02466 m += -2*pattern;
02467 ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
02468 if (gVerbose)
02469 cout << " subtract the trippled pattern with the negative sign" << endl;
02470 m -= -3*pattern;
02471 ok &= VerifyMatrixValue(m,2*pattern,gVerbose,EPSILON);
02472
02473 if (gVerbose)
02474 cout << "\nVerify comparison operations when all elems are the same" << endl;
02475 m = pattern;
02476 R__ASSERT( m == pattern && !(m != pattern) );
02477 R__ASSERT( m > 0 && m >= pattern && m <= pattern );
02478 R__ASSERT( m > -pattern && m >= -pattern );
02479 R__ASSERT( m <= pattern && !(m < pattern) );
02480 m -= 2*pattern;
02481 R__ASSERT( m < -pattern/2 && m <= -pattern/2 );
02482 R__ASSERT( m >= -pattern && !(m > -pattern) );
02483
02484 if (gVerbose)
02485 cout << "\nVerify comparison operations when not all elems are the same" << endl;
02486 {
02487 Int_t nr = rsize*csize;
02488 Int_t *irow = new Int_t[nr];
02489 Int_t *icol = new Int_t[nr];
02490 Double_t *val = new Double_t[nr];
02491
02492 Int_t n = 0;
02493 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
02494 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
02495 irow[n] = i;
02496 icol[n] = j;
02497 val[n] = pattern;
02498 n++;
02499 }
02500 }
02501 val[n-1] = pattern-1;
02502 m.SetMatrixArray(nr,irow,icol,val);
02503 delete [] irow;
02504 delete [] icol;
02505 delete [] val;
02506 }
02507
02508 R__ASSERT( !(m == pattern) && !(m != pattern) );
02509 R__ASSERT( m != 0 );
02510 R__ASSERT( !(m >= pattern) && m <= pattern && !(m<pattern) );
02511 R__ASSERT( !(m <= pattern-1) && m >= pattern-1 && !(m>pattern-1) );
02512 if (gVerbose)
02513 cout << "\nAssign 2*pattern to m by repeating additions" << endl;
02514 m = 0; m += pattern; m += pattern;
02515 if (gVerbose)
02516 cout << "Assign 2*pattern to m1 by multiplying by two " << endl;
02517 m1.SetSparseIndex(m);
02518 m1 = pattern; m1 *= 2;
02519 ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
02520 R__ASSERT( m == m1 );
02521 if (gVerbose)
02522 cout << "Multiply m1 by one half returning it to the 1*pattern" << endl;
02523 m1 *= 1/2.;
02524 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
02525
02526 if (gVerbose)
02527 cout << "\nAssign -pattern to m and m1" << endl;
02528 m = 0; m -= pattern; m1 = -pattern;
02529 ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
02530 R__ASSERT( m == m1 );
02531 if (gVerbose)
02532 cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same" << endl;
02533 m.Sqr();
02534 ok &= VerifyMatrixValue(m,pattern*pattern,gVerbose,EPSILON);
02535 m.Sqrt();
02536 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
02537 m1.Abs();
02538 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
02539 ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
02540
02541 if (gVerbose)
02542 cout << "\nDone\n" << endl;
02543
02544 StatusPrint(3,"Uniform matrix operations",ok);
02545 }
02546
02547
02548
02549
02550
02551 void spstress_binary_ebe_op(Int_t rsize, Int_t csize)
02552 {
02553 if (gVerbose)
02554 cout << "\n---> Test Binary Matrix element-by-element operations" << endl;
02555
02556 Bool_t ok = kTRUE;
02557
02558 const Double_t pattern = 4.25;
02559
02560 TMatrixD m_d(2,rsize+1,0,csize-1); m_d = 1;
02561 TMatrixDSparse m (TMatrixDSparse::kZero,m_d); m.SetSparseIndex (m_d);
02562 TMatrixDSparse m1(TMatrixDSparse::kZero,m); m1.SetSparseIndex(m_d);
02563
02564 TMatrixDSparse mp(TMatrixDSparse::kZero,m1); mp.SetSparseIndex(m_d);
02565 {
02566 for (Int_t i = mp.GetRowLwb(); i <= mp.GetRowUpb(); i++)
02567 for (Int_t j = mp.GetColLwb(); j <= mp.GetColUpb(); j++)
02568 mp(i,j) = TMath::Pi()*i+TMath::E()*(j+1);
02569 }
02570
02571 if (gVerbose)
02572 cout << "\nVerify assignment of a matrix to the matrix" << endl;
02573 m = pattern;
02574 m1 = m;
02575 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
02576 R__ASSERT( m1 == m );
02577
02578 if (gVerbose)
02579 cout << "\nAdding the matrix to itself, uniform pattern " << pattern << endl;
02580 m.Zero(); m.SetSparseIndex(m_d); m = pattern;
02581
02582 m1 = m; m1 += m1;
02583 ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
02584 if (gVerbose)
02585 cout << " subtracting two matrices ..." << endl;
02586 m1 -= m;
02587 ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
02588 if (gVerbose)
02589 cout << " subtracting the matrix from itself" << endl;
02590 m1 -= m1;
02591 ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
02592 m1.SetSparseIndex(m_d);
02593
02594 if (gVerbose) {
02595 cout << "\nArithmetic operations on matrices with not the same elements" << endl;
02596 cout << " adding mp to the zero matrix..." << endl;
02597 }
02598 m.Zero(); m.SetSparseIndex(m_d); m += mp;
02599 ok &= VerifyMatrixIdentity(m,mp,gVerbose,EPSILON);
02600 if (gVerbose)
02601 cout << " making m = 3*mp and m1 = 3*mp, via add() and succesive mult" << endl;
02602 m1 = m;
02603 Add(m,2.,mp);
02604 m1 += m1; m1 += mp;
02605 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
02606 if (gVerbose)
02607 cout << " clear both m and m1, by subtracting from itself and via add()" << endl;
02608 m1 -= m1;
02609 Add(m,-3.,mp);
02610 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
02611
02612 if (gVerbose) {
02613 cout << "\nTesting element-by-element multiplications and divisions" << endl;
02614 cout << " squaring each element with sqr() and via multiplication" << endl;
02615 }
02616 m.SetSparseIndex(m_d); m = mp;
02617 m1.SetSparseIndex(m_d); m1 = mp;
02618 m.Sqr();
02619 ElementMult(m1,m1);
02620 ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
02621 if (gVerbose)
02622 cout << " compare (m = pattern^2)/pattern with pattern" << endl;
02623 m = pattern; m1 = pattern;
02624 m.Sqr();
02625 ElementDiv(m,m1);
02626 ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
02627 if (gVerbose)
02628 Compare(m1,m);
02629
02630 if (gVerbose)
02631 cout << "\nDone\n" << endl;
02632
02633 StatusPrint(4,"Binary Matrix element-by-element operations",ok);
02634 }
02635
02636
02637
02638
02639
02640 void spstress_transposition(Int_t msize)
02641 {
02642 if (gVerbose) {
02643 cout << "\n---> Verify matrix transpose "
02644 "for matrices of a characteristic size " << msize << endl;
02645 }
02646
02647 Bool_t ok = kTRUE;
02648 {
02649 if (gVerbose)
02650 cout << "\nCheck to see that a square UnitMatrix stays the same";
02651 TMatrixDSparse m(msize,msize);
02652 m.UnitMatrix();
02653 TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
02654 ok &= ( m == mt ) ? kTRUE : kFALSE;
02655 }
02656
02657 {
02658 if (gVerbose)
02659 cout << "\nTest a non-square UnitMatrix";
02660 TMatrixDSparse m(msize,msize+1);
02661 m.UnitMatrix();
02662 TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
02663 R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows() );
02664 const Int_t rowlwb = m.GetRowLwb();
02665 const Int_t collwb = m.GetColLwb();
02666 const Int_t upb = TMath::Min(m.GetRowUpb(),m.GetColUpb());
02667 TMatrixDSparse m_part = m.GetSub(rowlwb,upb,collwb,upb);
02668 TMatrixDSparse mt_part = mt.GetSub(rowlwb,upb,collwb,upb);
02669 ok &= ( m_part == mt_part ) ? kTRUE : kFALSE;
02670 }
02671
02672 {
02673 if (gVerbose)
02674 cout << "\nCheck to see that a symmetric (Hilbert)Matrix stays the same";
02675 TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize,msize));
02676 TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
02677 ok &= ( m == mt ) ? kTRUE : kFALSE;
02678 }
02679
02680 {
02681 if (gVerbose)
02682 cout << "\nCheck transposing a non-symmetric matrix";
02683 TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize+1,msize));
02684 m(1,2) = TMath::Pi();
02685 TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
02686 R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows());
02687 R__ASSERT(mt(2,1) == (Double_t)TMath::Pi() && mt(1,2) != (Double_t)TMath::Pi());
02688 R__ASSERT(mt[2][1] == (Double_t)TMath::Pi() && mt[1][2] != (Double_t)TMath::Pi());
02689
02690 if (gVerbose)
02691 cout << "\nCheck double transposing a non-symmetric matrix" << endl;
02692 TMatrixDSparse mtt(TMatrixDSparse::kTransposed,mt);
02693 ok &= ( m == mtt ) ? kTRUE : kFALSE;
02694 }
02695
02696 if (gVerbose)
02697 cout << "\nDone\n" << endl;
02698
02699 StatusPrint(5,"Matrix transposition",ok);
02700 }
02701
02702
02703
02704
02705
02706 void spstress_norms(Int_t rsize_req,Int_t csize_req)
02707 {
02708 if (gVerbose)
02709 cout << "\n---> Verify norm calculations" << endl;
02710
02711 Bool_t ok = kTRUE;
02712 const Double_t pattern = 10.25;
02713
02714 Int_t rsize = rsize_req;
02715 if (rsize%2 != 0) rsize--;
02716 Int_t csize = csize_req;
02717 if (csize%2 != 0) csize--;
02718 if (rsize%2 == 1 || csize%2 == 1) {
02719 cout << "rsize: " << rsize <<endl;
02720 cout << "csize: " << csize <<endl;
02721 Fatal("spstress_norms","Sorry, size of the matrix to test must be even for this test\n");
02722 }
02723
02724 TMatrixD m_d(rsize,csize); m_d = 1;
02725 TMatrixDSparse m(rsize,csize); m.SetSparseIndex(m_d);
02726
02727 if (gVerbose)
02728 cout << "\nAssign " << pattern << " to all the elements and check norms" << endl;
02729 m = pattern;
02730 if (gVerbose)
02731 cout << " 1. (col) norm should be pattern*nrows" << endl;
02732 ok &= ( m.Norm1() == pattern*m.GetNrows() ) ? kTRUE : kFALSE;
02733 ok &= ( m.Norm1() == m.ColNorm() ) ? kTRUE : kFALSE;
02734 if (gVerbose)
02735 cout << " Inf (row) norm should be pattern*ncols" << endl;
02736 ok &= ( m.NormInf() == pattern*m.GetNcols() ) ? kTRUE : kFALSE;
02737 ok &= ( m.NormInf() == m.RowNorm() ) ? kTRUE : kFALSE;
02738 if (gVerbose)
02739 cout << " Square of the Eucl norm has got to be pattern^2 * no_elems" << endl;
02740 ok &= ( m.E2Norm() == (pattern*pattern)*m.GetNoElements() ) ? kTRUE : kFALSE;
02741 TMatrixDSparse m1(m); m1 = 1;
02742 ok &= ( m.E2Norm() == E2Norm(m+1.,m1) ) ? kTRUE : kFALSE;
02743
02744 if (gVerbose)
02745 cout << "\nDone\n" << endl;
02746
02747 StatusPrint(6,"Matrix Norms",ok);
02748 }
02749
02750
02751
02752
02753
02754 void spstress_mm_multiplications()
02755 {
02756 Bool_t ok = kTRUE;
02757 Int_t i;
02758
02759 Int_t iloop = 0;
02760 Int_t nr = 0;
02761 while (iloop <= gNrLoop) {
02762 const Int_t msize = gSizeA[iloop];
02763 const Double_t epsilon = EPSILON*msize/100;
02764
02765 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
02766
02767 if (msize <= 5) {
02768 iloop++;
02769 continue;
02770 }
02771
02772 if (verbose)
02773 cout << "\n---> Verify matrix multiplications "
02774 "for matrices of the characteristic size " << msize << endl;
02775
02776 {
02777 if (verbose)
02778 cout << "\nTest inline multiplications of the UnitMatrix" << endl;
02779 if (ok)
02780 {
02781 TMatrixDSparse m = TMatrixD(THilbertMatrixD(-1,msize,-1,msize));
02782 TMatrixDSparse ur(TMatrixDSparse::kUnit,m);
02783 TMatrixDSparse ul(TMatrixDSparse::kUnit,m);
02784 m(3,1) = TMath::Pi();
02785 ul *= m;
02786 m *= ur;
02787 ok &= VerifyMatrixIdentity(ul,m,verbose,epsilon);
02788 }
02789
02790 if (ok)
02791 {
02792 TMatrixD m_d = THilbertMatrixD(-1,msize,-1,msize);
02793 TMatrixDSparse ur(TMatrixDSparse::kUnit,m_d);
02794 TMatrixDSparse ul(TMatrixDSparse::kUnit,m_d);
02795 m_d(3,1) = TMath::Pi();
02796 ul *= m_d;
02797 m_d *= TMatrixD(ur);
02798 ok &= VerifyMatrixIdentity(ul,m_d,verbose,epsilon);
02799 }
02800 }
02801
02802 if (ok)
02803 {
02804 if (verbose)
02805 cout << "Test XPP = X where P is a permutation matrix" << endl;
02806 TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize-1,msize));
02807 m(2,3) = TMath::Pi();
02808 TMatrixDSparse eth = m;
02809 TMatrixDSparse p(msize,msize);
02810 {
02811 Int_t *irow = new Int_t[msize];
02812 Int_t *icol = new Int_t[msize];
02813 Double_t *val = new Double_t[msize];
02814
02815 Int_t n = 0;
02816 for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++) {
02817 irow[n] = p.GetRowUpb()+p.GetRowLwb()-i;
02818 icol[n] = i;
02819 val[n] = 1;
02820 n++;
02821 }
02822 p.SetMatrixArray(msize,irow,icol,val);
02823 delete [] irow;
02824 delete [] icol;
02825 delete [] val;
02826 }
02827 m *= p;
02828 m *= p;
02829 ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
02830 }
02831
02832 if (ok)
02833 {
02834 if (verbose)
02835 cout << "Test general matrix multiplication through inline mult" << endl;
02836 TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize-2,msize));
02837 m(3,3) = TMath::Pi();
02838 TMatrixD p_d = THilbertMatrixD(msize,msize);
02839 TMatrixDDiag(p_d) += 1;
02840 TMatrixDSparse mp(m,TMatrixDSparse::kMult,p_d);
02841 TMatrixDSparse m1 = m;
02842 m *= p_d;
02843 ok &= VerifyMatrixIdentity(m,mp,verbose,epsilon);
02844 TMatrixDSparse pt_d(TMatrixDSparse::kTransposed,p_d);
02845 TMatrixDSparse mp1(m1,TMatrixDSparse::kMultTranspose,pt_d);
02846 VerifyMatrixIdentity(m,mp1,verbose,epsilon);
02847
02848 ok &= ( !(m1 == m) );
02849 TMatrixDSparse mp2(TMatrixDSparse::kZero,m1);
02850 ok &= ( mp2 == 0 );
02851 mp2.SetSparseIndex(m1);
02852 mp2.Mult(m1,p_d);
02853 ok &= VerifyMatrixIdentity(m,mp2,verbose,epsilon);
02854
02855 if (verbose)
02856 cout << "Test XP=X*P vs XP=X;XP*=P" << endl;
02857 TMatrixDSparse mp3 = m1*p_d;
02858 ok &= VerifyMatrixIdentity(m,mp3,verbose,epsilon);
02859 }
02860
02861 #ifndef __CINT__
02862 if (nr >= Int_t(1.0e+5/msize/msize)) {
02863 #else
02864 if (nr >= Int_t(1.0e+3/msize/msize)) {
02865 #endif
02866 nr = 0;
02867 iloop++;
02868 } else
02869 nr++;
02870
02871 if (!ok) {
02872 if (gVerbose)
02873 cout << "General Matrix Multiplications failed for size " << msize << endl;
02874 break;
02875 }
02876 }
02877
02878 if (ok)
02879 {
02880 if (gVerbose)
02881 cout << "Check to see UU' = U'U = E when U is the Haar matrix" << endl;
02882 const Int_t order = 5;
02883 const Int_t no_sub_cols = (1<<order)-order;
02884 TMatrixDSparse haar_sub = TMatrixD(THaarMatrixD(order,no_sub_cols));
02885 TMatrixDSparse haar_sub_t(TMatrixDSparse::kTransposed,haar_sub);
02886 TMatrixDSparse hsths(haar_sub_t,TMatrixDSparse::kMult,haar_sub);
02887
02888 TMatrixDSparse square(no_sub_cols,no_sub_cols);
02889 TMatrixDSparse units(TMatrixDSparse::kUnit,square);
02890 ok &= ( hsths.GetNrows() == no_sub_cols && hsths.GetNcols() == no_sub_cols );
02891 ok &= VerifyMatrixIdentity(hsths,units,gVerbose,EPSILON);
02892
02893 TMatrixDSparse haar = TMatrixD(THaarMatrixD(order));
02894 TMatrixDSparse haar_t(TMatrixDSparse::kTransposed,haar);
02895 TMatrixDSparse unit(TMatrixDSparse::kUnit,haar);
02896
02897 TMatrixDSparse hth(haar_t,TMatrixDSparse::kMult,haar);
02898 ok &= VerifyMatrixIdentity(hth,unit,gVerbose,EPSILON);
02899
02900 TMatrixDSparse hht(haar,TMatrixDSparse::kMultTranspose,haar);
02901 ok &= VerifyMatrixIdentity(hht,unit,gVerbose,EPSILON);
02902
02903 TMatrixDSparse hht1 = haar; hht1 *= haar_t;
02904 ok &= VerifyMatrixIdentity(hht1,unit,gVerbose,EPSILON);
02905
02906 TMatrixDSparse hht2(TMatrixDSparse::kZero,haar);
02907 hht2.SetSparseIndex(hht1);
02908 hht2.Mult(haar,haar_t);
02909 ok &= VerifyMatrixIdentity(hht2,unit,gVerbose,EPSILON);
02910 }
02911
02912 if (gVerbose)
02913 cout << "\nDone\n" << endl;
02914
02915 StatusPrint(7,"General Matrix Multiplications",ok);
02916 }
02917
02918
02919
02920
02921
02922 void spstress_vm_multiplications()
02923 {
02924 Bool_t ok = kTRUE;
02925
02926 Int_t iloop = gNrLoop;
02927 Int_t nr = 0;
02928 while (iloop >= 0) {
02929 const Int_t msize = gSizeA[iloop];
02930
02931 const Double_t epsilon = EPSILON*msize;
02932
02933 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
02934
02935 if (verbose)
02936 cout << "\n---> Verify vector-matrix multiplications "
02937 "for matrices of the characteristic size " << msize << endl;
02938
02939 {
02940 if (verbose)
02941 cout << "\nCheck shrinking a vector by multiplying by a non-sq unit matrix" << endl;
02942 TVectorD vb(-2,msize);
02943 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
02944 vb(i) = TMath::Pi()-i;
02945 ok &= ( vb != 0 ) ? kTRUE : kFALSE;
02946 TMatrixDSparse mc(1,msize-2,-2,msize);
02947 mc.UnitMatrix();
02948 TVectorD v1 = vb;
02949 TVectorD v2 = vb;
02950 v1 *= mc;
02951 v2.ResizeTo(1,msize-2);
02952 ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
02953 }
02954
02955 if (ok)
02956 {
02957 if (verbose)
02958 cout << "Check expanding a vector by multiplying by a non-sq unit matrix" << endl;
02959 TVectorD vb(msize);
02960 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
02961 vb(i) = TMath::Pi()+i;
02962 ok &= ( vb != 0 ) ? kTRUE : kFALSE;
02963 TMatrixDSparse me(2,msize+5,0,msize-1);
02964 me.UnitMatrix();
02965 TVectorD v1 = vb;
02966 TVectorD v2 = vb;
02967 v1 *= me;
02968 v2.ResizeTo(v1);
02969 ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
02970 }
02971
02972 if (ok)
02973 {
02974 if (verbose)
02975 cout << "Check general matrix-vector multiplication" << endl;
02976 TVectorD vb(msize);
02977 for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
02978 vb(i) = TMath::Pi()+i;
02979 TMatrixD vm(msize,1);
02980 TMatrixDColumn(vm,0) = vb;
02981 TMatrixD hilbert_with_zeros = THilbertMatrixD(0,msize,0,msize-1);
02982 TMatrixDRow (hilbert_with_zeros,3) = 0.0;
02983 TMatrixDColumn(hilbert_with_zeros,3) = 0.0;
02984 const TMatrixDSparse m = hilbert_with_zeros;
02985 vb *= m;
02986 ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
02987 TMatrixDSparse mvm(m,TMatrixDSparse::kMult,vm);
02988 TMatrixD mvb(msize+1,1);
02989 TMatrixDColumn(mvb,0) = vb;
02990 ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
02991 }
02992
02993 #ifndef __CINT__
02994 if (nr >= Int_t(1.0e+5/msize/msize)) {
02995 #else
02996 if (nr >= Int_t(1.0e+3/msize/msize)) {
02997 #endif
02998 nr = 0;
02999 iloop--;
03000 } else
03001 nr++;
03002
03003 if (!ok) {
03004 if (gVerbose)
03005 cout << "Vector Matrix Multiplications failed for size " << msize << endl;
03006 break;
03007 }
03008 }
03009
03010 if (gVerbose)
03011 cout << "\nDone\n" << endl;
03012
03013 StatusPrint(8,"Matrix Vector Multiplications",ok);
03014 }
03015
03016
03017
03018
03019
03020 void spstress_matrix_slices(Int_t vsize)
03021 {
03022 if (gVerbose)
03023 cout << "\n---> Test operations with vectors and sparse matrix slices" << endl;
03024
03025 Bool_t ok = kTRUE;
03026 const Double_t pattern = 8.625;
03027
03028 TVectorD vc(0,vsize);
03029 TVectorD vr(0,vsize+1);
03030 TMatrixD m_d(0,vsize,0,vsize+1); m_d = pattern;
03031 TMatrixDSparse m(m_d);
03032
03033 Int_t i,j;
03034 if (gVerbose)
03035 cout << "\nCheck modifying the matrix row-by-row" << endl;
03036 m = pattern;
03037 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03038 for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
03039 TMatrixDSparseRow(m,i) = pattern+2;
03040 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03041 vr = TMatrixDSparseRow(m,i);
03042 ok &= VerifyVectorValue(vr,pattern+2,gVerbose,EPSILON);
03043 vr = TMatrixDSparseRow(m,i+1 > m.GetRowUpb() ? m.GetRowLwb() : i+1);
03044 ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
03045 TMatrixDSparseRow(m,i) += -2;
03046 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03047 }
03048
03049 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03050 for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
03051 vr = pattern-2;
03052 TMatrixDSparseRow(m,i) = vr;
03053 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03054 {
03055 TMatrixDSparseRow mr(m,i);
03056 for (j = m.GetColLwb(); j <= m.GetColUpb(); j++)
03057 mr(j) *= 8;
03058 }
03059 vr = TMatrixDSparseRow(m,i);
03060 ok &= VerifyVectorValue(vr,8*(pattern-2),gVerbose,EPSILON);
03061 TMatrixDSparseRow(m,i) *= 1./8;
03062 TMatrixDSparseRow(m,i) += 2;
03063 vr = TMatrixDSparseRow(m,i);
03064 ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
03065 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03066 }
03067
03068 if (gVerbose)
03069 cout << "\nCheck modifying the matrix diagonal" << endl;
03070 m = pattern;
03071 TMatrixDSparseDiag td = m;
03072 td = pattern-3;
03073 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03074 vc = TMatrixDSparseDiag(m);
03075 ok &= VerifyVectorValue(vc,pattern-3,gVerbose,EPSILON);
03076 td += 3;
03077 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03078 vc = pattern+3;
03079 td = vc;
03080 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03081 {
03082 TMatrixDSparseDiag md(m);
03083 for (j = 0; j < md.GetNdiags(); j++)
03084 md(j) /= 1.5;
03085 }
03086 vc = TMatrixDSparseDiag(m);
03087 ok &= VerifyVectorValue(vc,(pattern+3)/1.5,gVerbose,EPSILON);
03088 TMatrixDSparseDiag(m) *= 1.5;
03089 TMatrixDSparseDiag(m) += -3;
03090 vc = TMatrixDSparseDiag(m);
03091 ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
03092 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03093
03094 if (gVerbose)
03095 cout << "\nDone\n" << endl;
03096 StatusPrint(9,"Matrix Slices to Vectors",ok);
03097 }
03098
03099
03100
03101
03102
03103 void spstress_matrix_io()
03104 {
03105 if (gVerbose)
03106 cout << "\n---> Test matrix I/O" << endl;
03107
03108 Bool_t ok = kTRUE;
03109 const Double_t pattern = TMath::Pi();
03110
03111 TFile *f = new TFile("vmatrix.root", "RECREATE");
03112
03113 Char_t name[80];
03114 Int_t iloop = gNrLoop;
03115 while (iloop >= 0) {
03116 const Int_t msize = gSizeA[iloop];
03117
03118 const Int_t verbose = (gVerbose && iloop==gNrLoop);
03119
03120 TMatrixD m_d(msize,msize); m_d = pattern;
03121 TMatrixDSparse m(m_d);
03122 TMatrixDSparse ma;
03123 ma.Use(m);
03124
03125 if (verbose)
03126 cout << "\nWrite matrix m to database" << endl;
03127 sprintf(name,"m_%d",msize);
03128 m.Write(name);
03129
03130 if (verbose)
03131 cout << "\nWrite matrix ma which adopts to database" << endl;
03132 sprintf(name,"ma_%d",msize);
03133 ma.Write(name);
03134
03135 iloop--;
03136 }
03137
03138 if (gVerbose)
03139 cout << "\nClose database" << endl;
03140 delete f;
03141
03142 if (gVerbose)
03143 cout << "\nOpen database in read-only mode and read matrix" << endl;
03144 TFile *f1 = new TFile("vmatrix.root");
03145
03146 iloop = gNrLoop;
03147 while (iloop >= 0) {
03148 const Int_t msize = gSizeA[iloop];
03149
03150 const Int_t verbose = (gVerbose && iloop==gNrLoop);
03151
03152 TMatrixD m_d(msize,msize); m_d = pattern;
03153 TMatrixDSparse m(m_d);
03154 sprintf(name,"m_%d",msize);
03155 TMatrixDSparse *mr = (TMatrixDSparse*) f1->Get(name);
03156 sprintf(name,"ma_%d",msize);
03157 TMatrixDSparse *mar = (TMatrixDSparse*) f1->Get(name);
03158 sprintf(name,"ms_%d",msize);
03159
03160 if (verbose)
03161 cout << "\nRead matrix should be same as original" << endl;
03162 ok &= ((*mr) == m) ? kTRUE : kFALSE;
03163 ok &= ((*mar) == m) ? kTRUE : kFALSE;
03164
03165 iloop--;
03166 }
03167
03168 delete f1;
03169
03170 if (gVerbose)
03171 cout << "\nDone\n" << endl;
03172
03173 StatusPrint(10,"Matrix Persistence",ok);
03174 }
03175
03176
03177
03178
03179 void vstress_allocation(Int_t msize)
03180 {
03181 if (gVerbose)
03182 cout << "\n\n---> Test allocation and compatibility check" << endl;
03183
03184 Int_t i;
03185 Bool_t ok = kTRUE;
03186 TVectorD v1(msize);
03187 TVectorD v2(0,msize-1);
03188 TVectorD v3(1,msize);
03189 TVectorD v4(v1);
03190
03191 if (gVerbose) {
03192 cout << "\nStatus information reported for vector v3:" << endl;
03193 cout << " Lower bound ... " << v3.GetLwb() << endl;
03194 cout << " Upper bound ... " << v3.GetUpb() << endl;
03195 cout << " No. of elements " << v3.GetNoElements() << endl;
03196 }
03197
03198 if (gVerbose)
03199 cout << "\nCheck vectors 1 & 2 for compatibility" << endl;
03200 ok &= AreCompatible(v1,v2,gVerbose);
03201
03202 if (gVerbose)
03203 cout << "Check vectors 1 & 4 for compatibility" << endl;
03204 ok &= AreCompatible(v1,v4,gVerbose);
03205
03206 if (gVerbose)
03207 cout << "v2 has to be compatible with v3 after resizing to v3" << endl;
03208 v2.ResizeTo(v3);
03209 ok &= AreCompatible(v2,v3,gVerbose);
03210
03211 TVectorD v5(v1.GetUpb()+5);
03212 if (gVerbose)
03213 cout << "v1 has to be compatible with v5 after resizing to v5.upb" << endl;
03214 v1.ResizeTo(v5.GetNoElements());
03215 ok &= AreCompatible(v1,v5,gVerbose);
03216
03217 {
03218 if (gVerbose)
03219 cout << "Check that shrinking does not change remaining elements" << endl;
03220 TVectorD vb(-1,msize);
03221 for (i = vb.GetLwb(); i <= vb.GetUpb(); i++)
03222 vb(i) = i+TMath::Pi();
03223 TVectorD v = vb;
03224 ok &= ( v == vb ) ? kTRUE : kFALSE;
03225 ok &= ( v != 0 ) ? kTRUE : kFALSE;
03226 v.ResizeTo(0,msize/2);
03227 for (i = v.GetLwb(); i <= v.GetUpb(); i++)
03228 ok &= ( v(i) == vb(i) ) ? kTRUE : kFALSE;
03229 if (gVerbose)
03230 cout << "Check that expansion expands by zeros" << endl;
03231 const Int_t old_nelems = v.GetNoElements();
03232 const Int_t old_lwb = v.GetLwb();
03233 v.ResizeTo(vb);
03234 ok &= ( !(v == vb) ) ? kTRUE : kFALSE;
03235 for (i = old_lwb; i < old_lwb+old_nelems; i++)
03236 ok &= ( v(i) == vb(i) ) ? kTRUE : kFALSE;
03237 for (i = v.GetLwb(); i < old_lwb; i++)
03238 ok &= ( v(i) == 0 ) ? kTRUE : kFALSE;
03239 for (i = old_lwb+old_nelems; i <= v.GetUpb(); i++)
03240 ok &= ( v(i) == 0 ) ? kTRUE : kFALSE;
03241 }
03242
03243 if (gVerbose)
03244 cout << "\nDone\n" << endl;
03245 StatusPrint(1,"Allocation, Filling, Resizing",ok);
03246 }
03247
03248
03249
03250
03251
03252 class SinAction : public TElementActionD {
03253 void Operation(Double_t &element) const { element = TMath::Sin(element); }
03254 public:
03255 SinAction() { }
03256 };
03257
03258 class CosAction : public TElementPosActionD {
03259 Double_t factor;
03260 void Operation(Double_t &element) const { element = TMath::Cos(factor*fI); }
03261 public:
03262 CosAction(Int_t no_elems): factor(2*TMath::Pi()/no_elems) { }
03263 };
03264
03265 void vstress_element_op(Int_t vsize)
03266 {
03267 if (gVerbose)
03268 cout << "\n---> Test operations that treat each element uniformly" << endl;
03269
03270 Bool_t ok = kTRUE;
03271 const Double_t pattern = TMath::Pi();
03272
03273 TVectorD v(-1,vsize-2);
03274 TVectorD v1(v);
03275
03276 if (gVerbose)
03277 cout << "\nWriting zeros to v..." << endl;
03278 for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
03279 v(i) = 0;
03280 ok &= VerifyVectorValue(v,0.0,gVerbose,EPSILON);
03281
03282 if (gVerbose)
03283 cout << "Clearing v1 ..." << endl;
03284 v1.Zero();
03285 ok &= VerifyVectorValue(v1,0.0,gVerbose,EPSILON);
03286
03287 if (gVerbose)
03288 cout << "Comparing v1 with 0 ..." << endl;
03289 ok &= (v1 == 0) ? kTRUE : kFALSE;
03290
03291 if (gVerbose)
03292 cout << "Writing a pattern " << pattern << " by assigning to v(i)..." << endl;
03293 {
03294 for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
03295 v(i) = pattern;
03296 ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
03297 }
03298
03299 if (gVerbose)
03300 cout << "Writing the pattern by assigning to v1 as a whole ..." << endl;
03301 v1 = pattern;
03302 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03303
03304 if (gVerbose)
03305 cout << "Comparing v and v1 ..." << endl;
03306 ok &= (v == v1) ? kTRUE : kFALSE;
03307 if (gVerbose)
03308 cout << "Comparing (v=0) and v1 ..." << endl;
03309 ok &= (!(v.Zero() == v1)) ? kTRUE : kFALSE;
03310
03311 if (gVerbose)
03312 cout << "\nClear v and add the pattern" << endl;
03313 v.Zero();
03314 v += pattern;
03315 ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
03316 if (gVerbose)
03317 cout << " add the doubled pattern with the negative sign" << endl;
03318 v += -2*pattern;
03319 ok &= VerifyVectorValue(v,-pattern,gVerbose,EPSILON);
03320 if (gVerbose)
03321 cout << " subtract the trippled pattern with the negative sign" << endl;
03322 v -= -3*pattern;
03323 ok &= VerifyVectorValue(v,2*pattern,gVerbose,EPSILON);
03324
03325 if (gVerbose)
03326 cout << "\nVerify comparison operations" << endl;
03327 v = pattern;
03328 ok &= ( v == pattern && !(v != pattern) && v >= pattern && v <= pattern ) ? kTRUE : kFALSE;
03329 ok &= ( v > 0 && v >= 0 ) ? kTRUE : kFALSE;
03330 ok &= ( v > -pattern && v >= -pattern ) ? kTRUE : kFALSE;
03331 ok &= ( v < pattern+1 && v <= pattern+1 ) ? kTRUE : kFALSE;
03332 v(v.GetUpb()) += 1;
03333 ok &= ( !(v==pattern) && !(v != pattern) && v != pattern-1 ) ? kTRUE : kFALSE;
03334 ok &= ( v >= pattern && !(v > pattern) && !(v >= pattern+1) ) ? kTRUE : kFALSE;
03335 ok &= ( v <= pattern+1.001 && !(v < pattern+1) && !(v <= pattern) ) ? kTRUE : kFALSE;
03336
03337 if (gVerbose)
03338 cout << "\nAssign 2*pattern to v by repeating additions" << endl;
03339 v = 0; v += pattern; v += pattern;
03340 if (gVerbose)
03341 cout << "Assign 2*pattern to v1 by multiplying by two" << endl;
03342 v1 = pattern; v1 *= 2;
03343 ok &= VerifyVectorValue(v1,2*pattern,gVerbose,EPSILON);
03344 ok &= ( v == v1 ) ? kTRUE : kFALSE;
03345 if (gVerbose)
03346 cout << "Multiply v1 by one half returning it to the 1*pattern" << endl;
03347 v1 *= 1/2.;
03348 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03349
03350 if (gVerbose)
03351 cout << "\nAssign -pattern to v and v1" << endl;
03352 v.Zero(); v -= pattern; v1 = -pattern;
03353 ok &= VerifyVectorValue(v,-pattern,gVerbose,EPSILON);
03354 ok &= ( v == v1 ) ? kTRUE : kFALSE;
03355 if (gVerbose)
03356 cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same" << endl;
03357 v.Sqr();
03358 ok &= VerifyVectorValue(v,pattern*pattern,gVerbose,EPSILON);
03359 v.Sqrt();
03360 ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
03361 v1.Abs();
03362 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03363 ok &= ( v == v1 ) ? kTRUE : kFALSE;
03364
03365 {
03366 if (gVerbose)
03367 cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1" << endl;
03368 for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
03369 v(i) = 2*TMath::Pi()/v.GetNoElements()*i;
03370 #ifndef __CINT__
03371 SinAction s;
03372 v.Apply(s);
03373 CosAction c(v.GetNoElements());
03374 v1.Apply(c);
03375 #else
03376 for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++) {
03377 v(i) = TMath::Sin(v(i));
03378 v1(i) = TMath::Cos(2*TMath::Pi()/v1.GetNrows()*i);
03379 }
03380 #endif
03381 TVectorD v2 = v;
03382 TVectorD v3 = v1;
03383 v.Sqr();
03384 v1.Sqr();
03385 v += v1;
03386 ok &= VerifyVectorValue(v,1.,gVerbose,EPSILON);
03387 }
03388
03389 if (gVerbose)
03390 cout << "\nVerify constructor with initialization" << endl;
03391 #ifndef __CINT__
03392 TVectorD vi(0,4,0.0,1.0,2.0,3.0,4.0,"END");
03393 #else
03394 Double_t vval[] = {0.0,1.0,2.0,3.0,4.0};
03395 TVectorD vi(5,vval);
03396 #endif
03397 TVectorD vit(5);
03398 {
03399 for (Int_t i = vit.GetLwb(); i <= vit.GetUpb(); i++)
03400 vit(i) = Double_t(i);
03401 ok &= VerifyVectorIdentity(vi,vit,gVerbose,EPSILON);
03402 }
03403
03404 if (gVerbose)
03405 cout << "\nDone\n" << endl;
03406 StatusPrint(2,"Uniform vector operations",ok);
03407 }
03408
03409
03410
03411
03412
03413 void vstress_binary_op(Int_t vsize)
03414 {
03415 if (gVerbose)
03416 cout << "\n---> Test Binary Vector operations" << endl;
03417
03418 Bool_t ok = kTRUE;
03419 const Double_t pattern = TMath::Pi();
03420
03421 const Double_t epsilon = EPSILON*vsize/10;
03422
03423 TVectorD v(2,vsize+1);
03424 TVectorD v1(v);
03425
03426 if (gVerbose)
03427 cout << "\nVerify assignment of a vector to the vector" << endl;
03428 v = pattern;
03429 v1.Zero();
03430 v1 = v;
03431 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03432 ok &= ( v1 == v ) ? kTRUE : kFALSE;
03433
03434 if (gVerbose)
03435 cout << "\nAdding one vector to itself, uniform pattern " << pattern << endl;
03436 v.Zero(); v = pattern;
03437 v1 = v; v1 += v1;
03438 ok &= VerifyVectorValue(v1,2*pattern,gVerbose,EPSILON);
03439 if (gVerbose)
03440 cout << " subtracting two vectors ..." << endl;
03441 v1 -= v;
03442 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03443 if (gVerbose)
03444 cout << " subtracting the vector from itself" << endl;
03445 v1 -= v1;
03446 ok &= VerifyVectorValue(v1,0.,gVerbose,EPSILON);
03447 if (gVerbose)
03448 cout << " adding two vectors together" << endl;
03449 v1 += v;
03450 ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
03451
03452 TVectorD vp(2,vsize+1);
03453 {
03454 for (Int_t i = vp.GetLwb(); i <= vp.GetUpb(); i++)
03455 vp(i) = (i-vp.GetNoElements()/2.)*pattern;
03456 }
03457
03458 if (gVerbose) {
03459 cout << "\nArithmetic operations on vectors with not the same elements" << endl;
03460 cout << " adding vp to the zero vector..." << endl;
03461 }
03462 v.Zero();
03463 ok &= ( v == 0.0 ) ? kTRUE : kFALSE;
03464 v += vp;
03465 ok &= VerifyVectorIdentity(v,vp,gVerbose,epsilon);
03466 v1 = v;
03467 if (gVerbose)
03468 cout << " making v = 3*vp and v1 = 3*vp, via add() and succesive mult" << endl;
03469 Add(v,2.,vp);
03470 v1 += v1; v1 += vp;
03471 ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
03472 if (gVerbose)
03473 cout << " clear both v and v1, by subtracting from itself and via add()" << endl;
03474 v1 -= v1;
03475 Add(v,-3.,vp);
03476 ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
03477
03478 if (gVerbose) {
03479 cout << "\nTesting element-by-element multiplications and divisions" << endl;
03480 cout << " squaring each element with sqr() and via multiplication" << endl;
03481 }
03482 v = vp; v1 = vp;
03483 v.Sqr();
03484 ElementMult(v1,v1);
03485 ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
03486 if (gVerbose)
03487 cout << " compare (v = pattern^2)/pattern with pattern" << endl;
03488 v = pattern; v1 = pattern;
03489 v.Sqr();
03490 ElementDiv(v,v1);
03491 ok &= VerifyVectorValue(v,pattern,gVerbose,epsilon);
03492 if (gVerbose)
03493 Compare(v1,v);
03494
03495 if (gVerbose)
03496 cout << "\nDone\n" << endl;
03497 StatusPrint(3,"Binary vector element-by-element operations",ok);
03498 }
03499
03500
03501
03502
03503
03504 void vstress_norms(Int_t vsize)
03505 {
03506 if (gVerbose)
03507 cout << "\n---> Verify norm calculations" << endl;
03508
03509 Bool_t ok = kTRUE;
03510 const Double_t pattern = 10.25;
03511
03512 if ( vsize % 2 == 1 )
03513 Fatal("vstress_norms", "size of the vector to test must be even for this test\n");
03514
03515 TVectorD v(vsize);
03516 TVectorD v1(v);
03517
03518 if (gVerbose)
03519 cout << "\nAssign " << pattern << " to all the elements and check norms" << endl;
03520 v = pattern;
03521 if (gVerbose)
03522 cout << " 1. norm should be pattern*no_elems" << endl;
03523 ok &= ( v.Norm1() == pattern*v.GetNoElements() ) ? kTRUE : kFALSE;
03524 if (gVerbose)
03525 cout << " Square of the 2. norm has got to be pattern^2 * no_elems" << endl;
03526 ok &= ( v.Norm2Sqr() == (pattern*pattern)*v.GetNoElements() ) ? kTRUE : kFALSE;
03527 if (gVerbose)
03528 cout << " Inf norm should be pattern itself" << endl;
03529 ok &= ( v.NormInf() == pattern ) ? kTRUE : kFALSE;
03530 if (gVerbose)
03531 cout << " Scalar product of vector by itself is the sqr(2. vector norm)" << endl;
03532 ok &= ( v.Norm2Sqr() == v*v ) ? kTRUE : kFALSE;
03533
03534 Double_t ap_step = 1;
03535 Double_t ap_a0 = -pattern;
03536 Int_t n = v.GetNoElements();
03537 if (gVerbose) {
03538 cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
03539 "\nand the difference " << ap_step << endl;
03540 }
03541 {
03542 for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
03543 v(i) = (i-v.GetLwb())*ap_step + ap_a0;
03544 }
03545 Int_t l = TMath::Min(TMath::Max((int)TMath::Ceil(-ap_a0/ap_step),0),n);
03546 Double_t norm = (2*ap_a0+(l+n-1)*ap_step)/2*(n-l) +
03547 (-2*ap_a0-(l-1)*ap_step)/2*l;
03548 if (gVerbose)
03549 cout << " 1. norm should be " << norm << endl;
03550 ok &= ( v.Norm1() == norm ) ? kTRUE : kFALSE;
03551 norm = n*( (ap_a0*ap_a0)+ap_a0*ap_step*(n-1)+(ap_step*ap_step)*(n-1)*(2*n-1)/6);
03552 if (gVerbose) {
03553 cout << " Square of the 2. norm has got to be "
03554 "n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << endl;
03555 }
03556 ok &= ( TMath::Abs( (v.Norm2Sqr()-norm)/norm ) < EPSILON ) ? kTRUE : kFALSE;
03557
03558 norm = TMath::Max(TMath::Abs(v(v.GetLwb())),TMath::Abs(v(v.GetUpb())));
03559 if (gVerbose)
03560 cout << " Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm << endl;
03561 ok &= ( v.NormInf() == norm ) ? kTRUE : kFALSE;
03562 if (gVerbose)
03563 cout << " Scalar product of vector by itself is the sqr(2. vector norm)" << endl;
03564 ok &= ( v.Norm2Sqr() == v*v ) ? kTRUE : kFALSE;
03565
03566 #if 0
03567 v1.Zero();
03568 Compare(v,v1);
03569 #endif
03570
03571 if (gVerbose)
03572 cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)..." << endl;
03573 {
03574 for (Int_t i = 0; i < v1.GetNoElements(); i++)
03575 v1(i+v1.GetLwb()) = v(v.GetUpb()-i) * ( i % 2 == 1 ? -1 : 1 );
03576 }
03577 if (gVerbose)
03578 cout << "||v1|| has got to be equal ||v|| regardless of the norm def" << endl;
03579 ok &= ( v1.Norm1() == v.Norm1() ) ? kTRUE : kFALSE;
03580 ok &= ( v1.Norm2Sqr() == v.Norm2Sqr() ) ? kTRUE : kFALSE;
03581 ok &= ( v1.NormInf() == v.NormInf() ) ? kTRUE : kFALSE;
03582 if (gVerbose)
03583 cout << "But the scalar product has to be zero" << endl;
03584 ok &= ( v1 * v == 0 ) ? kTRUE : kFALSE;
03585
03586 if (gVerbose)
03587 cout << "\nDone\n" << endl;
03588 StatusPrint(4,"Vector Norms",ok);
03589 }
03590
03591
03592
03593
03594
03595 void vstress_matrix_slices(Int_t vsize)
03596 {
03597 if (gVerbose)
03598 cout << "\n---> Test operations with vectors and matrix slices" << endl;
03599
03600 Bool_t ok = kTRUE;
03601 const Double_t pattern = 8.625;
03602
03603 TVectorD vc(0,vsize);
03604 TVectorD vr(0,vsize+1);
03605 TMatrixD m(0,vsize,0,vsize+1);
03606
03607 Int_t i,j;
03608 if (gVerbose)
03609 cout << "\nCheck modifying the matrix column-by-column" << endl;
03610 m = pattern;
03611 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03612 for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
03613 TMatrixDColumn(m,i) = pattern-1;
03614 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03615 TMatrixDColumn(m,i) *= 2;
03616 vc = TMatrixDColumn(m,i);
03617 ok &= VerifyVectorValue(vc,2*(pattern-1),gVerbose);
03618 vc = TMatrixDColumn(m,i+1 > m.GetColUpb() ? m.GetColLwb() : i+1);
03619 ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
03620 TMatrixDColumn(m,i) *= 0.5;
03621 TMatrixDColumn(m,i) += 1;
03622 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03623 }
03624
03625 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03626 for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
03627 vc = pattern+1;
03628 TMatrixDColumn(m,i) = vc;
03629 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03630 {
03631 TMatrixDColumn mc(m,i);
03632 for (j = m.GetRowLwb(); j <= m.GetRowUpb(); j++)
03633 mc(j) *= 4;
03634 }
03635 vc = TMatrixDColumn(m,i);
03636 ok &= VerifyVectorValue(vc,4*(pattern+1),gVerbose,EPSILON);
03637 TMatrixDColumn(m,i) *= 0.25;
03638 TMatrixDColumn(m,i) += -1;
03639 vc = TMatrixDColumn(m,i);
03640 ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
03641 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03642 }
03643
03644 if (gVerbose)
03645 cout << "\nCheck modifying the matrix row-by-row" << endl;
03646 m = pattern;
03647 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03648 for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
03649 TMatrixDRow(m,i) = pattern+2;
03650 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03651 vr = TMatrixDRow(m,i);
03652 ok &= VerifyVectorValue(vr,pattern+2,gVerbose,EPSILON);
03653 vr = TMatrixDRow(m,i+1 > m.GetRowUpb() ? m.GetRowLwb() : i+1);
03654 ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
03655 TMatrixDRow(m,i) += -2;
03656 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03657 }
03658
03659 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03660 for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
03661 vr = pattern-2;
03662 TMatrixDRow(m,i) = vr;
03663 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03664 {
03665 TMatrixDRow mr(m,i);
03666 for (j = m.GetColLwb(); j <= m.GetColUpb(); j++)
03667 mr(j) *= 8;
03668 }
03669 vr = TMatrixDRow(m,i);
03670 ok &= VerifyVectorValue(vr,8*(pattern-2),gVerbose,EPSILON);
03671 TMatrixDRow(m,i) *= 1./8;
03672 TMatrixDRow(m,i) += 2;
03673 vr = TMatrixDRow(m,i);
03674 ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
03675 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03676 }
03677
03678 if (gVerbose)
03679 cout << "\nCheck modifying the matrix diagonal" << endl;
03680 m = pattern;
03681 TMatrixDDiag td = m;
03682 td = pattern-3;
03683 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03684 vc = TMatrixDDiag(m);
03685 ok &= VerifyVectorValue(vc,pattern-3,gVerbose,EPSILON);
03686 td += 3;
03687 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03688 vc = pattern+3;
03689 td = vc;
03690 ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
03691 {
03692 TMatrixDDiag md(m);
03693 for (j = 0; j < md.GetNdiags(); j++)
03694 md(j) /= 1.5;
03695 }
03696 vc = TMatrixDDiag(m);
03697 ok &= VerifyVectorValue(vc,(pattern+3)/1.5,gVerbose,EPSILON);
03698 TMatrixDDiag(m) *= 1.5;
03699 TMatrixDDiag(m) += -3;
03700 vc = TMatrixDDiag(m);
03701 ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
03702 ok &= ( m == pattern ) ? kTRUE : kFALSE;
03703
03704 if (gVerbose) {
03705 cout << "\nCheck out to see that multiplying by diagonal is column-wise"
03706 "\nmatrix multiplication" << endl;
03707 }
03708 TMatrixD mm(m);
03709 TMatrixD m1(m.GetRowLwb(),TMath::Max(m.GetRowUpb(),m.GetColUpb()),
03710 m.GetColLwb(),TMath::Max(m.GetRowUpb(),m.GetColUpb()));
03711 TVectorD vc1(vc),vc2(vc);
03712 for (i = m.GetRowLwb(); i < m.GetRowUpb(); i++)
03713 TMatrixDRow(m,i) = pattern+i;
03714 mm = m;
03715
03716 m1 = pattern+10;
03717 for (i = vr.GetLwb(); i <= vr.GetUpb(); i++)
03718 vr(i) = i+2;
03719 TMatrixDDiag td2 = m1;
03720 td2 = vr;
03721 ok &= ( !(m1 == pattern+10) ) ? kTRUE : kFALSE;
03722
03723 m *= TMatrixDDiag(m1);
03724 for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
03725 vc1 = TMatrixDColumn(mm,i);
03726 vc1 *= vr(i);
03727 vc2 = TMatrixDColumn(m,i);
03728 ok &= VerifyVectorIdentity(vc1,vc2,gVerbose,EPSILON);
03729 }
03730
03731 if (gVerbose)
03732 cout << "\nDone\n" << endl;
03733 StatusPrint(5,"Matrix Slices to Vectors",ok);
03734 }
03735
03736
03737
03738
03739
03740 void vstress_vector_io()
03741 {
03742 if (gVerbose)
03743 cout << "\n---> Test vector I/O" << endl;
03744
03745 Bool_t ok = kTRUE;
03746 const Double_t pattern = TMath::Pi();
03747
03748 TFile *f = new TFile("vvector.root","RECREATE");
03749
03750 Char_t name[80];
03751 Int_t iloop = gNrLoop;
03752 while (iloop >= 0) {
03753 const Int_t msize = gSizeA[iloop];
03754
03755 const Int_t verbose = (gVerbose && iloop==gNrLoop);
03756
03757 TVectorD v(msize);
03758 v = pattern;
03759
03760 Double_t *pattern_array = new Double_t[msize];
03761 for (Int_t i = 0; i < msize; i++)
03762 pattern_array[i] = pattern;
03763 TVectorD va;
03764 va.Use(msize,pattern_array);
03765
03766 if (verbose)
03767 cout << "\nWrite vector v to database" << endl;
03768
03769 sprintf(name,"v_%d",msize);
03770 v.Write(name);
03771 sprintf(name,"va_%d",msize);
03772 va.Write(name);
03773
03774 delete [] pattern_array;
03775
03776 iloop--;
03777 }
03778
03779 if (gVerbose)
03780 cout << "\nClose database" << endl;
03781 delete f;
03782
03783 if (gVerbose)
03784 cout << "\nOpen database in read-only mode and read vector" << endl;
03785 TFile *f1 = new TFile("vvector.root");
03786
03787 iloop = gNrLoop;
03788 while (iloop >= 0) {
03789 const Int_t msize = gSizeA[iloop];
03790
03791 const Int_t verbose = (gVerbose && iloop==gNrLoop);
03792
03793 TVectorD v(msize);
03794 v = pattern;
03795
03796 sprintf(name,"v_%d",msize);
03797 TVectorD *vr = (TVectorD*) f1->Get(name);
03798 sprintf(name,"va_%d",msize);
03799 TVectorD *var = (TVectorD*) f1->Get(name);
03800
03801 if (verbose)
03802 cout << "\nRead vector should be same as original still in memory" << endl;
03803 ok &= ((*vr) == v) ? kTRUE : kFALSE;
03804 ok &= ((*var) == v) ? kTRUE : kFALSE;
03805
03806 iloop--;
03807 }
03808
03809 delete f1;
03810
03811 if (gVerbose)
03812 cout << "\nDone\n" << endl;
03813 StatusPrint(6,"Vector Persistence",ok);
03814 }
03815
03816 Bool_t test_svd_expansion(const TMatrixD &A)
03817 {
03818 if (gVerbose)
03819 cout << "\n\nSVD-decompose matrix A and check if we can compose it back\n" <<endl;
03820
03821 Bool_t ok = kTRUE;
03822
03823 TDecompSVD svd(A);
03824 if (gVerbose) {
03825 cout << "left factor U" <<endl;
03826 svd.GetU().Print();
03827 cout << "Vector of Singular values" <<endl;
03828 svd.GetSig().Print();
03829 cout << "right factor V" <<endl;
03830 svd.GetV().Print();
03831 }
03832
03833 {
03834 if (gVerbose)
03835 cout << "\tchecking that U is orthogonal indeed, i.e., U'U=E and UU'=E" <<endl;
03836 const Int_t nRows = svd.GetU().GetNrows();
03837 const Int_t nCols = svd.GetU().GetNcols();
03838 TMatrixD E1(nRows,nRows); E1.UnitMatrix();
03839 TMatrixD E2(nCols,nCols); E2.UnitMatrix();
03840 TMatrixD ut(TMatrixD::kTransposed,svd.GetU());
03841 ok &= VerifyMatrixIdentity(ut * svd.GetU(),E2,gVerbose,100*EPSILON);
03842 ok &= VerifyMatrixIdentity(svd.GetU() * ut,E1,gVerbose,100*EPSILON);
03843 }
03844
03845 {
03846 if (gVerbose)
03847 cout << "\tchecking that V is orthogonal indeed, i.e., V'V=E and VV'=E" <<endl;
03848 const Int_t nRows = svd.GetV().GetNrows();
03849 const Int_t nCols = svd.GetV().GetNcols();
03850 TMatrixD E1(nRows,nRows); E1.UnitMatrix();
03851 TMatrixD E2(nCols,nCols); E2.UnitMatrix();
03852 TMatrixD vt(TMatrixD::kTransposed,svd.GetV());
03853 ok &= VerifyMatrixIdentity(vt * svd.GetV(),E2,gVerbose,100*EPSILON);
03854 ok &= VerifyMatrixIdentity(svd.GetV() * vt,E1,gVerbose,100*EPSILON);
03855 }
03856
03857 {
03858 if (gVerbose)
03859 cout << "\tchecking that U*Sig*V' is indeed A" <<endl;
03860 const Int_t nRows = svd.GetU().GetNrows();
03861 const Int_t nCols = svd.GetV().GetNcols();
03862 TMatrixD s(nRows,nCols);
03863 TMatrixDDiag diag(s); diag = svd.GetSig();
03864 TMatrixD vt(TMatrixD::kTransposed,svd.GetV());
03865 TMatrixD tmp = s * vt;
03866 ok &= VerifyMatrixIdentity(A,svd.GetU() * tmp,gVerbose,100*EPSILON);
03867 if (gVerbose) {
03868 cout << "U*Sig*V'" <<endl;
03869 (svd.GetU()*tmp).Print();
03870 }
03871 }
03872
03873 return ok;
03874 }
03875
03876 #ifndef __CINT__
03877
03878 class MakeMatrix : public TMatrixDLazy {
03879 const Double_t *array;
03880 Int_t no_elems;
03881 void FillIn(TMatrixD& m) const {
03882 R__ASSERT( m.GetNrows() * m.GetNcols() == no_elems );
03883 const Double_t *ap = array;
03884 Double_t *mp = m.GetMatrixArray();
03885 for (Int_t i = 0; i < no_elems; i++)
03886 *mp++ = *ap++;
03887 }
03888
03889 public:
03890 MakeMatrix(Int_t nrows,Int_t ncols,
03891 const Double_t *_array,Int_t _no_elems)
03892 :TMatrixDLazy(nrows,ncols), array(_array), no_elems(_no_elems) {}
03893 MakeMatrix(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
03894 const Double_t *_array,Int_t _no_elems)
03895 : TMatrixDLazy(row_lwb,row_upb,col_lwb,col_upb),
03896 array(_array), no_elems(_no_elems) {}
03897 };
03898 #else
03899 TMatrixD MakeMatrix(Int_t nrows,Int_t ncols,
03900 const Double_t *_array,Int_t _no_elems)
03901 {
03902 TMatrixD m(nrows,ncols,_array);
03903 return m;
03904 }
03905 #endif
03906
03907 void astress_decomp()
03908 {
03909 Bool_t ok = kTRUE;
03910
03911 {
03912 if (gVerbose)
03913 cout << "\nBrandt example page 503\n" <<endl;
03914
03915 Double_t array0[] = { -2,1,0,0, 2,1,0,0, 0,0,0,0, 0,0,0,0 };
03916 ok &= test_svd_expansion(MakeMatrix(4,4,array0,sizeof(array0)/sizeof(array0[0])));
03917 }
03918
03919 {
03920 if (gVerbose)
03921 cout << "\nRotated by PI/2 Matrix Diag(1,4,9)\n" <<endl;
03922
03923 Double_t array1[] = {0,-4,0, 1,0,0, 0,0,9 };
03924 ok &= test_svd_expansion(MakeMatrix(3,3,array1,sizeof(array1)/sizeof(array1[0])));
03925 }
03926
03927 {
03928 if (gVerbose)
03929 cout << "\nExample from the Forsythe, Malcolm, Moler's book\n" <<endl;
03930
03931 Double_t array2[] =
03932 { 1,6,11, 2,7,12, 3,8,13, 4,9,14, 5,10,15};
03933 ok &= test_svd_expansion(MakeMatrix(5,3,array2,sizeof(array2)/sizeof(array2[0])));
03934 }
03935
03936 {
03937 if (gVerbose)
03938 cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
03939 "Singular numbers are 0, 19.5959, 20, 0, 35.3270\n" <<endl;
03940
03941 Double_t array3[] =
03942 { 22, 10, 2, 3, 7, 14, 7, 10, 0, 8,
03943 -1, 13, -1, -11, 3, -3, -2, 13, -2, 4,
03944 9, 8, 1, -2, 4, 9, 1, -7, 5, -1,
03945 2, -6, 6, 5, 1, 4, 5, 0, -2, 2 };
03946 ok &= test_svd_expansion(MakeMatrix(8,5,array3,sizeof(array3)/sizeof(array3[0])));
03947 }
03948
03949 {
03950 if (gVerbose)
03951 cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
03952 "Ordered singular numbers are Sig[21-k] = sqrt(k*(k-1))\n" <<endl;
03953 TMatrixD A(21,20);
03954 for (Int_t irow = A.GetRowLwb(); irow <= A.GetRowUpb(); irow++)
03955 for (Int_t icol = A.GetColLwb(); icol <= A.GetColUpb(); icol++)
03956 A(irow,icol) = (irow>icol ? 0 : ((irow==icol) ? 20-icol : -1));
03957
03958 ok &= test_svd_expansion(A);
03959 }
03960
03961 if (0)
03962 {
03963 if (gVerbose) {
03964 cout << "\nTest by W. Meier <wmeier@manu.com> to catch an obscure "
03965 << "bug in QR\n" <<endl;
03966 cout << "expect singular values to be\n"
03967 << "1.4666e-024 1.828427 3.828427 4.366725 7.932951\n" <<endl;
03968 }
03969
03970 Double_t array4[] =
03971 { 1, 2, 0, 0, 0,
03972 0, 2, 3, 0, 0,
03973 0, 0, 0, 4, 0,
03974 0, 0, 0, 4, 5,
03975 0, 0, 0, 0, 5 };
03976 ok &= test_svd_expansion(MakeMatrix(5,5,array4,sizeof(array4)/sizeof(array4[0])));
03977 }
03978
03979 {
03980 const TMatrixD m1 = THilbertMatrixD(5,5);
03981 TDecompLU lu(m1);
03982 ok &= VerifyMatrixIdentity(lu.GetMatrix(),m1,gVerbose,100*EPSILON);
03983 }
03984
03985 {
03986 const TMatrixD m2 = THilbertMatrixD(5,5);
03987 const TMatrixDSym mtm(TMatrixDSym::kAtA,m2);
03988 TDecompChol chol(mtm);
03989 ok &= VerifyMatrixIdentity(chol.GetMatrix(),mtm,gVerbose,100*EPSILON);
03990 }
03991
03992 if (gVerbose)
03993 cout << "\nDone" <<endl;
03994
03995 StatusPrint(1,"Decomposition / Reconstruction",ok);
03996 }
03997
03998 void astress_lineqn()
03999 {
04000 if (gVerbose)
04001 cout << "\nSolve Ax=b where A is a Hilbert matrix and b(i) = sum_j Aij\n" <<endl;
04002
04003 Bool_t ok = kTRUE;
04004
04005 Int_t iloop = gNrLoop;
04006 Int_t nr = 0;
04007
04008 while (iloop >= 0) {
04009 const Int_t msize = gSizeA[iloop];
04010
04011 const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
04012
04013 if (verbose)
04014 cout << "\nSolve Ax=b for size = " << msize <<endl;
04015
04016
04017
04018
04019 TMatrixDSym m1 = THilbertMatrixDSym(-1,msize-2);
04020 TMatrixDDiag diag1(m1);
04021 diag1 += 1.;
04022
04023 TVectorD rowsum1(-1,msize-2); rowsum1.Zero();
04024 TVectorD colsum1(-1,msize-2); colsum1.Zero();
04025 for (Int_t irow = m1.GetRowLwb(); irow <= m1.GetColUpb(); irow++) {
04026 for (Int_t icol = m1.GetColLwb(); icol <= m1.GetColUpb(); icol++) {
04027 rowsum1(irow) += m1(irow,icol);
04028 colsum1(icol) += m1(irow,icol);
04029 }
04030 }
04031
04032 TMatrixDSym m2 = THilbertMatrixDSym(msize);
04033 TMatrixDDiag diag2(m2);
04034 diag2 += 1.;
04035
04036 TVectorD rowsum2(msize); rowsum2.Zero();
04037 TVectorD colsum2(msize); colsum2.Zero();
04038 for (Int_t irow = m2.GetRowLwb(); irow <= m2.GetColUpb(); irow++) {
04039 for (Int_t icol = m2.GetColLwb(); icol <= m2.GetColUpb(); icol++) {
04040 rowsum2(irow) += m2(irow,icol);
04041 colsum2(icol) += m2(irow,icol);
04042 }
04043 }
04044
04045 TVectorD b1(-1,msize-2);
04046 TVectorD b2(msize);
04047 {
04048 TDecompLU lu(m1,1.0e-20);
04049 b1 = rowsum1;
04050 lu.Solve(b1);
04051 if (msize < 10)
04052 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04053 b1 = colsum1;
04054 lu.TransSolve(b1);
04055 if (msize < 10)
04056 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04057 }
04058
04059 {
04060 TDecompChol chol(m1,1.0e-20);
04061 b1 = rowsum1;
04062 chol.Solve(b1);
04063 if (msize < 10)
04064 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04065 b1 = colsum1;
04066 chol.TransSolve(b1);
04067 if (msize < 10)
04068 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04069 }
04070
04071 {
04072 TDecompQRH qrh1(m1,1.0e-20);
04073 b1 = rowsum1;
04074 qrh1.Solve(b1);
04075 if (msize < 10)
04076 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04077
04078 TDecompQRH qrh2(m2,1.0e-20);
04079 b2 = colsum2;
04080 qrh2.TransSolve(b2);
04081 if (msize < 10)
04082 ok &= VerifyVectorValue(b2,1.0,verbose,msize*EPSILON);
04083 }
04084
04085 {
04086 TDecompSVD svd1(m1);
04087 b1 = rowsum1;
04088 svd1.Solve(b1);
04089 if (msize < 10)
04090 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04091
04092 TDecompSVD svd2(m2);
04093 b2 = colsum2;
04094 svd2.TransSolve(b2);
04095 if (msize < 10)
04096 ok &= VerifyVectorValue(b2,1.0,verbose,msize*EPSILON);
04097 }
04098
04099 {
04100 TDecompBK bk(m1,1.0e-20);
04101 b1 = rowsum1;
04102 bk.Solve(b1);
04103 if (msize < 10)
04104 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04105 b1 = colsum1;
04106 bk.TransSolve(b1);
04107 if (msize < 10)
04108 ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
04109 }
04110
04111 #ifndef __CINT__
04112 if (nr >= Int_t(1.0e+5/msize/msize)) {
04113 #else
04114 if (nr >= Int_t(1.0e+3/msize/msize)) {
04115 #endif
04116 nr = 0;
04117 iloop--;
04118 } else
04119 nr++;
04120
04121 if (!ok) {
04122 if (gVerbose)
04123 cout << "Linear Equations failed for size " << msize << endl;
04124 break;
04125 }
04126 }
04127
04128 if (gVerbose)
04129 cout << "\nDone" <<endl;
04130
04131 StatusPrint(2,"Linear Equations",ok);
04132 }
04133
04134 void astress_pseudo()
04135 {
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145
04146
04147 Bool_t ok = kTRUE;
04148
04149
04150 enum {nrows = 4, ncols = 3};
04151 #ifndef __CINT__
04152 const Double_t A_data [nrows*ncols] =
04153 #else
04154 const Double_t A_data [12] =
04155 #endif
04156 {0, 0, 0,
04157 0, 0, 0,
04158 1, 1, 0,
04159 4, 3, 0};
04160 TMatrixD A(nrows,ncols,A_data);
04161
04162
04163
04164 TDecompSVD svd(A);
04165 if (gVerbose)
04166 cout << "\ncondition number " << svd.Condition() << endl;
04167
04168
04169 TMatrixD Ainv(nrows,nrows); Ainv.UnitMatrix();
04170 svd.MultiSolve(Ainv);
04171 Ainv.ResizeTo(ncols,nrows);
04172 TMatrixD Ainv2(nrows,nrows);
04173 svd.Invert(Ainv2);
04174 Ainv2.ResizeTo(ncols,nrows);
04175 ok &= VerifyMatrixIdentity(Ainv,Ainv2,gVerbose,EPSILON);
04176
04177 if (gVerbose) {
04178 cout << "\nChecking the Moore-Penrose conditions for the Ainv" << endl;
04179 cout << "\nVerify that Ainv * A * Ainv is Ainv" << endl;
04180 }
04181 ok &= VerifyMatrixIdentity(Ainv * A * Ainv, Ainv,gVerbose,EPSILON);
04182
04183 if (gVerbose)
04184 cout << "\nVerify that A * Ainv * A is A" << endl;
04185 ok &= VerifyMatrixIdentity(A * Ainv * A, A,gVerbose,EPSILON);
04186
04187 if (gVerbose)
04188 cout << "\nVerify that Ainv * A is symmetric" << endl;
04189 ok &= VerifyMatrixIdentity(Ainv * A, (Ainv*A).T(),gVerbose,EPSILON);
04190
04191 if (gVerbose)
04192 cout << "\nVerify that A * Ainv is symmetric" << endl;
04193 ok &= VerifyMatrixIdentity(A * Ainv, (A*Ainv).T(),gVerbose,EPSILON);
04194
04195 StatusPrint(3,"Pseudo-Inverse, Moore-Penrose",ok);
04196 }
04197
04198 void astress_eigen(Int_t msize)
04199 {
04200 Bool_t ok = kTRUE;
04201
04202
04203
04204
04205
04206
04207
04208 const TMatrixD m = THilbertMatrixD(msize,msize);
04209
04210 TDecompSVD svd(m);
04211 TVectorD sig = svd.GetSig(); sig.Sqr();
04212
04213
04214 TMatrixDSym mtm1(TMatrixDSym::kAtA,m);
04215 const TMatrixDSymEigen eigen1(mtm1);
04216 const TVectorD eigenVal1 = eigen1.GetEigenValues();
04217
04218 ok &= VerifyVectorIdentity(sig,eigenVal1,gVerbose,EPSILON);
04219
04220
04221 const TMatrixD mtm1x = mtm1 * eigen1.GetEigenVectors();
04222 TMatrixD lam_x1 = eigen1.GetEigenVectors();
04223 lam_x1.NormByRow(eigen1.GetEigenValues(),"");
04224
04225 ok &= VerifyMatrixIdentity(mtm1x,lam_x1,gVerbose,EPSILON);
04226
04227
04228 TMatrixD mtm2(m,TMatrixD::kTransposeMult,m);
04229 const TMatrixDEigen eigen2(mtm2);
04230 const TVectorD eigenVal2 = eigen2.GetEigenValuesRe();
04231
04232 ok &= VerifyVectorIdentity(sig,eigenVal2,gVerbose,EPSILON);
04233
04234
04235 const TMatrixD mtm2x = mtm2 * eigen2.GetEigenVectors();
04236 TMatrixD lam_x2 = eigen2.GetEigenVectors();
04237 lam_x2.NormByRow(eigen2.GetEigenValuesRe(),"");
04238
04239 ok &= VerifyMatrixIdentity(mtm2x,lam_x2,gVerbose,EPSILON);
04240
04241
04242 const TVectorD eigenValIm = eigen2.GetEigenValuesIm();
04243 TVectorD epsVec(msize); epsVec = EPSILON;
04244 ok &= VerifyVectorIdentity(epsVec,eigenValIm,gVerbose,EPSILON);
04245
04246 StatusPrint(4,"Eigen - Values/Vectors",ok);
04247 }
04248
04249 void astress_decomp_io(Int_t msize)
04250 {
04251
04252
04253
04254
04255 if (gVerbose)
04256 cout << "\n---> Test decomp I/O" << endl;
04257
04258 Bool_t ok = kTRUE;
04259
04260 const TMatrixDSym m = THilbertMatrixDSym(msize);
04261 TVectorD rowsum(msize); rowsum.Zero();
04262 TVectorD colsum(msize); colsum.Zero();
04263
04264 for (Int_t irow = 0; irow < m.GetNrows(); irow++) {
04265 for (Int_t icol = 0; icol < m.GetNcols(); icol++) {
04266 rowsum(irow) += m(irow,icol);
04267 colsum(icol) += m(irow,icol);
04268 }
04269 }
04270
04271 if (gVerbose)
04272 cout << "\nWrite decomp m to database" << endl;
04273
04274 TFile *f = new TFile("vdecomp.root", "RECREATE");
04275
04276 TDecompLU lu(m,1.0e-20);
04277 TDecompQRH qrh(m,1.0e-20);
04278 TDecompChol chol(m,1.0e-20);
04279 TDecompSVD svd(m);
04280 TDecompBK bk(m,1.0e-20);
04281
04282 lu.Write("lu");
04283 qrh.Write("qrh");
04284 chol.Write("chol");
04285 svd.Write("svd");
04286 bk.Write("bk");
04287
04288 if (gVerbose)
04289 cout << "\nClose database" << endl;
04290 delete f;
04291
04292 if (gVerbose)
04293 cout << "\nOpen database in read-only mode and read matrix" << endl;
04294 TFile *f1 = new TFile("vdecomp.root");
04295
04296 if (gVerbose)
04297 cout << "\nRead decompositions should create same solutions" << endl;
04298 {
04299 TDecompLU *rlu = (TDecompLU*) f1->Get("lu");
04300 TVectorD b1(rowsum);
04301 lu.Solve(b1);
04302 TVectorD b2(rowsum);
04303 rlu->Solve(b2);
04304 ok &= (b1 == b2);
04305 b1 = colsum;
04306 lu.TransSolve(b1);
04307 b2 = colsum;
04308 rlu->TransSolve(b2);
04309 ok &= (b1 == b2);
04310 }
04311
04312 {
04313 TDecompChol *rchol = (TDecompChol*) f1->Get("chol");
04314 TVectorD b1(rowsum);
04315 chol.Solve(b1);
04316 TVectorD b2(rowsum);
04317 rchol->Solve(b2);
04318 ok &= (b1 == b2);
04319 b1 = colsum;
04320 chol.TransSolve(b1);
04321 b2 = colsum;
04322 rchol->TransSolve(b2);
04323 ok &= (b1 == b2);
04324 }
04325
04326 {
04327 TDecompQRH *rqrh = (TDecompQRH*) f1->Get("qrh");
04328 TVectorD b1(rowsum);
04329 qrh.Solve(b1);
04330 TVectorD b2(rowsum);
04331 rqrh->Solve(b2);
04332 ok &= (b1 == b2);
04333 b1 = colsum;
04334 qrh.TransSolve(b1);
04335 b2 = colsum;
04336 rqrh->TransSolve(b2);
04337 ok &= (b1 == b2);
04338 }
04339
04340 {
04341 TDecompSVD *rsvd = (TDecompSVD*) f1->Get("svd");
04342 TVectorD b1(rowsum);
04343 svd.Solve(b1);
04344 TVectorD b2(rowsum);
04345 rsvd->Solve(b2);
04346 ok &= (b1 == b2);
04347 b1 = colsum;
04348 svd.TransSolve(b1);
04349 b2 = colsum;
04350 rsvd->TransSolve(b2);
04351 ok &= (b1 == b2);
04352 }
04353
04354 {
04355 TDecompBK *rbk = (TDecompBK*) f1->Get("bk");
04356 TVectorD b1(rowsum);
04357 bk.Solve(b1);
04358 TVectorD b2(rowsum);
04359 rbk->Solve(b2);
04360 ok &= (b1 == b2);
04361 b1 = colsum;
04362 bk.TransSolve(b1);
04363 b2 = colsum;
04364 rbk->TransSolve(b2);
04365 ok &= (b1 == b2);
04366 }
04367
04368 delete f1;
04369
04370 if (gVerbose)
04371 cout << "\nDone\n" << endl;
04372
04373 StatusPrint(5,"Decomposition Persistence",ok);
04374 }
04375
04376 void stress_backward_io()
04377 {
04378 TFile::SetCacheFileDir(".");
04379 TFile *f = TFile::Open("http://root.cern.ch/files/linearIO.root","CACHEREAD");
04380
04381 TMatrixF mf1 = THilbertMatrixF(-5,5,-5,5);
04382 mf1[1][2] = TMath::Pi();
04383 TMatrixFSym mf2 = THilbertMatrixFSym(-5,5);
04384 TVectorF vf_row(mf1.GetRowLwb(),mf1.GetRowUpb()); vf_row = TMatrixFRow(mf1,3);
04385
04386 TMatrixF *mf1_r = (TMatrixF*) f->Get("mf1");
04387 TMatrixFSym *mf2_r = (TMatrixFSym*) f->Get("mf2");
04388 TVectorF *vf_row_r = (TVectorF*) f->Get("vf_row");
04389
04390 Bool_t ok = kTRUE;
04391
04392 ok &= ((*mf1_r) == mf1) ? kTRUE : kFALSE;
04393 ok &= ((*mf2_r) == mf2) ? kTRUE : kFALSE;
04394 ok &= ((*vf_row_r) == vf_row) ? kTRUE : kFALSE;
04395
04396 TMatrixD md1 = THilbertMatrixD(-5,5,-5,5);
04397 md1[1][2] = TMath::Pi();
04398 TMatrixDSym md2 = THilbertMatrixDSym(-5,5);
04399 TMatrixDSparse md3 = md1;
04400 TVectorD vd_row(md1.GetRowLwb(),md1.GetRowUpb()); vd_row = TMatrixDRow(md1,3);
04401
04402 TMatrixD *md1_r = (TMatrixD*) f->Get("md1");
04403 TMatrixDSym *md2_r = (TMatrixDSym*) f->Get("md2");
04404 TMatrixDSparse *md3_r = (TMatrixDSparse*) f->Get("md3");
04405 TVectorD *vd_row_r = (TVectorD*) f->Get("vd_row");
04406
04407 ok &= ((*md1_r) == md1) ? kTRUE : kFALSE;
04408 ok &= ((*md2_r) == md2) ? kTRUE : kFALSE;
04409 ok &= ((*md3_r) == md3) ? kTRUE : kFALSE;
04410 ok &= ((*vd_row_r) == vd_row) ? kTRUE : kFALSE;
04411
04412 StatusPrint(1,"Streamers",ok);
04413 }
04414
04415 void cleanup()
04416 {
04417 gSystem->Unlink("vmatrix.root");
04418 gSystem->Unlink("vvector.root");
04419 gSystem->Unlink("vdecomp.root");
04420 }
