Revision 736681ec963848c66b75f8250ee2f1bc70027cfa authored by Philippe Canal on 18 February 2016, 22:45:16 UTC, committed by Philippe Canal on 18 February 2016, 22:45:16 UTC
1 parent ce57c08
stressLinear.cxx
// @(#)root/test:$name: $:$id: stressLinear.cxx,v 1.15 2002/10/25 10:47:51 rdm exp $
// Authors: Fons Rademakers, Eddy Offermann Jan 2004
/////////////////////////////////////////////////////////////////
//
// R O O T LINEAR ALGEBRA T E S T S U I T E and B E N C H M A R K S
// =====================================================================
//
// The suite of programs below tests many elements of the vector, matrix
// and matrix decomposition classes.
// To run in batch, do
// stressLinear : with no parameters, run standard test with matrices upto 100x100
// stressLinear 30 : run test with matrices upto 30x30
// stressLinear 30 1 : run test with matrices upto 30x30 in verbose mode
//
// To run interactively, do
// root -b
// Root > .x stressLinear.cxx(100) run standard test with matrices upto 100x100
// Root > .x stressLinear.cxx(30) run with matrices upto 30x30
// Root > .x stressLinear.cxx(30,1) run with matrices upto 30x30 in verbose mode
//
// Several tests are run sequentially. Each test will produce one line (Test OK or Test FAILED) .
// At the end of the test a table is printed showing the global results
// Real Time and Cpu Time.
// One single number (ROOTMARKS) is also calculated showing the relative
// performance of your machine compared to a reference machine
// a Pentium IV 2.4 Ghz) with 512 MBytes of memory
// and 120 GBytes IDE disk.
//
// An example of output when all the tests run OK is shown below:
//////////////////////////////////////////////////////////////////////////
// //
// Linear Algebra Package -- Matrix Verifications. //
// //
// This file implements a large set of TMat operation tests. //
// ******************************************************************* //
// * Linear Algebra - S T R E S S suite * //
// ******************************************************************* //
// ******************************************************************* //
// * Starting Matrix - S T R E S S * //
// ******************************************************************* //
// Test 1 : Allocation, Resizing.................................. OK //
// Test 2 : Filling, Inserting, Using............................. OK //
// Test 3 : Uniform matrix operations............................. OK //
// Test 4 : Binary Matrix element-by-element operations............OK //
// Test 5 : Matrix transposition...................................OK //
// Test 6 : Haar/Hilbert Matrix....................................OK //
// Test 7 : Matrix promises........................................OK //
// Test 8 : Matrix Norms...........................................OK //
// Test 9 : Matrix Determinant.....................................OK //
// Test 10 : General Matrix Multiplications.........................OK //
// Test 11 : Symmetric Matrix Multiplications.......................OK //
// Test 12 : Matrix Vector Multiplications..........................OK //
// Test 13 : Matrix Inversion.......................................OK //
// Test 14 : Matrix Persistence.....................................OK //
// ****************************************************************** //
// * Starting Sparse Matrix - S T R E S S * //
// ****************************************************************** //
// Test 1 : Allocation, Resizing...................................OK //
// Test 2 : Filling, Inserting, Using..............................OK //
// Test 3 : Uniform matrix operations..............................OK //
// Test 4 : Binary Matrix element-by-element operations............OK //
// Test 5 : Matrix transposition...................................OK //
// Test 6 : Matrix Norms...........................................OK //
// Test 7 : General Matrix Multiplications.........................OK //
// Test 8 : Matrix Vector Multiplications..........................OK //
// Test 9 : Matrix Slices to Vectors...............................OK //
// Test 10 : Matrix Persistence.....................................OK //
// ******************************************************************* //
// * Starting Vector - S T R E S S * //
// ******************************************************************* //
// Test 1 : Allocation, Filling, Resizing......................... OK //
// Test 2 : Uniform vector operations............................. OK //
// Test 3 : Binary vector element-by-element operations............OK //
// Test 4 : Vector Norms...........................................OK //
// Test 5 : Matrix Slices to Vectors...............................OK //
// Test 6 : Vector Persistence.....................................OK //
// ******************************************************************* //
// * Starting Linear Algebra - S T R E S S * //
// ******************************************************************* //
// Test 1 : Decomposition / Reconstruction........................ OK //
// Test 2 : Linear Equations...................................... OK //
// Test 3 : Pseudo-Inverse, Moore-Penrose......................... OK //
// Test 4 : Eigen - Values/Vectors.................................OK //
// Test 5 : Decomposition Persistence..............................OK //
// ******************************************************************* //
// //
//////////////////////////////////////////////////////////////////////////
#include <stdlib.h>
#include <Riostream.h>
#include <TSystem.h>
#include <TFile.h>
#include <TBenchmark.h>
#include <TArrayD.h>
#include <TF1.h>
#include <TGraph.h>
#include <TROOT.h>
#include "TMath.h"
#include "TMatrixF.h"
#include "TMatrixFSym.h"
#include "TMatrixFSparse.h"
#include "TMatrixFLazy.h"
#include "TVectorF.h"
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixDSparse.h"
#include "TMatrixDLazy.h"
#include "TMatrixDUtils.h"
#include "TVectorD.h"
#include "TDecompLU.h"
#include "TDecompChol.h"
#include "TDecompQRH.h"
#include "TDecompSVD.h"
#include "TDecompBK.h"
#include "TMatrixDEigen.h"
#include "TMatrixDSymEigen.h"
void stressLinear (Int_t maxSizeReq=100,Int_t verbose=0);
void StatusPrint (Int_t id,const TString &title,Bool_t status);
void mstress_allocation (Int_t msize);
void mstress_matrix_fill (Int_t rsize,Int_t csize);
void mstress_element_op (Int_t rsize,Int_t csize);
void mstress_binary_ebe_op (Int_t rsize, Int_t csize);
void mstress_transposition (Int_t msize);
void mstress_special_creation (Int_t dim);
void mstress_matrix_fill (Int_t rsize,Int_t csize);
void mstress_matrix_promises (Int_t dim);
void mstress_norms (Int_t rsize,Int_t csize);
void mstress_determinant (Int_t msize);
void mstress_mm_multiplications ();
void mstress_sym_mm_multiplications(Int_t msize);
void mstress_vm_multiplications ();
void mstress_inversion ();
void mstress_matrix_io ();
void spstress_allocation (Int_t msize);
void spstress_matrix_fill (Int_t rsize,Int_t csize);
void spstress_element_op (Int_t rsize,Int_t csize);
void spstress_binary_ebe_op (Int_t rsize, Int_t csize);
void spstress_transposition (Int_t msize);
void spstress_norms (Int_t rsize,Int_t csize);
void spstress_mm_multiplications ();
void spstress_vm_multiplications ();
void spstress_matrix_slices (Int_t vsize);
void spstress_matrix_io ();
void vstress_allocation (Int_t msize);
void vstress_element_op (Int_t vsize);
void vstress_binary_op (Int_t vsize);
void vstress_norms (Int_t vsize);
void vstress_matrix_slices (Int_t vsize);
void vstress_vector_io ();
Bool_t test_svd_expansion (const TMatrixD &A);
void astress_decomp ();
void astress_lineqn ();
void astress_pseudo ();
void astress_eigen (Int_t msize);
void astress_decomp_io (Int_t msize);
void stress_backward_io ();
void cleanup ();
//_____________________________batch only_____________________
#ifndef __CINT__
int main(int argc,const char *argv[])
{
gROOT->SetBatch();
Int_t maxSizeReq = 100;
if (argc > 1) maxSizeReq = atoi(argv[1]);
Int_t verbose = 0;
if (argc > 2) verbose = (atoi(argv[2]) > 0 ? 1 : 0);
gBenchmark = new TBenchmark();
stressLinear(maxSizeReq,verbose);
cleanup();
return 0;
}
#endif
#define EPSILON 1.0e-14
Int_t gVerbose = 0;
Int_t gNrLoop;
const Int_t nrSize = 20;
const Int_t gSizeA[] = {5,6,7,8,9,10,15,20,30,40,50,60,70,80,90,100,300,500,700,1000};
//________________________________common part_________________________
void stressLinear(Int_t maxSizeReq,Int_t verbose)
{
std::cout << "******************************************************************" <<std::endl;
std::cout << "* Starting Linear Algebra - S T R E S S suite *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
gVerbose = verbose;
gBenchmark->Start("stressLinear");
gNrLoop = nrSize-1;
while (gNrLoop > 0 && maxSizeReq < gSizeA[gNrLoop])
gNrLoop--;
const Int_t maxSize = gSizeA[gNrLoop];
// Matrix
{
std::cout << "* Starting Matrix - S T R E S S *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
mstress_allocation(maxSize);
mstress_matrix_fill(maxSize,maxSize/2);
mstress_element_op(maxSize,maxSize/2);
mstress_binary_ebe_op(maxSize/2,maxSize);
mstress_transposition(maxSize);
mstress_special_creation(maxSize);
#ifndef __CINT__
mstress_matrix_promises(maxSize);
#endif
mstress_norms(maxSize,maxSize/2);
mstress_determinant(maxSize);
mstress_mm_multiplications();
mstress_sym_mm_multiplications(maxSize);
mstress_vm_multiplications();
mstress_inversion();
mstress_matrix_io();
std::cout << "******************************************************************" <<std::endl;
}
// Sparse Matrix
{
std::cout << "* Starting Sparse Matrix - S T R E S S *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
spstress_allocation(maxSize);
spstress_matrix_fill(maxSize,maxSize/2);
spstress_element_op(maxSize,maxSize/2);
spstress_binary_ebe_op(maxSize/2,maxSize);
spstress_transposition(maxSize);
spstress_norms(maxSize,maxSize/2);
spstress_mm_multiplications();
spstress_vm_multiplications();
spstress_matrix_slices(maxSize);
spstress_matrix_io();
std::cout << "******************************************************************" <<std::endl;
}
{
std::cout << "* Starting Vector - S T R E S S *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
vstress_allocation(maxSize);
vstress_element_op(maxSize);
vstress_binary_op(maxSize);
vstress_norms(maxSize);
vstress_matrix_slices(maxSize);
vstress_vector_io();
std::cout << "******************************************************************" <<std::endl;
}
// Linear Algebra
{
std::cout << "* Starting Linear Algebra - S T R E S S *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
astress_decomp();
astress_lineqn();
astress_pseudo();
astress_eigen(5);
astress_decomp_io(10);
std::cout << "******************************************************************" <<std::endl;
}
//Backward Compatibility of Streamers
{
std::cout << "* Starting Backward IO compatibility - S T R E S S *" <<std::endl;
std::cout << "******************************************************************" <<std::endl;
stress_backward_io();
std::cout << "******************************************************************" <<std::endl;
}
gBenchmark->Stop("stressLinear");
//Print table with results
Bool_t UNIX = strcmp(gSystem->GetName(), "Unix") == 0;
printf("******************************************************************\n");
if (UNIX) {
TString sp = gSystem->GetFromPipe("uname -a");
sp.Resize(60);
printf("* SYS: %s\n",sp.Data());
if (strstr(gSystem->GetBuildNode(),"Linux")) {
sp = gSystem->GetFromPipe("lsb_release -d -s");
printf("* SYS: %s\n",sp.Data());
}
if (strstr(gSystem->GetBuildNode(),"Darwin")) {
sp = gSystem->GetFromPipe("sw_vers -productVersion");
sp += " Mac OS X ";
printf("* SYS: %s\n",sp.Data());
}
} else {
const Char_t *os = gSystem->Getenv("OS");
if (!os) printf("* SYS: Windows 95\n");
else printf("* SYS: %s %s \n",os,gSystem->Getenv("PROCESSOR_IDENTIFIER"));
}
printf("******************************************************************\n");
gBenchmark->Print("stressLinear");
const Int_t nr = 7;
const Double_t x_b12[] = { 10., 30., 50., 100., 300., 500., 700.};
const Double_t y_b12[] = {10.74, 15.72, 20.00, 35.79, 98.77, 415.34, 1390.33};
TGraph gr(nr,x_b12,y_b12);
Double_t ct = gBenchmark->GetCpuTime("stressLinear");
const Double_t rootmarks = 600*gr.Eval(maxSize)/ct;
printf("******************************************************************\n");
printf("* ROOTMARKS =%6.1f * Root%-8s %d/%d\n",rootmarks,gROOT->GetVersion(),
gROOT->GetVersionDate(),gROOT->GetVersionTime());
printf("******************************************************************\n");
}
////////////////////////////////////////////////////////////////////////////////
/// Print test program number and its title
void StatusPrint(Int_t id,const TString &title,Bool_t status)
{
const Int_t kMAX = 65;
Char_t number[4];
snprintf(number,4,"%2d",id);
TString header = TString("Test ")+number+" : "+title;
const Int_t nch = header.Length();
for (Int_t i = nch; i < kMAX; i++) header += '.';
std::cout << header << (status ? "OK" : "FAILED") << std::endl;
}
//------------------------------------------------------------------------
// Test allocation functions and compatibility check
//
void mstress_allocation(Int_t msize)
{
if (gVerbose)
std::cout << "\n\n---> Test allocation and compatibility check" << std::endl;
Int_t i,j;
Bool_t ok = kTRUE;
TMatrixD m1(4,msize);
for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
m1(i,j) = TMath::Pi()*i+TMath::E()*j;
TMatrixD m2(0,3,0,msize-1);
TMatrixD m3(1,4,0,msize-1);
TMatrixD m4(m1);
if (gVerbose) {
std::cout << "\nStatus information reported for matrix m3:" << std::endl;
std::cout << " Row lower bound ... " << m3.GetRowLwb() << std::endl;
std::cout << " Row upper bound ... " << m3.GetRowUpb() << std::endl;
std::cout << " Col lower bound ... " << m3.GetColLwb() << std::endl;
std::cout << " Col upper bound ... " << m3.GetColUpb() << std::endl;
std::cout << " No. rows ..........." << m3.GetNrows() << std::endl;
std::cout << " No. cols ..........." << m3.GetNcols() << std::endl;
std::cout << " No. of elements ...." << m3.GetNoElements() << std::endl;
}
if (gVerbose)
std::cout << "\nCheck matrices 1 & 2 for compatibility" << std::endl;
ok &= AreCompatible(m1,m2,gVerbose);
if (gVerbose)
std::cout << "Check matrices 1 & 4 for compatibility" << std::endl;
ok &= AreCompatible(m1,m4,gVerbose);
if (gVerbose)
std::cout << "m2 has to be compatible with m3 after resizing to m3" << std::endl;
m2.ResizeTo(m3);
ok &= AreCompatible(m2,m3,gVerbose);
TMatrixD m5(m1.GetNrows()+1,m1.GetNcols()+5);
for (i = m5.GetRowLwb(); i <= m5.GetRowUpb(); i++)
for (j = m5.GetColLwb(); j <= m5.GetColUpb(); j++)
m5(i,j) = TMath::Pi()*i+TMath::E()*j;
if (gVerbose)
std::cout << "m1 has to be compatible with m5 after resizing to m5" << std::endl;
m1.ResizeTo(m5.GetNrows(),m5.GetNcols());
ok &= AreCompatible(m1,m5,gVerbose);
if (gVerbose)
std::cout << "m1 has to be equal to m4 after stretching and shrinking" << std::endl;
m1.ResizeTo(m4.GetNrows(),m4.GetNcols());
ok &= VerifyMatrixIdentity(m1,m4,gVerbose,EPSILON);
if (gVerbose)
std::cout << "m5 has to be equal to m1 after shrinking" << std::endl;
m5.ResizeTo(m1.GetNrows(),m1.GetNcols());
ok &= VerifyMatrixIdentity(m1,m5,gVerbose,EPSILON);
if (gVerbose)
std::cout << "stretching and shrinking for small matrices (stack)" << std::endl;
if (gVerbose)
std::cout << "m8 has to be equal to m7 after stretching and shrinking" << std::endl;
TMatrixD m6(4,4);
for (i = m6.GetRowLwb(); i <= m6.GetRowUpb(); i++)
for (j = m6.GetColLwb(); j <= m6.GetColUpb(); j++)
m6(i,j) = TMath::Pi()*i+TMath::E()*j;
TMatrixD m8(3,3);
for (i = m8.GetRowLwb(); i <= m8.GetRowUpb(); i++)
for (j = m8.GetColLwb(); j <= m8.GetColUpb(); j++)
m8(i,j) = TMath::Pi()*i+TMath::E()*j;
TMatrixD m7(m8);
m8.ResizeTo(4,4);
m8.ResizeTo(3,3);
ok &= VerifyMatrixIdentity(m7,m8,gVerbose,EPSILON);
if (gVerbose)
std::cout << "m6 has to be equal to m8 after shrinking" << std::endl;
m6.ResizeTo(3,3);
ok &= VerifyMatrixIdentity(m6,m8,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(1,"Allocation, Resizing",ok);
}
class FillMatrix : public TElementPosActionD {
Int_t no_elems,no_cols;
void Operation(Double_t &element) const
{ element = 4*TMath::Pi()/no_elems * (fI*no_cols+fJ); }
public:
FillMatrix() {}
FillMatrix(const TMatrixD &m) :
no_elems(m.GetNoElements()),no_cols(m.GetNcols()) { }
};
//
//------------------------------------------------------------------------
// Test Filling of matrix
//
void mstress_matrix_fill(Int_t rsize,Int_t csize)
{
if (gVerbose)
std::cout << "\n\n---> Test different matrix filling methods\n" << std::endl;
Bool_t ok = kTRUE;
if (gVerbose)
std::cout << "Creating m with Apply function..." << std::endl;
TMatrixD m(-1,rsize-2,1,csize);
#ifndef __CINT__
FillMatrix f(m);
m.Apply(f);
#else
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = 4*TMath::Pi()/m.GetNoElements() * (i*m.GetNcols()+j);
#endif
{
if (gVerbose)
std::cout << "Check identity between m and matrix filled through (i,j)" << std::endl;
TMatrixD m_overload1(-1,rsize-2,1,csize);
TMatrixD m_overload2(-1,rsize-2,1,csize);
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
{
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
{
const Double_t val = 4*TMath::Pi()/rsize/csize*(i*csize+j);
m_overload1(i,j) = val;
m_overload2[i][j] = val;
}
}
ok &= VerifyMatrixIdentity(m,m_overload1,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Check identity between m and matrix filled through [i][j]" << std::endl;
ok &= VerifyMatrixIdentity(m,m_overload2,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Check identity between matrix filled through [i][j] and (i,j)" << std::endl;
ok &= VerifyMatrixIdentity(m_overload1,m_overload2,gVerbose,EPSILON);
}
{
TArrayD a_fortran(rsize*csize);
TArrayD a_c (rsize*csize);
for (Int_t i = 0; i < rsize; i++)
{
for (Int_t j = 0; j < csize; j++)
{
a_c[i*csize+j] = 4*TMath::Pi()/rsize/csize*((i-1)*csize+j+1);
a_fortran[i+rsize*j] = a_c[i*csize+j];
}
}
if (gVerbose)
std::cout << "Creating m_fortran by filling with fortran stored matrix" << std::endl;
TMatrixD m_fortran(-1,rsize-2,1,csize,a_fortran.GetArray(),"F");
if (gVerbose)
std::cout << "Check identity between m and m_fortran" << std::endl;
ok &= VerifyMatrixIdentity(m,m_fortran,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Creating m_c by filling with c stored matrix" << std::endl;
TMatrixD m_c(-1,rsize-2,1,csize,a_c.GetArray());
if (gVerbose)
std::cout << "Check identity between m and m_c" << std::endl;
ok &= VerifyMatrixIdentity(m,m_c,gVerbose,EPSILON);
}
{
if (gVerbose)
std::cout << "Check insertion/extraction of sub-matrices" << std::endl;
{
TMatrixD m_sub1 = m;
m_sub1.ResizeTo(0,rsize-2,2,csize);
TMatrixD m_sub2 = m.GetSub(0,rsize-2,2,csize,"");
ok &= VerifyMatrixIdentity(m_sub1,m_sub2,gVerbose,EPSILON);
}
{
TMatrixD m2(-1,rsize-2,1,csize);
TMatrixD m_part1 = m.GetSub(0,rsize-2,2,csize,"");
TMatrixD m_part2 = m.GetSub(0,rsize-2,1,1,"");
TMatrixD m_part3 = m.GetSub(-1,-1,2,csize,"");
TMatrixD m_part4 = m.GetSub(-1,-1,1,1,"");
m2.SetSub(0,2,m_part1);
m2.SetSub(0,1,m_part2);
m2.SetSub(-1,2,m_part3);
m2.SetSub(-1,1,m_part4);
ok &= VerifyMatrixIdentity(m,m2,gVerbose,EPSILON);
}
{
TMatrixD m2(-1,rsize-2,1,csize);
TMatrixD m_part1 = m.GetSub(0,rsize-2,2,csize,"S");
TMatrixD m_part2 = m.GetSub(0,rsize-2,1,1,"S");
TMatrixD m_part3 = m.GetSub(-1,-1,2,csize,"S");
TMatrixD m_part4 = m.GetSub(-1,-1,1,1,"S");
m2.SetSub(0,2,m_part1);
m2.SetSub(0,1,m_part2);
m2.SetSub(-1,2,m_part3);
m2.SetSub(-1,1,m_part4);
ok &= VerifyMatrixIdentity(m,m2,gVerbose,EPSILON);
}
}
{
if (gVerbose)
std::cout << "Check sub-matrix views" << std::endl;
{
TMatrixD m3(-1,rsize-2,1,csize);
TMatrixDSub(m3,0,rsize-2,2,csize) = TMatrixDSub(m,0,rsize-2,2,csize);
TMatrixDSub(m3,0,rsize-2,1,1) = TMatrixDSub(m,0,rsize-2,1,1);
TMatrixDSub(m3,-1,-1,2,csize) = TMatrixDSub(m,-1,-1,2,csize);
TMatrixDSub(m3,-1,-1,1,1) = TMatrixDSub(m,-1,-1,1,1);
ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
TMatrixD unit(3,3);
TMatrixDSub(m3,1,3,1,3) = unit.UnitMatrix();
TMatrixDSub(m3,1,3,1,3) *= m.GetSub(1,3,1,3);
ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
TMatrixDSub(m3,0,rsize-2,2,csize) = 1.0;
TMatrixDSub(m3,0,rsize-2,1,1) = 1.0;
TMatrixDSub(m3,-1,-1,2,csize) = 1.0;
TMatrixDSub(m3,-1,-1,1,1) = 1.0;
ok &= (m3 == 1.0);
}
}
{
if (gVerbose)
std::cout << "Check array Use" << std::endl;
{
TMatrixD *m1a = new TMatrixD(m);
TMatrixD *m2a = new TMatrixD();
m2a->Use(m1a->GetRowLwb(),m1a->GetRowUpb(),m1a->GetColLwb(),m1a->GetColUpb(),m1a->GetMatrixArray());
ok &= VerifyMatrixIdentity(m,*m2a,gVerbose,EPSILON);
m2a->Sqr();
TMatrixD m4 = m; m4.Sqr();
ok &= VerifyMatrixIdentity(m4,*m1a,gVerbose,EPSILON);
delete m1a;
delete m2a;
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(2,"Filling, Inserting, Using",ok);
}
//
//------------------------------------------------------------------------
// Test uniform element operations
//
typedef Double_t (*dfunc)(Double_t);
class ApplyFunction : public TElementActionD {
dfunc fFunc;
void Operation(Double_t &element) const { element = fFunc(Double_t(element)); }
public:
ApplyFunction(dfunc func) : fFunc(func) { }
};
void mstress_element_op(Int_t rsize,Int_t csize)
{
Bool_t ok = kTRUE;
const Double_t pattern = 8.625;
TMatrixD m(-1,rsize-2,1,csize);
if (gVerbose)
std::cout << "\nWriting zeros to m..." << std::endl;
{
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for(Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = 0;
ok &= VerifyMatrixValue(m,0.,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "Creating zero m1 ..." << std::endl;
TMatrixD m1(TMatrixD::kZero, m);
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing m1 with 0 ..." << std::endl;
R__ASSERT(m1 == 0);
R__ASSERT(!(m1 != 0));
if (gVerbose)
std::cout << "Writing a pattern " << pattern << " by assigning to m(i,j)..." << std::endl;
{
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = pattern;
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "Writing the pattern by assigning to m1 as a whole ..." << std::endl;
m1 = pattern;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing m and m1 ..." << std::endl;
R__ASSERT(m == m1);
if (gVerbose)
std::cout << "Comparing (m=0) and m1 ..." << std::endl;
R__ASSERT(!(m.Zero() == m1));
if (gVerbose)
std::cout << "Clearing m1 ..." << std::endl;
m1.Zero();
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nClear m and add the pattern" << std::endl;
m.Zero();
m += pattern;
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " add the doubled pattern with the negative sign" << std::endl;
m += -2*pattern;
ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtract the trippled pattern with the negative sign" << std::endl;
m -= -3*pattern;
ok &= VerifyMatrixValue(m,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify comparison operations when all elems are the same" << std::endl;
m = pattern;
R__ASSERT( m == pattern && !(m != pattern) );
R__ASSERT( m > 0 && m >= pattern && m <= pattern );
R__ASSERT( m > -pattern && m >= -pattern );
R__ASSERT( m <= pattern && !(m < pattern) );
m -= 2*pattern;
R__ASSERT( m < -pattern/2 && m <= -pattern/2 );
R__ASSERT( m >= -pattern && !(m > -pattern) );
if (gVerbose)
std::cout << "\nVerify comparison operations when not all elems are the same" << std::endl;
m = pattern; m(m.GetRowUpb(),m.GetColUpb()) = pattern-1;
R__ASSERT( !(m == pattern) && !(m != pattern) );
R__ASSERT( m != 0 ); // none of elements are 0
R__ASSERT( !(m >= pattern) && m <= pattern && !(m<pattern) );
R__ASSERT( !(m <= pattern-1) && m >= pattern-1 && !(m>pattern-1) );
if (gVerbose)
std::cout << "\nAssign 2*pattern to m by repeating additions" << std::endl;
m = 0; m += pattern; m += pattern;
if (gVerbose)
std::cout << "Assign 2*pattern to m1 by multiplying by two " << std::endl;
m1 = pattern; m1 *= 2;
ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
R__ASSERT( m == m1 );
if (gVerbose)
std::cout << "Multiply m1 by one half returning it to the 1*pattern" << std::endl;
m1 *= 1/2.;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nAssign -pattern to m and m1" << std::endl;
m.Zero(); m -= pattern; m1 = -pattern;
ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
R__ASSERT( m == m1 );
if (gVerbose)
std::cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same" << std::endl;
m.Sqr();
ok &= VerifyMatrixValue(m,pattern*pattern,gVerbose,EPSILON);
m.Sqrt();
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
m1.Abs();
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1" << std::endl;
{
#ifndef __CINT__
FillMatrix f(m);
m.Apply(f);
#else
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = 4*TMath::Pi()/m.GetNoElements() * (i*m.GetNcols()+j);
#endif
}
m1 = m;
{
#ifndef __CINT__
ApplyFunction s(&TMath::Sin);
ApplyFunction c(&TMath::Cos);
m.Apply(s);
m1.Apply(c);
#else
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
m(i,j) = TMath::Sin(m(i,j));
m1(i,j) = TMath::Cos(m1(i,j));
}
}
#endif
}
m.Sqr();
m1.Sqr();
m += m1;
ok &= VerifyMatrixValue(m,1.,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(3,"Uniform matrix operations",ok);
}
//
//------------------------------------------------------------------------
// Test binary matrix element-by-element operations
//
void mstress_binary_ebe_op(Int_t rsize, Int_t csize)
{
if (gVerbose)
std::cout << "\n---> Test Binary Matrix element-by-element operations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 4.25;
TMatrixD m(2,rsize+1,0,csize-1);
TMatrixD m1(TMatrixD::kZero,m);
TMatrixD mp(TMatrixD::kZero,m);
{
for (Int_t i = mp.GetRowLwb(); i <= mp.GetRowUpb(); i++)
for (Int_t j = mp.GetColLwb(); j <= mp.GetColUpb(); j++)
mp(i,j) = (i-m.GetNrows()/2.)*j*pattern;
}
if (gVerbose)
std::cout << "\nVerify assignment of a matrix to the matrix" << std::endl;
m = pattern;
m1.Zero();
m1 = m;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
R__ASSERT( m1 == m );
if (gVerbose)
std::cout << "\nAdding the matrix to itself, uniform pattern " << pattern << std::endl;
m.Zero(); m = pattern;
m1 = m; m1 += m1;
ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting two matrices ..." << std::endl;
m1 -= m;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting the matrix from itself" << std::endl;
m1 -= m1;
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << " adding two matrices together" << std::endl;
m1 += m;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose) {
std::cout << "\nArithmetic operations on matrices with not the same elements" << std::endl;
std::cout << " adding mp to the zero matrix..." << std::endl;
}
m.Zero(); m += mp;
ok &= VerifyMatrixIdentity(m,mp,gVerbose,EPSILON);
m1 = m;
if (gVerbose)
std::cout << " making m = 3*mp and m1 = 3*mp, via add() and succesive mult" << std::endl;
Add(m,2.,mp);
m1 += m1; m1 += mp;
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose)
std::cout << " clear both m and m1, by subtracting from itself and via add()" << std::endl;
m1 -= m1;
Add(m,-3.,mp);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose) {
std::cout << "\nTesting element-by-element multiplications and divisions" << std::endl;
std::cout << " squaring each element with sqr() and via multiplication" << std::endl;
}
m = mp; m1 = mp;
m.Sqr();
ElementMult(m1,m1);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose)
std::cout << " compare (m = pattern^2)/pattern with pattern" << std::endl;
m = pattern; m1 = pattern;
m.Sqr();
ElementDiv(m,m1);
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
if (gVerbose)
Compare(m1,m);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(4,"Binary Matrix element-by-element operations",ok);
}
//
//------------------------------------------------------------------------
// Verify matrix transposition
//
void mstress_transposition(Int_t msize)
{
if (gVerbose) {
std::cout << "\n---> Verify matrix transpose "
"for matrices of a characteristic size " << msize << std::endl;
}
Bool_t ok = kTRUE;
{
if (gVerbose)
std::cout << "\nCheck to see that a square UnitMatrix stays the same";
TMatrixD m(msize,msize);
m.UnitMatrix();
TMatrixD mt(TMatrixD::kTransposed,m);
ok &= ( m == mt ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nTest a non-square UnitMatrix";
TMatrixD m(msize,msize+1);
m.UnitMatrix();
TMatrixD mt(TMatrixD::kTransposed,m);
R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows() );
for (Int_t i = m.GetRowLwb(); i <= TMath::Min(m.GetRowUpb(),m.GetColUpb()); i++)
for (Int_t j = m.GetColLwb(); j <= TMath::Min(m.GetRowUpb(),m.GetColUpb()); j++)
ok &= ( m(i,j) == mt(i,j) ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nCheck to see that a symmetric (Hilbert)Matrix stays the same";
TMatrixD m = THilbertMatrixD(msize,msize);
TMatrixD mt(TMatrixD::kTransposed,m);
ok &= ( m == mt ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nCheck transposing a non-symmetric matrix";
TMatrixD m = THilbertMatrixD(msize+1,msize);
m(1,2) = TMath::Pi();
TMatrixD mt(TMatrixD::kTransposed,m);
R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows());
R__ASSERT(mt(2,1) == (Double_t)TMath::Pi() && mt(1,2) != (Double_t)TMath::Pi());
R__ASSERT(mt[2][1] == (Double_t)TMath::Pi() && mt[1][2] != (Double_t)TMath::Pi());
if (gVerbose)
std::cout << "\nCheck double transposing a non-symmetric matrix" << std::endl;
TMatrixD mtt(TMatrixD::kTransposed,mt);
ok &= ( m == mtt ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(5,"Matrix transposition",ok);
}
//
//------------------------------------------------------------------------
// Test special matrix creation
//
class MakeHilbert : public TElementPosActionD {
void Operation(Double_t &element) const { element = 1./(fI+fJ+1); }
public:
MakeHilbert() { }
};
#if !defined (__CINT__) || defined (__MAKECINT__)
class TestUnit : public TElementPosActionD {
mutable Int_t fIsUnit;
void Operation(Double_t &element) const
{ if (fIsUnit)
fIsUnit = ((fI==fJ) ? (element == 1.0) : (element == 0)); }
public:
TestUnit() : fIsUnit(0==0) { }
Int_t is_indeed_unit() const { return fIsUnit; }
};
#else
Bool_t is_indeed_unit(TMatrixD &m) {
Bool_t isUnit = kTRUE;
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
if (isUnit)
isUnit = ((i==j) ? (m(i,j) == 1.0) : (m(i,j) == 0));
}
return isUnit;
}
#endif
void mstress_special_creation(Int_t dim)
{
if (gVerbose)
std::cout << "\n---> Check creating some special matrices of dimension " << dim << std::endl;
Bool_t ok = kTRUE;
{
if (gVerbose)
std::cout << "\ntest creating Hilbert matrices" << std::endl;
TMatrixD m = THilbertMatrixD(dim+1,dim);
TMatrixD m1(TMatrixD::kZero,m);
ok &= ( !(m == m1) ) ? kTRUE : kFALSE;
ok &= ( m != 0 ) ? kTRUE : kFALSE;
#ifndef __CINT__
MakeHilbert mh;
m1.Apply(mh);
#else
for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
m1(i,j) = 1./(i+j+1);
#endif
ok &= ( m1 != 0 ) ? kTRUE : kFALSE;
ok &= ( m == m1 ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "test creating zero matrix and copy constructor" << std::endl;
TMatrixD m = THilbertMatrixD(dim,dim+1);
ok &= ( m != 0 ) ? kTRUE : kFALSE;
TMatrixD m1(m); // Applying the copy constructor
ok &= ( m1 == m ) ? kTRUE : kFALSE;
TMatrixD m2(TMatrixD::kZero,m);
ok &= ( m2 == 0 ) ? kTRUE : kFALSE;
ok &= ( m != 0 ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "test creating unit matrices" << std::endl;
TMatrixD m(dim,dim);
#ifndef __CINT__
{
TestUnit test_unit;
m.Apply(test_unit);
ok &= ( !test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
}
#else
ok &= ( !is_indeed_unit(m) ) ? kTRUE : kFALSE;
#endif
m.UnitMatrix();
#ifndef __CINT__
{
TestUnit test_unit;
m.Apply(test_unit);
ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
}
#else
ok &= ( is_indeed_unit(m) ) ? kTRUE : kFALSE;
#endif
m.ResizeTo(dim-1,dim);
TMatrixD m2(TMatrixD::kUnit,m);
#ifndef __CINT__
{
TestUnit test_unit;
m2.Apply(test_unit);
ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
}
#else
ok &= ( is_indeed_unit(m2) ) ? kTRUE : kFALSE;
#endif
m.ResizeTo(dim,dim-2);
m.UnitMatrix();
#ifndef __CINT__
{
TestUnit test_unit;
m.Apply(test_unit);
ok &= ( test_unit.is_indeed_unit() ) ? kTRUE : kFALSE;
}
#else
ok &= ( is_indeed_unit(m) ) ? kTRUE : kFALSE;
#endif
}
{
if (gVerbose)
std::cout << "check to see that Haar matrix has *exactly* orthogonal columns" << std::endl;
Int_t j;
const Int_t order = 5;
const TMatrixD haar = THaarMatrixD(order);
ok &= ( haar.GetNrows() == (1<<order) &&
haar.GetNrows() == haar.GetNcols() ) ? kTRUE : kFALSE;
TVectorD colj(1<<order);
TVectorD coll(1<<order);
for (j = haar.GetColLwb(); j <= haar.GetColUpb(); j++) {
colj = TMatrixDColumn_const(haar,j);
ok &= (TMath::Abs(colj*colj-1.0) <= 1.0e-15 ) ? kTRUE : kFALSE;
for (Int_t l = j+1; l <= haar.GetColUpb(); l++) {
coll = TMatrixDColumn_const(haar,l);
const Double_t val = colj*coll;
ok &= ( TMath::Abs(val) <= 1.0e-15 ) ? kTRUE : kFALSE;
}
}
if (gVerbose)
std::cout << "make Haar (sub)matrix and test it *is* a submatrix" << std::endl;
const Int_t no_sub_cols = (1<<order) - 3;
const TMatrixD haar_sub = THaarMatrixD(order,no_sub_cols);
ok &= ( haar_sub.GetNrows() == (1<<order) &&
haar_sub.GetNcols() == no_sub_cols ) ? kTRUE : kFALSE;
for (j = haar_sub.GetColLwb(); j <= haar_sub.GetColUpb(); j++) {
colj = TMatrixDColumn_const(haar,j);
coll = TMatrixDColumn_const(haar_sub,j);
ok &= VerifyVectorIdentity(colj,coll,gVerbose);
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(6,"Haar/Hilbert Matrix",ok);
}
#ifndef __CINT__
//
//------------------------------------------------------------------------
// Test matrix promises
//
class hilbert_matrix_promise : public TMatrixDLazy {
void FillIn(TMatrixD &m) const { m = THilbertMatrixD(m.GetRowLwb(),m.GetRowUpb(),
m.GetColLwb(),m.GetColUpb()); }
public:
hilbert_matrix_promise(Int_t nrows,Int_t ncols)
: TMatrixDLazy(nrows,ncols) {}
hilbert_matrix_promise(Int_t row_lwb,Int_t row_upb,
Int_t col_lwb,Int_t col_upb)
: TMatrixDLazy(row_lwb,row_upb,col_lwb,col_upb) { }
};
void mstress_matrix_promises(Int_t dim)
{
if (gVerbose)
std::cout << "\n---> Check making/forcing promises, (lazy)matrices of dimension " << dim << std::endl;
Bool_t ok = kTRUE;
{
if (gVerbose)
std::cout << "\nmake a promise and force it by a constructor" << std::endl;
TMatrixD m = hilbert_matrix_promise(dim,dim+1);
TMatrixD m1 = THilbertMatrixD(dim,dim+1);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
}
{
if (gVerbose)
std::cout << "make a promise and force it by an assignment" << std::endl;
TMatrixD m(-1,dim,0,dim);
m = hilbert_matrix_promise(-1,dim,0,dim);
TMatrixD m1 = THilbertMatrixD(-1,dim,0,dim);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(7,"Matrix promises",ok);
}
#endif
//
//------------------------------------------------------------------------
// Verify the norm calculation
//
void mstress_norms(Int_t rsize_req,Int_t csize_req)
{
if (gVerbose)
std::cout << "\n---> Verify norm calculations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 10.25;
Int_t rsize = rsize_req;
if (rsize%2 != 0) rsize--;
Int_t csize = csize_req;
if (csize%2 != 0) csize--;
if (rsize%2 == 1 || csize%2 == 1) {
std::cout << "rsize: " << rsize <<std::endl;
std::cout << "csize: " << csize <<std::endl;
Fatal("mstress_norms","Sorry, size of the matrix to test must be even for this test\n");
}
TMatrixD m(rsize,csize);
if (gVerbose)
std::cout << "\nAssign " << pattern << " to all the elements and check norms" << std::endl;
m = pattern;
if (gVerbose)
std::cout << " 1. (col) norm should be pattern*nrows" << std::endl;
ok &= ( m.Norm1() == pattern*m.GetNrows() ) ? kTRUE : kFALSE;
ok &= ( m.Norm1() == m.ColNorm() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Inf (row) norm should be pattern*ncols" << std::endl;
ok &= ( m.NormInf() == pattern*m.GetNcols() ) ? kTRUE : kFALSE;
ok &= ( m.NormInf() == m.RowNorm() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Square of the Eucl norm has got to be pattern^2 * no_elems" << std::endl;
ok &= ( m.E2Norm() == (pattern*pattern)*m.GetNoElements() ) ? kTRUE : kFALSE;
TMatrixD m1(TMatrixD::kZero,m);
ok &= ( m.E2Norm() == E2Norm(m,m1) ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(8,"Matrix Norms",ok);
}
//
//------------------------------------------------------------------------
// Verify the determinant evaluation
//
void mstress_determinant(Int_t msize)
{
if (gVerbose)
std::cout << "\n---> Verify determinant evaluation for a square matrix of size " << msize << std::endl;
Bool_t ok = kTRUE;
TMatrixD m(msize,msize);
const Double_t pattern = 2.5;
{
if (gVerbose)
std::cout << "\nCheck to see that the determinant of the unit matrix is one" <<std::endl;
m.UnitMatrix();
Double_t d1,d2;
m.Determinant(d1,d2);
const Double_t det = d1*TMath::Power(2.,d2);
if (gVerbose) {
std::cout << "det = " << det << " deviation= " << TMath::Abs(det-1);
std::cout << ( (TMath::Abs(det-1) < EPSILON) ? " OK" : " too large") <<std::endl;
}
ok &= (TMath::Abs(det-1) < EPSILON) ? kTRUE : kFALSE;
}
if (ok)
{
if (gVerbose)
std::cout << "\nCheck the determinant for the matrix with " << pattern << " at the diagonal" << std::endl;
{
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = ( i==j ? pattern : 0 );
}
Double_t d1,d2;
m.Determinant(d1,d2);
const Double_t d1_abs = TMath::Abs(d1);
const Double_t log_det1 = d2*TMath::Log(2.)+TMath::Log(d1_abs);
const Double_t log_det2 = ((Double_t)m.GetNrows())*TMath::Log(pattern);
if (gVerbose) {
std::cout << "log_det1 = " << log_det1 << " deviation= " << TMath::Abs(log_det1-log_det2);
std::cout << ( (TMath::Abs(log_det1-log_det2) < msize*EPSILON) ? " OK" : " too large") <<std::endl;
}
ok &= ( TMath::Abs(log_det1-log_det2) < msize*EPSILON ) ? kTRUE : kFALSE;
}
if (ok)
{
if (gVerbose)
std::cout << "\nCheck the determinant of the transposed matrix" << std::endl;
m.UnitMatrix();
m(1,2) = pattern;
TMatrixD m_tran(TMatrixD::kTransposed,m);
ok &= ( !(m == m_tran) ) ? kTRUE : kFALSE;
const Double_t det1 = m.Determinant();
const Double_t det2 = m_tran.Determinant();
if (gVerbose) {
std::cout << "det1 = " << det1 << " deviation= " << TMath::Abs(det1-det2);
std::cout << ( (TMath::Abs(det1-det2) < msize*EPSILON) ? " OK" : " too large") <<std::endl;
}
ok &= ( TMath::Abs(det1-det2) < msize*EPSILON ) ? kTRUE : kFALSE;
}
if (ok)
{
if (gVerbose)
std::cout << "\nswap two rows/cols of a matrix through method 1 and watch det's sign" << std::endl;
m.UnitMatrix();
TMatrixDRow(m,3) = pattern;
Double_t d1,d2;
m.Determinant(d1,d2);
TMatrixDRow row1(m,1);
TVectorD vrow1(m.GetRowLwb(),m.GetRowUpb()); vrow1 = row1;
TVectorD vrow3(m.GetRowLwb(),m.GetRowUpb()); vrow3 = TMatrixDRow(m,3);
row1 = vrow3; TMatrixDRow(m,3) = vrow1;
Double_t d1_p,d2_p;
m.Determinant(d1_p,d2_p);
if (gVerbose) {
std::cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1+d1_p);
std::cout << ( ( d1 == -d1_p ) ? " OK" : " too large") <<std::endl;
}
ok &= ( d1 == -d1_p ) ? kTRUE : kFALSE;
TMatrixDColumn col2(m,2);
TVectorD vcol2(m.GetRowLwb(),m.GetRowUpb()); vcol2 = col2;
TVectorD vcol4(m.GetRowLwb(),m.GetRowUpb()); vcol4 = TMatrixDColumn(m,4);
col2 = vcol4; TMatrixDColumn(m,4) = vcol2;
m.Determinant(d1_p,d2_p);
if (gVerbose) {
std::cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1-d1_p);
std::cout << ( ( d1 == d1_p ) ? " OK" : " too large") <<std::endl;
}
ok &= ( d1 == d1_p ) ? kTRUE : kFALSE;
}
if (ok)
{
if (gVerbose)
std::cout << "\nswap two rows/cols of a matrix through method 2 and watch det's sign" << std::endl;
m.UnitMatrix();
TMatrixDRow(m,3) = pattern;
Double_t d1,d2;
m.Determinant(d1,d2);
TMatrixD m_save( m);
TMatrixDRow(m,1) = TMatrixDRow(m_save,3);
TMatrixDRow(m,3) = TMatrixDRow(m_save,1);
Double_t d1_p,d2_p;
m.Determinant(d1_p,d2_p);
if (gVerbose) {
std::cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1+d1_p);
std::cout << ( ( d1 == -d1_p ) ? " OK" : " too large") <<std::endl;
}
ok &= ( d1 == -d1_p ) ? kTRUE : kFALSE;
m_save = m;
TMatrixDColumn(m,2) = TMatrixDColumn(m_save,4);
TMatrixDColumn(m,4) = TMatrixDColumn(m_save,2);
m.Determinant(d1_p,d2_p);
if (gVerbose) {
std::cout << "d1 = " << d1 << " deviation= " << TMath::Abs(d1-d1_p);
std::cout << ( ( d1 == d1_p ) ? " OK" : " too large") <<std::endl;
}
ok &= ( d1 == d1_p ) ? kTRUE : kFALSE;
}
if (ok)
{
if (gVerbose)
std::cout << "\nCheck the determinant for the matrix with " << pattern << " at the anti-diagonal" << std::endl;
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = ( i==(m.GetColUpb()+m.GetColLwb()-j) ? pattern : 0 );
Double_t d1,d2;
m.Determinant(d1,d2);
const Double_t d1_abs = TMath::Abs(d1);
const Double_t log_det1 = d2*TMath::Log(2.)+TMath::Log(d1_abs);
const Double_t log_det2 = ((Double_t)m.GetNrows())*TMath::Log(pattern);
const Double_t sign = ( m.GetNrows()*(m.GetNrows()-1)/2 & 1 ? -1 : 1 );
if (gVerbose) {
std::cout << "log_det1 = " << log_det1 << " deviation= " << TMath::Abs(log_det1-log_det2);
std::cout << ( (TMath::Abs(log_det1-log_det2) < msize*EPSILON) ? " OK" : " too large") <<std::endl;
std::cout << ( sign * d1 > 0. ? " OK sign" : " wrong sign") <<std::endl;
}
ok &= ( TMath::Abs(log_det1-log_det2) < msize*EPSILON ) ? kTRUE : kFALSE;
ok &= ( sign * d1 > 0. ) ? kTRUE : kFALSE;
}
// Next test disabled because it produces (of course) a Warning
if (0) {
if (gVerbose)
std::cout << "\nCheck the determinant for the singular matrix"
"\n\tdefined as above with zero first row" << std::endl;
m.Zero();
{
for (Int_t i = m.GetRowLwb()+1; i <= m.GetRowUpb(); i++)
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++)
m(i,j) = ( i==(m.GetColUpb()+m.GetColLwb()-j) ? pattern : 0 );
}
std::cout << "\n\tdeterminant is " << m.Determinant();
ok &= ( m.Determinant() == 0 ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nCheck out the determinant of the Hilbert matrix";
TMatrixDSym H = THilbertMatrixDSym(3);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 3x3 Hilbert matrix: exact determinant 1/2160 ";
std::cout << "\n computed 1/"<< 1/H.Determinant();
}
H.ResizeTo(4,4);
H = THilbertMatrixDSym(4);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 4x4 Hilbert matrix: exact determinant 1/6048000 ";
std::cout << "\n computed 1/"<< 1/H.Determinant();
}
H.ResizeTo(5,5);
H = THilbertMatrixDSym(5);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 5x5 Hilbert matrix: exact determinant 3.749295e-12";
std::cout << "\n computed "<< H.Determinant();
}
if (gVerbose) {
TDecompQRH qrh(H);
Double_t d1,d2;
qrh.Det(d1,d2);
std::cout << "\n qrh det = " << d1*TMath::Power(2.0,d2) <<std::endl;
}
if (gVerbose) {
TDecompChol chol(H);
Double_t d1,d2;
chol.Det(d1,d2);
std::cout << "\n chol det = " << d1*TMath::Power(2.0,d2) <<std::endl;
}
if (gVerbose) {
TDecompSVD svd(H);
Double_t d1,d2;
svd.Det(d1,d2);
std::cout << "\n svd det = " << d1*TMath::Power(2.0,d2) <<std::endl;
}
H.ResizeTo(7,7);
H = THilbertMatrixDSym(7);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 7x7 Hilbert matrix: exact determinant 4.8358e-25";
std::cout << "\n computed "<< H.Determinant();
}
H.ResizeTo(9,9);
H = THilbertMatrixDSym(9);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 9x9 Hilbert matrix: exact determinant 9.72023e-43";
std::cout << "\n computed "<< H.Determinant();
}
H.ResizeTo(10,10);
H = THilbertMatrixDSym(10);
H.SetTol(1.0e-20);
if (gVerbose) {
std::cout << "\n 10x10 Hilbert matrix: exact determinant 2.16418e-53";
std::cout << "\n computed "<< H.Determinant();
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(9,"Matrix Determinant",ok);
}
//
//------------------------------------------------------------------------
// Verify matrix multiplications
//
void mstress_mm_multiplications()
{
Bool_t ok = kTRUE;
Int_t iloop = 0;
Int_t nr = 0;
while (iloop <= gNrLoop) {
const Int_t msize = gSizeA[iloop];
const Double_t epsilon = EPSILON*msize/100;
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
Int_t i,j;
if (msize <= 5) {
iloop++;
continue; // some references to m(3,..)
}
if (verbose)
std::cout << "\n---> Verify matrix multiplications "
"for matrices of the characteristic size " << msize << std::endl;
{
if (verbose)
std::cout << "\nTest inline multiplications of the UnitMatrix" << std::endl;
TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
TMatrixD u(TMatrixD::kUnit,m);
m(3,1) = TMath::Pi();
u *= m;
ok &= VerifyMatrixIdentity(u,m,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test inline multiplications by a DiagMat" << std::endl;
TMatrixD m = THilbertMatrixD(msize+3,msize);
m(1,3) = TMath::Pi();
TVectorD v(msize);
for (i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = 1+i;
TMatrixD diag(msize,msize);
TMatrixDDiag d(diag);
d = v;
TMatrixD eth = m;
for (i = eth.GetRowLwb(); i <= eth.GetRowUpb(); i++)
for (j = eth.GetColLwb(); j <= eth.GetColUpb(); j++)
eth(i,j) *= v(j);
m *= diag;
ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test XPP = X where P is a permutation matrix" << std::endl;
TMatrixD m = THilbertMatrixD(msize-1,msize);
m(2,3) = TMath::Pi();
TMatrixD eth = m;
TMatrixD p(msize,msize);
for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++)
p(p.GetRowUpb()+p.GetRowLwb()-i,i) = 1;
m *= p;
m *= p;
ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test general matrix multiplication through inline mult" << std::endl;
TMatrixD m = THilbertMatrixD(msize-2,msize);
m(3,3) = TMath::Pi();
TMatrixD mt(TMatrixD::kTransposed,m);
TMatrixD p = THilbertMatrixD(msize,msize);
TMatrixDDiag(p) += 1;
TMatrixD mp(m,TMatrixD::kMult,p);
TMatrixD m1 = m;
m *= p;
ok &= VerifyMatrixIdentity(m,mp,verbose,epsilon);
TMatrixD mp1(mt,TMatrixD::kTransposeMult,p);
VerifyMatrixIdentity(m,mp1,verbose,epsilon);
ok &= ( !(m1 == m) );
TMatrixD mp2(TMatrixD::kZero,m1);
ok &= ( mp2 == 0 );
mp2.Mult(m1,p);
ok &= VerifyMatrixIdentity(m,mp2,verbose,epsilon);
if (verbose)
std::cout << "Test XP=X*P vs XP=X;XP*=P" << std::endl;
TMatrixD mp3 = m1*p;
ok &= VerifyMatrixIdentity(m,mp3,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check n * m == n * m_sym; m_sym * n == m * n; m_sym * m_sym == m * m" <<std::endl;
const TMatrixD n = THilbertMatrixD(0,msize-1,0,msize-1);
const TMatrixD m = n;
const TMatrixDSym m_sym = THilbertMatrixDSym(0,msize-1);
const TMatrixD nm1 = n * m_sym;
const TMatrixD nm2 = n * m;
const TMatrixD mn1 = m_sym * n;
const TMatrixD mn2 = m * n;
const TMatrixD mm1 = m_sym * m_sym;
const TMatrixD mm2 = m * m;
ok &= VerifyMatrixIdentity(nm1,nm2,verbose,epsilon);
ok &= VerifyMatrixIdentity(mn1,mn2,verbose,epsilon);
ok &= VerifyMatrixIdentity(mm1,mm2,verbose,epsilon);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop++;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "General Matrix Multiplications failed for size " << msize << std::endl;
break;
}
}
if (ok)
{
if (gVerbose)
std::cout << "Check to see UU' = U'U = E when U is the Haar matrix" << std::endl;
const Int_t order = 5;
const Int_t no_sub_cols = (1<<order)-5;
TMatrixD haar_sub = THaarMatrixD(5,no_sub_cols);
TMatrixD haar_sub_t(TMatrixD::kTransposed,haar_sub);
TMatrixD hsths(haar_sub_t,TMatrixD::kMult,haar_sub);
TMatrixD hsths1(TMatrixD::kZero,hsths); hsths1.Mult(haar_sub_t,haar_sub);
TMatrixD hsths_eth(TMatrixD::kUnit,hsths);
ok &= ( hsths.GetNrows() == no_sub_cols && hsths.GetNcols() == no_sub_cols );
ok &= VerifyMatrixIdentity(hsths,hsths_eth,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(hsths1,hsths_eth,gVerbose,EPSILON);
TMatrixD haar = THaarMatrixD(order);
TMatrixD unit(TMatrixD::kUnit,haar);
TMatrixD haar_t(TMatrixD::kTransposed,haar);
TMatrixD hth(haar,TMatrixD::kTransposeMult,haar);
TMatrixD hht(haar,TMatrixD::kMult,haar_t);
TMatrixD hht1 = haar; hht1 *= haar_t;
TMatrixD hht2(TMatrixD::kZero,haar); hht2.Mult(haar,haar_t);
ok &= VerifyMatrixIdentity(unit,hth,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(unit,hht,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(unit,hht1,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(unit,hht2,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(10,"General Matrix Multiplications",ok);
}
//
//------------------------------------------------------------------------
// Verify symmetric matrix multiplications
//
void mstress_sym_mm_multiplications(Int_t msize)
{
if (gVerbose)
std::cout << "\n---> Verify symmetric matrix multiplications "
"for matrices of the characteristic size " << msize << std::endl;
Int_t i,j;
Bool_t ok = kTRUE;
const Double_t epsilon = EPSILON*msize/100;
{
if (gVerbose)
std::cout << "\nTest inline multiplications of the UnitMatrix" << std::endl;
TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
TMatrixDSym u(TMatrixDSym::kUnit,m_sym);
TMatrixD u2 = u * m_sym;
ok &= VerifyMatrixIdentity(u2,m_sym,gVerbose,epsilon);
}
if (ok)
{
if (gVerbose)
std::cout << "\nTest symmetric multiplications" << std::endl;
{
if (gVerbose)
std::cout << "\n Test m * m_sym == m_sym * m == m_sym * m_sym multiplications" << std::endl;
TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
TMatrixD mm = m * m;
TMatrixD mm_sym1 = m_sym * m_sym;
TMatrixD mm_sym2 = m * m_sym;
TMatrixD mm_sym3 = m_sym * m;
ok &= VerifyMatrixIdentity(mm,mm_sym1,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mm,mm_sym2,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mm,mm_sym3,gVerbose,epsilon);
}
{
if (gVerbose)
std::cout << "\n Test m^T * m_sym == m_sym^T * m == m_sym^T * m_sym multiplications" << std::endl;
TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
TMatrixD mtm = TMatrixD(m ,TMatrixD::kTransposeMult,m);
TMatrixD mtm_sym1 = TMatrixD(m_sym,TMatrixD::kTransposeMult,m_sym);
TMatrixD mtm_sym2 = TMatrixD(m ,TMatrixD::kTransposeMult,m_sym);
TMatrixD mtm_sym3 = TMatrixD(m_sym,TMatrixD::kTransposeMult,m);
ok &= VerifyMatrixIdentity(mtm,mtm_sym1,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mtm,mtm_sym2,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mtm,mtm_sym3,gVerbose,epsilon);
}
{
if (gVerbose)
std::cout << "\n Test m * m_sym^T == m_sym * m^T == m_sym * m_sym^T multiplications" << std::endl;
TMatrixD m = THilbertMatrixD(-1,msize,-1,msize);
TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
TMatrixD mmt = TMatrixD(m ,TMatrixD::kMultTranspose,m);
TMatrixD mmt_sym1 = TMatrixD(m_sym,TMatrixD::kMultTranspose,m_sym);
TMatrixD mmt_sym2 = TMatrixD(m ,TMatrixD::kMultTranspose,m_sym);
TMatrixD mmt_sym3 = TMatrixD(m_sym,TMatrixD::kMultTranspose,m);
ok &= VerifyMatrixIdentity(mmt,mmt_sym1,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mmt,mmt_sym2,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(mmt,mmt_sym3,gVerbose,epsilon);
}
{
if (gVerbose)
std::cout << "\n Test n * m_sym == n * m multiplications" << std::endl;
TMatrixD n = THilbertMatrixD(-1,msize,-1,msize);
TMatrixD m = n;
n(1,3) = TMath::Pi();
n(3,1) = TMath::Pi();
TMatrixDSym m_sym(-1,msize,m.GetMatrixArray());
TMatrixD nm1 = n * m_sym;
TMatrixD nm2 = n * m;
ok &= VerifyMatrixIdentity(nm1,nm2,gVerbose,epsilon);
}
}
if (ok)
{
if (gVerbose)
std::cout << "Test inline multiplications by a DiagMatrix" << std::endl;
TMatrixDSym ms = THilbertMatrixDSym(msize);
ms(1,3) = TMath::Pi();
ms(3,1) = TMath::Pi();
TVectorD v(msize);
for (i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = 1+i;
TMatrixDSym diag(msize);
TMatrixDDiag d(diag); d = v;
TMatrixDSym eth = ms;
for (i = eth.GetRowLwb(); i <= eth.GetRowUpb(); i++)
for (j = eth.GetColLwb(); j <= eth.GetColUpb(); j++)
eth(i,j) *= v(j);
TMatrixD m2 = ms * diag;
ok &= VerifyMatrixIdentity(m2,eth,gVerbose,epsilon);
}
if (ok)
{
if (gVerbose)
std::cout << "Test XPP = X where P is a permutation matrix" << std::endl;
TMatrixDSym ms = THilbertMatrixDSym(msize);
ms(2,3) = TMath::Pi();
ms(3,2) = TMath::Pi();
TMatrixDSym eth = ms;
TMatrixDSym p(msize);
for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++)
p(p.GetRowUpb()+p.GetRowLwb()-i,i) = 1;
TMatrixD m2 = ms * p;
m2 *= p;
ok &= VerifyMatrixIdentity(m2,eth,gVerbose,epsilon);
}
if (ok)
{
if (gVerbose)
std::cout << "Test general matrix multiplication through inline mult" << std::endl;
TMatrixDSym ms = THilbertMatrixDSym(msize);
ms(2,3) = TMath::Pi();
ms(3,2) = TMath::Pi();
TMatrixDSym mt(TMatrixDSym::kTransposed,ms);
TMatrixDSym p = THilbertMatrixDSym(msize);
TMatrixDDiag(p) += 1;
TMatrixD mp(ms,TMatrixD::kMult,p);
TMatrixDSym m1 = ms;
TMatrixD m3(ms,TMatrixD::kMult,p);
memcpy(ms.GetMatrixArray(),m3.GetMatrixArray(),msize*msize*sizeof(Double_t));
ok &= VerifyMatrixIdentity(ms,mp,gVerbose,epsilon);
TMatrixD mp1(mt,TMatrixD::kTransposeMult,p);
ok &= VerifyMatrixIdentity(ms,mp1,gVerbose,epsilon);
ok &= ( !(m1 == ms) ) ? kTRUE : kFALSE;
TMatrixDSym mp2(TMatrixDSym::kZero,ms);
ok &= ( mp2 == 0 ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "Test XP=X*P vs XP=X;XP*=P" << std::endl;
TMatrixD mp3 = m1*p;
ok &= VerifyMatrixIdentity(ms,mp3,gVerbose,epsilon);
}
if (ok)
{
if (gVerbose)
std::cout << "Check to see UU' = U'U = E when U is the Haar matrix" << std::endl;
{
const Int_t order = 5;
const Int_t no_sub_cols = (1<<order)-5;
TMatrixD haarb = THaarMatrixD(5,no_sub_cols);
TMatrixD haarb_t(TMatrixD::kTransposed,haarb);
TMatrixD hth(haarb_t,TMatrixD::kMult,haarb);
TMatrixDSym hth1(TMatrixDSym::kAtA,haarb);
ok &= VerifyMatrixIdentity(hth,hth1,gVerbose,epsilon);
}
{
TMatrixD haar = THaarMatrixD(5);
TMatrixD unit(TMatrixD::kUnit,haar);
TMatrixD haar_t(TMatrixD::kTransposed,haar);
TMatrixDSym hths(TMatrixDSym::kAtA,haar);
TMatrixD hht(haar,TMatrixD::kMult,haar_t);
TMatrixD hht1 = haar; hht1 *= haar_t;
ok &= VerifyMatrixIdentity(unit,hths,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(unit,hht,gVerbose,epsilon);
ok &= VerifyMatrixIdentity(unit,hht1,gVerbose,epsilon);
}
}
if (ok)
{
if (gVerbose)
std::cout << "Check to see A.Similarity(Haar) = Haar * A * Haar^T" <<
" and A.SimilarityT(Haar) = Haar^T * A * Haar" << std::endl;
{
TMatrixD h = THaarMatrixD(5);
TMatrixDSym ms = THilbertMatrixDSym(1<<5);
TMatrixD hmht = h*TMatrixD(ms,TMatrixD::kMultTranspose,h);
ok &= VerifyMatrixIdentity(ms.Similarity(h),hmht,gVerbose,epsilon);
TMatrixD htmh = TMatrixD(h,TMatrixD::kTransposeMult,ms)*h;
ok &= VerifyMatrixIdentity(ms.SimilarityT(h),htmh,gVerbose,epsilon);
}
if (gVerbose)
std::cout << "Check to see A.Similarity(B_sym) = A.Similarity(B)" << std::endl;
{
TMatrixDSym nsym = THilbertMatrixDSym(5);
TMatrixD n = THilbertMatrixD(5,5);
TMatrixDSym ms = THilbertMatrixDSym(5);
ok &= VerifyMatrixIdentity(ms.Similarity(nsym),ms.Similarity(n),gVerbose,epsilon);
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(11,"Symmetric Matrix Multiplications",ok);
}
//
//------------------------------------------------------------------------
// Verify vector-matrix multiplications
//
void mstress_vm_multiplications()
{
Bool_t ok = kTRUE;
Int_t iloop = gNrLoop;
Int_t nr = 0;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Double_t epsilon = EPSILON*msize;
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
if (verbose)
std::cout << "\n---> Verify vector-matrix multiplications "
"for matrices of the characteristic size " << msize << std::endl;
{
if (verbose)
std::cout << "\nCheck shrinking a vector by multiplying by a non-sq unit matrix" << std::endl;
TVectorD vb(-2,msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()-i;
ok &= ( vb != 0 ) ? kTRUE : kFALSE;
TMatrixD mc(1,msize-2,-2,msize); // contracting matrix
mc.UnitMatrix();
TVectorD v1 = vb;
TVectorD v2 = vb;
v1 *= mc;
v2.ResizeTo(1,msize-2);
ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check expanding a vector by multiplying by a non-sq unit matrix" << std::endl;
TVectorD vb(msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()+i;
ok &= ( vb != 0 ) ? kTRUE : kFALSE;
TMatrixD me(2,msize+5,0,msize-1); // expanding matrix
me.UnitMatrix();
TVectorD v1 = vb;
TVectorD v2 = vb;
v1 *= me;
v2.ResizeTo(v1);
ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check general matrix-vector multiplication" << std::endl;
TVectorD vb(msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()+i;
TMatrixD vm(msize,1);
TMatrixDColumn(vm,0) = vb;
TMatrixD m = THilbertMatrixD(0,msize,0,msize-1);
vb *= m;
ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
TMatrixD mvm(m,TMatrixD::kMult,vm);
TMatrixD mvb(msize+1,1);
TMatrixDColumn(mvb,0) = vb;
ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check symmetric matrix-vector multiplication" << std::endl;
TVectorD vb(msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()+i;
TMatrixD vm(msize,1);
TMatrixDColumn(vm,0) = vb;
TMatrixDSym ms = THilbertMatrixDSym(0,msize-1);
vb *= ms;
ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
TMatrixD mvm(ms,TMatrixD::kMult,vm);
TMatrixD mvb(msize,1);
TMatrixDColumn(mvb,0) = vb;
ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop--;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "Vector Matrix Multiplications failed for size " << msize << std::endl;
break;
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(12,"Matrix Vector Multiplications",ok);
}
//
//------------------------------------------------------------------------
// Verify matrix inversion
//
void mstress_inversion()
{
Bool_t ok = kTRUE;
Int_t iloop = gNrLoop;
Int_t nr = 0;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Double_t epsilon = EPSILON*msize/10;
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
if (verbose)
std::cout << "\n---> Verify matrix inversion for square matrices of size " << msize << std::endl;
{
if (verbose)
std::cout << "\nTest inversion of a diagonal matrix" << std::endl;
TMatrixD m(-1,msize,-1,msize);
TMatrixD mi(TMatrixD::kZero,m);
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
m(i,i)=i-m.GetRowLwb()+1;
mi(i,i) = 1/m(i,i);
}
TMatrixD mi1(TMatrixD::kInverted,m);
m.Invert();
ok &= VerifyMatrixIdentity(m,mi,verbose,epsilon);
ok &= VerifyMatrixIdentity(mi1,mi,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test inversion of an orthonormal (Haar) matrix" << std::endl;
const Int_t order = Int_t(TMath::Log(msize)/TMath::Log(2));
TMatrixD mr = THaarMatrixD(order);
TMatrixD morig = mr;
TMatrixD mt(TMatrixD::kTransposed,mr);
Double_t det = -1; // init to a wrong val to see if it's changed
mr.Invert(&det);
ok &= VerifyMatrixIdentity(mr,mt,verbose,epsilon);
ok &= ( TMath::Abs(det-1) <= msize*epsilon ) ? kTRUE : kFALSE;
if (verbose) {
std::cout << "det = " << det << " deviation= " << TMath::Abs(det-1);
std::cout << ( (TMath::Abs(det-1) < msize*epsilon) ? " OK" : " too large") <<std::endl;
}
TMatrixD mti(TMatrixD::kInverted,mt);
ok &= VerifyMatrixIdentity(mti,morig,verbose,msize*epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test inversion of a good matrix with diagonal dominance" << std::endl;
TMatrixD mh = THilbertMatrixD(msize,msize);
TMatrixDDiag(mh) += 1;
TMatrixD morig = mh;
Double_t det_inv = 0;
const Double_t det_comp = mh.Determinant();
mh.Invert(&det_inv);
if (verbose) {
std::cout << "\tcomputed determinant " << det_comp << std::endl;
std::cout << "\tdeterminant returned by Invert() " << det_inv << std::endl;
}
if (verbose)
std::cout << "\tcheck to see M^(-1) * M is E" << std::endl;
TMatrixD mim(mh,TMatrixD::kMult,morig);
TMatrixD unit(TMatrixD::kUnit,mh);
ok &= VerifyMatrixIdentity(mim,unit,verbose,epsilon);
if (verbose)
std::cout << "Test inversion through the matrix decompositions" << std::endl;
TMatrixDSym ms = THilbertMatrixDSym(msize);
TMatrixDDiag(ms) += 1;
if (verbose)
std::cout << "Test inversion through SVD" << std::endl;
TMatrixD inv_svd (msize,msize); TDecompSVD svd (ms); svd.Invert(inv_svd);
ok &= VerifyMatrixIdentity(inv_svd,mh,verbose,epsilon);
if (verbose)
std::cout << "Test inversion through LU" << std::endl;
TMatrixD inv_lu (msize,msize); TDecompLU lu (ms); lu.Invert(inv_lu);
ok &= VerifyMatrixIdentity(inv_lu,mh,verbose,epsilon);
if (verbose)
std::cout << "Test inversion through Cholesky" << std::endl;
TMatrixDSym inv_chol(msize); TDecompChol chol(ms); chol.Invert(inv_chol);
ok &= VerifyMatrixIdentity(inv_chol,mh,verbose,epsilon);
if (verbose)
std::cout << "Test inversion through QRH" << std::endl;
TMatrixD inv_qrh (msize,msize); TDecompQRH qrh (ms); qrh.Invert(inv_qrh);
ok &= VerifyMatrixIdentity(inv_qrh,mh,verbose,epsilon);
if (verbose)
std::cout << "Test inversion through Bunch-Kaufman" << std::endl;
TMatrixDSym inv_bk(msize); TDecompBK bk(ms); chol.Invert(inv_bk);
ok &= VerifyMatrixIdentity(inv_bk,mh,verbose,epsilon);
if (verbose)
std::cout << "\tcheck to see M * M^(-1) is E" << std::endl;
TMatrixD mmi = morig; mmi *= mh;
ok &= VerifyMatrixIdentity(mmi,unit,verbose,epsilon);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop--;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "Matrix Inversion failed for size " << msize << std::endl;
break;
}
}
{
if (gVerbose) {
std::cout << "Check to see that Invert() and InvertFast() give identical results" <<std::endl;
std::cout << " for size < (7x7)" <<std::endl;
}
Int_t size;
for (size = 2; size < 7; size++) {
TMatrixD m1 = THilbertMatrixD(size,size);
m1(0,1) = TMath::Pi();
TMatrixDDiag(m1) += 1;
TMatrixD m2 = m1;
Double_t det1 = 0.0;
Double_t det2 = 0.0;
m1.Invert(&det1);
m2.InvertFast(&det2);
ok &= VerifyMatrixIdentity(m1,m2,gVerbose,EPSILON);
ok &= (TMath::Abs(det1-det2) < EPSILON);
if (gVerbose) {
std::cout << "det(Invert)= " << det1 << " det(InvertFast)= " << det2 <<std::endl;
std::cout << " deviation= " << TMath::Abs(det1-det2);
std::cout << ( (TMath::Abs(det1-det2) < EPSILON) ? " OK" : " too large") <<std::endl;
}
}
for (size = 2; size < 7; size++) {
TMatrixDSym ms1 = THilbertMatrixDSym(size);
TMatrixDDiag(ms1) += 1;
TMatrixDSym ms2 = ms1;
Double_t det1 = 0.0;
Double_t det2 = 0.0;
ms1.Invert(&det1);
ms2.InvertFast(&det2);
ok &= VerifyMatrixIdentity(ms1,ms2,gVerbose,EPSILON);
ok &= (TMath::Abs(det1-det2) < EPSILON);
if (gVerbose) {
std::cout << "det(Invert)= " << det1 << " det(InvertFast)= " << det2 <<std::endl;
std::cout << " deviation= " << TMath::Abs(det1-det2);
std::cout << ( (TMath::Abs(det1-det2) < EPSILON) ? " OK" : " too large") <<std::endl;
}
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(13,"Matrix Inversion",ok);
}
//
//------------------------------------------------------------------------
// Test matrix I/O
//
void mstress_matrix_io()
{
if (gVerbose)
std::cout << "\n---> Test matrix I/O" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = TMath::Pi();
TFile *f = new TFile("stress-vmatrix.root", "RECREATE");
Char_t name[80];
Int_t iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TMatrixD m(msize,msize);
m = pattern;
Double_t *pattern_array = new Double_t[msize*msize];
for (Int_t i = 0; i < msize*msize; i++)
pattern_array[i] = pattern;
TMatrixD ma;
ma.Use(msize,msize,pattern_array);
TMatrixDSym ms(msize);
ms = pattern;
if (verbose)
std::cout << "\nWrite matrix m to database" << std::endl;
snprintf(name,80,"m_%d",msize);
m.Write(name);
if (verbose)
std::cout << "\nWrite matrix ma which adopts to database" << std::endl;
snprintf(name,80,"ma_%d",msize);
ma.Write(name);
if (verbose)
std::cout << "\nWrite symmetric matrix ms to database" << std::endl;
snprintf(name,80,"ms_%d",msize);
ms.Write(name);
delete [] pattern_array;
iloop--;
}
if (gVerbose)
std::cout << "\nClose database" << std::endl;
delete f;
if (gVerbose)
std::cout << "\nOpen database in read-only mode and read matrix" << std::endl;
TFile *f1 = new TFile("stress-vmatrix.root");
iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TMatrixD m(msize,msize);
m = pattern;
snprintf(name,80,"m_%d",msize);
TMatrixD *mr = (TMatrixD*) f1->Get(name);
snprintf(name,80,"ma_%d",msize);
TMatrixD *mar = (TMatrixD*) f1->Get(name);
snprintf(name,80,"ms_%d",msize);
TMatrixDSym *msr = (TMatrixDSym*) f1->Get(name);
if (verbose)
std::cout << "\nRead matrix should be same as original" << std::endl;
ok &= ((*mr) == m) ? kTRUE : kFALSE;
ok &= ((*mar) == m) ? kTRUE : kFALSE;
ok &= ((*msr) == m) ? kTRUE : kFALSE;
iloop--;
}
delete f1;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(14,"Matrix Persistence",ok);
}
//------------------------------------------------------------------------
// Test allocation functions and compatibility check
//
void spstress_allocation(Int_t msize)
{
if (gVerbose)
std::cout << "\n\n---> Test allocation and compatibility check" << std::endl;
Int_t i,j;
Bool_t ok = kTRUE;
TMatrixDSparse m1(0,3,0,msize-1);
{
Int_t nr = 4*msize;
Int_t *irow = new Int_t[nr];
Int_t *icol = new Int_t[nr];
Double_t *val = new Double_t[nr];
Int_t n = 0;
for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
irow[n] = i;
icol[n] = j;
val[n] = TMath::Pi()*i+TMath::E()*j;
n++;
}
}
m1.SetMatrixArray(nr,irow,icol,val);
delete [] irow;
delete [] icol;
delete [] val;
}
TMatrixDSparse m2(0,3,0,msize-1);
TMatrixDSparse m3(1,4,0,msize-1);
TMatrixDSparse m4(m1);
if (gVerbose) {
std::cout << "\nStatus information reported for matrix m3:" << std::endl;
std::cout << " Row lower bound ... " << m3.GetRowLwb() << std::endl;
std::cout << " Row upper bound ... " << m3.GetRowUpb() << std::endl;
std::cout << " Col lower bound ... " << m3.GetColLwb() << std::endl;
std::cout << " Col upper bound ... " << m3.GetColUpb() << std::endl;
std::cout << " No. rows ..........." << m3.GetNrows() << std::endl;
std::cout << " No. cols ..........." << m3.GetNcols() << std::endl;
std::cout << " No. of elements ...." << m3.GetNoElements() << std::endl;
}
if (gVerbose)
std::cout << "Check matrices 1 & 4 for compatibility" << std::endl;
ok &= AreCompatible(m1,m4,gVerbose);
if (gVerbose)
std::cout << "m2 has to be compatible with m3 after resizing to m3" << std::endl;
m2.ResizeTo(m3);
ok &= AreCompatible(m2,m3,gVerbose);
TMatrixD m5_d(m1.GetNrows()+1,m1.GetNcols()+5);
for (i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
for (j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
m5_d(i,j) = TMath::Pi()*i+TMath::E()*j;
TMatrixDSparse m5(m5_d);
if (gVerbose)
std::cout << "m1 has to be compatible with m5 after resizing to m5" << std::endl;
m1.ResizeTo(m5.GetNrows(),m5.GetNcols());
ok &= AreCompatible(m1,m5,gVerbose);
if (gVerbose)
std::cout << "m1 has to be equal to m4 after stretching and shrinking" << std::endl;
m1.ResizeTo(m4.GetNrows(),m4.GetNcols());
ok &= VerifyMatrixIdentity(m1,m4,gVerbose,EPSILON);
if (gVerbose)
std::cout << "m5 has to be equal to m1 after shrinking" << std::endl;
m5.ResizeTo(m1.GetNrows(),m1.GetNcols());
ok &= VerifyMatrixIdentity(m1,m5,gVerbose,msize*EPSILON);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(1,"Allocation, Resizing",ok);
}
//
//------------------------------------------------------------------------
// Test Filling of matrix
//
void spstress_matrix_fill(Int_t rsize,Int_t csize)
{
Bool_t ok = kTRUE;
if (csize < 4) {
Error("spstress_matrix_fill","rsize should be >= 4");
ok = kFALSE;
StatusPrint(2,"Filling, Inserting, Using",ok);
return;
}
if (csize < 4) {
Error("spstress_matrix_fill","csize should be >= 4");
ok = kFALSE;
StatusPrint(2,"Filling, Inserting, Using",ok);
return;
}
TMatrixD m_d(-1,rsize-2,1,csize);
for (Int_t i = m_d.GetRowLwb(); i <= m_d.GetRowUpb(); i++)
for (Int_t j = m_d.GetColLwb(); j <= m_d.GetColUpb(); j++)
// Keep column 2 on zero
if (j != 2)
m_d(i,j) = TMath::Pi()*i+TMath::E()*j;
TMatrixDSparse m(m_d);
{
if (gVerbose)
std::cout << "Check filling through operator(i,j) without setting sparse index" << std::endl;
TMatrixDSparse m1(-1,rsize-2,1,csize);
for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++)
for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++)
if (j != 2)
m1(i,j) = TMath::Pi()*i+TMath::E()*j;
ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
}
{
if (gVerbose)
std::cout << "Check filling through operator(i,j)" << std::endl;
TMatrixDSparse m2(-1,rsize-2,1,csize);
m2.SetSparseIndex(m);
for (Int_t i = m2.GetRowLwb(); i <= m2.GetRowUpb(); i++)
for (Int_t j = m2.GetColLwb(); j <= m2.GetColUpb(); j++)
if (j != 2)
m2(i,j) = TMath::Pi()*i+TMath::E()*j;
ok &= VerifyMatrixIdentity(m2,m,gVerbose,EPSILON);
}
{
if (gVerbose)
std::cout << "Check insertion/extraction of sub-matrices" << std::endl;
{
TMatrixDSparse m_sub1 = m;
m_sub1.ResizeTo(0,rsize-2,2,csize);
TMatrixDSparse m_sub2 = m.GetSub(0,rsize-2,2,csize,"");
ok &= VerifyMatrixIdentity(m_sub1,m_sub2,gVerbose,EPSILON);
}
{
TMatrixDSparse m3(-1,rsize-2,1,csize);
TMatrixDSparse m_part1 = m.GetSub(-1,rsize-2,1,csize,"");
m3.SetSub(-1,1,m_part1);
ok &= VerifyMatrixIdentity(m,m3,gVerbose,EPSILON);
}
{
TMatrixDSparse m4(-1,rsize-2,1,csize);
TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"");
TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"");
TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"");
TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"");
m4.SetSub(0,2,m_part1);
m4.SetSub(0,1,m_part2);
m4.SetSub(-1,2,m_part3);
m4.SetSub(-1,1,m_part4);
ok &= VerifyMatrixIdentity(m,m4,gVerbose,EPSILON);
}
{
// change the insertion order
TMatrixDSparse m5(-1,rsize-2,1,csize);
TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"");
TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"");
TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"");
TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"");
m5.SetSub(-1,1,m_part4);
m5.SetSub(-1,2,m_part3);
m5.SetSub(0,1,m_part2);
m5.SetSub(0,2,m_part1);
ok &= VerifyMatrixIdentity(m,m5,gVerbose,EPSILON);
}
{
TMatrixDSparse m6(-1,rsize-2,1,csize);
TMatrixDSparse m_part1 = m.GetSub(0,rsize-2,2,csize,"S");
TMatrixDSparse m_part2 = m.GetSub(0,rsize-2,1,1,"S");
TMatrixDSparse m_part3 = m.GetSub(-1,-1,2,csize,"S");
TMatrixDSparse m_part4 = m.GetSub(-1,-1,1,1,"S");
m6.SetSub(0,2,m_part1);
m6.SetSub(0,1,m_part2);
m6.SetSub(-1,2,m_part3);
m6.SetSub(-1,1,m_part4);
ok &= VerifyMatrixIdentity(m,m6,gVerbose,EPSILON);
}
}
{
if (gVerbose)
std::cout << "Check insertion/extraction of rows" << std::endl;
TMatrixDSparse m1 = m;
TVectorD v1(1,csize);
TVectorD v2(1,csize);
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
v1 = TMatrixDSparseRow(m,i);
m1.InsertRow(i,1,v1.GetMatrixArray());
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
m1.InsertRow(i,3,v1.GetMatrixArray()+2,v1.GetNrows()-2);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
m1.ExtractRow(i,1,v2.GetMatrixArray());
ok &= VerifyVectorIdentity(v1,v2,gVerbose,EPSILON);
m1.ExtractRow(i,3,v2.GetMatrixArray()+2,v1.GetNrows()-2);
ok &= VerifyVectorIdentity(v1,v2,gVerbose,EPSILON);
}
}
{
if (gVerbose)
std::cout << "Check array Use" << std::endl;
{
TMatrixDSparse *m1a = new TMatrixDSparse(m);
TMatrixDSparse *m2a = new TMatrixDSparse();
m2a->Use(*m1a);
m2a->Sqr();
TMatrixDSparse m7 = m; m7.Sqr();
ok &= VerifyMatrixIdentity(m7,*m1a,gVerbose,EPSILON);
delete m1a;
delete m2a;
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(2,"Filling, Inserting, Using",ok);
}
//
//------------------------------------------------------------------------
// Test uniform element operations
//
void spstress_element_op(Int_t rsize,Int_t csize)
{
Bool_t ok = kTRUE;
const Double_t pattern = 8.625;
TMatrixDSparse m(-1,rsize-2,1,csize);
if (gVerbose)
std::cout << "Creating zero m1 ..." << std::endl;
TMatrixDSparse m1(TMatrixDSparse::kZero, m);
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing m1 with 0 ..." << std::endl;
R__ASSERT(m1 == 0);
R__ASSERT(!(m1 != 0));
if (gVerbose)
std::cout << "Writing a pattern " << pattern << " by assigning through SetMatrixArray..." << std::endl;
{
const Int_t nr = rsize*csize;
Int_t *irow = new Int_t[nr];
Int_t *icol = new Int_t[nr];
Double_t *val = new Double_t[nr];
Int_t n = 0;
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
irow[n] = i;
icol[n] = j;
val[n] = pattern;
n++;
}
}
m.SetMatrixArray(nr,irow,icol,val);
delete [] irow;
delete [] icol;
delete [] val;
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "Writing the pattern by assigning to m1 as a whole ..." << std::endl;
m1.SetSparseIndex(m);
m1 = pattern;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing m and m1 ..." << std::endl;
R__ASSERT(m == m1);
if (gVerbose)
std::cout << "Comparing (m=0) and m1 ..." << std::endl;
R__ASSERT(!((m=0) == m1));
if (gVerbose)
std::cout << "Clearing m1 ..." << std::endl;
m1.Zero();
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nSet m = pattern" << std::endl;
m = pattern;
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " add the doubled pattern with the negative sign" << std::endl;
m += -2*pattern;
ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtract the trippled pattern with the negative sign" << std::endl;
m -= -3*pattern;
ok &= VerifyMatrixValue(m,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify comparison operations when all elems are the same" << std::endl;
m = pattern;
R__ASSERT( m == pattern && !(m != pattern) );
R__ASSERT( m > 0 && m >= pattern && m <= pattern );
R__ASSERT( m > -pattern && m >= -pattern );
R__ASSERT( m <= pattern && !(m < pattern) );
m -= 2*pattern;
R__ASSERT( m < -pattern/2 && m <= -pattern/2 );
R__ASSERT( m >= -pattern && !(m > -pattern) );
if (gVerbose)
std::cout << "\nVerify comparison operations when not all elems are the same" << std::endl;
{
Int_t nr = rsize*csize;
Int_t *irow = new Int_t[nr];
Int_t *icol = new Int_t[nr];
Double_t *val = new Double_t[nr];
Int_t n = 0;
for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
irow[n] = i;
icol[n] = j;
val[n] = pattern;
n++;
}
}
val[n-1] = pattern-1;
m.SetMatrixArray(nr,irow,icol,val);
delete [] irow;
delete [] icol;
delete [] val;
}
R__ASSERT( !(m == pattern) && !(m != pattern) );
R__ASSERT( m != 0 ); // none of elements are 0
R__ASSERT( !(m >= pattern) && m <= pattern && !(m<pattern) );
R__ASSERT( !(m <= pattern-1) && m >= pattern-1 && !(m>pattern-1) );
if (gVerbose)
std::cout << "\nAssign 2*pattern to m by repeating additions" << std::endl;
m = 0; m += pattern; m += pattern;
if (gVerbose)
std::cout << "Assign 2*pattern to m1 by multiplying by two " << std::endl;
m1.SetSparseIndex(m);
m1 = pattern; m1 *= 2;
ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
R__ASSERT( m == m1 );
if (gVerbose)
std::cout << "Multiply m1 by one half returning it to the 1*pattern" << std::endl;
m1 *= 1/2.;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nAssign -pattern to m and m1" << std::endl;
m = 0; m -= pattern; m1 = -pattern;
ok &= VerifyMatrixValue(m,-pattern,gVerbose,EPSILON);
R__ASSERT( m == m1 );
if (gVerbose)
std::cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same" << std::endl;
m.Sqr();
ok &= VerifyMatrixValue(m,pattern*pattern,gVerbose,EPSILON);
m.Sqrt();
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
m1.Abs();
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
ok &= VerifyMatrixIdentity(m1,m,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(3,"Uniform matrix operations",ok);
}
//
//------------------------------------------------------------------------
// Test binary matrix element-by-element operations
//
void spstress_binary_ebe_op(Int_t rsize, Int_t csize)
{
if (gVerbose)
std::cout << "\n---> Test Binary Matrix element-by-element operations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 4.25;
TMatrixD m_d(2,rsize+1,0,csize-1); m_d = 1;
TMatrixDSparse m (TMatrixDSparse::kZero,m_d); m.SetSparseIndex (m_d);
TMatrixDSparse m1(TMatrixDSparse::kZero,m); m1.SetSparseIndex(m_d);
TMatrixDSparse mp(TMatrixDSparse::kZero,m1); mp.SetSparseIndex(m_d);
{
for (Int_t i = mp.GetRowLwb(); i <= mp.GetRowUpb(); i++)
for (Int_t j = mp.GetColLwb(); j <= mp.GetColUpb(); j++)
mp(i,j) = TMath::Pi()*i+TMath::E()*(j+1);
}
if (gVerbose)
std::cout << "\nVerify assignment of a matrix to the matrix" << std::endl;
m = pattern;
m1 = m;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
R__ASSERT( m1 == m );
if (gVerbose)
std::cout << "\nAdding the matrix to itself, uniform pattern " << pattern << std::endl;
m.Zero(); m.SetSparseIndex(m_d); m = pattern;
m1 = m; m1 += m1;
ok &= VerifyMatrixValue(m1,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting two matrices ..." << std::endl;
m1 -= m;
ok &= VerifyMatrixValue(m1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting the matrix from itself" << std::endl;
m1 -= m1;
ok &= VerifyMatrixValue(m1,0.,gVerbose,EPSILON);
m1.SetSparseIndex(m_d);
if (gVerbose) {
std::cout << "\nArithmetic operations on matrices with not the same elements" << std::endl;
std::cout << " adding mp to the zero matrix..." << std::endl;
}
m.Zero(); m.SetSparseIndex(m_d); m += mp;
ok &= VerifyMatrixIdentity(m,mp,gVerbose,EPSILON);
if (gVerbose)
std::cout << " making m = 3*mp and m1 = 3*mp, via add() and succesive mult" << std::endl;
m1 = m;
Add(m,2.,mp);
m1 += m1; m1 += mp;
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose)
std::cout << " clear both m and m1, by subtracting from itself and via add()" << std::endl;
m1 -= m1;
Add(m,-3.,mp);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose) {
std::cout << "\nTesting element-by-element multiplications and divisions" << std::endl;
std::cout << " squaring each element with sqr() and via multiplication" << std::endl;
}
m.SetSparseIndex(m_d); m = mp;
m1.SetSparseIndex(m_d); m1 = mp;
m.Sqr();
ElementMult(m1,m1);
ok &= VerifyMatrixIdentity(m,m1,gVerbose,EPSILON);
if (gVerbose)
std::cout << " compare (m = pattern^2)/pattern with pattern" << std::endl;
m = pattern; m1 = pattern;
m.Sqr();
ElementDiv(m,m1);
ok &= VerifyMatrixValue(m,pattern,gVerbose,EPSILON);
if (gVerbose)
Compare(m1,m);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(4,"Binary Matrix element-by-element operations",ok);
}
//
//------------------------------------------------------------------------
// Verify matrix transposition
//
void spstress_transposition(Int_t msize)
{
if (gVerbose) {
std::cout << "\n---> Verify matrix transpose "
"for matrices of a characteristic size " << msize << std::endl;
}
Bool_t ok = kTRUE;
{
if (gVerbose)
std::cout << "\nCheck to see that a square UnitMatrix stays the same";
TMatrixDSparse m(msize,msize);
m.UnitMatrix();
TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
ok &= ( m == mt ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nTest a non-square UnitMatrix";
TMatrixDSparse m(msize,msize+1);
m.UnitMatrix();
TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows() );
const Int_t rowlwb = m.GetRowLwb();
const Int_t collwb = m.GetColLwb();
const Int_t upb = TMath::Min(m.GetRowUpb(),m.GetColUpb());
TMatrixDSparse m_part = m.GetSub(rowlwb,upb,collwb,upb);
TMatrixDSparse mt_part = mt.GetSub(rowlwb,upb,collwb,upb);
ok &= ( m_part == mt_part ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nCheck to see that a symmetric (Hilbert)Matrix stays the same";
TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize,msize));
TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
ok &= ( m == mt ) ? kTRUE : kFALSE;
}
{
if (gVerbose)
std::cout << "\nCheck transposing a non-symmetric matrix";
TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize+1,msize));
m(1,2) = TMath::Pi();
TMatrixDSparse mt(TMatrixDSparse::kTransposed,m);
R__ASSERT(m.GetNrows() == mt.GetNcols() && m.GetNcols() == mt.GetNrows());
R__ASSERT(mt(2,1) == (Double_t)TMath::Pi() && mt(1,2) != (Double_t)TMath::Pi());
R__ASSERT(mt[2][1] == (Double_t)TMath::Pi() && mt[1][2] != (Double_t)TMath::Pi());
if (gVerbose)
std::cout << "\nCheck double transposing a non-symmetric matrix" << std::endl;
TMatrixDSparse mtt(TMatrixDSparse::kTransposed,mt);
ok &= ( m == mtt ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(5,"Matrix transposition",ok);
}
//
//------------------------------------------------------------------------
// Verify the norm calculation
//
void spstress_norms(Int_t rsize_req,Int_t csize_req)
{
if (gVerbose)
std::cout << "\n---> Verify norm calculations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 10.25;
Int_t rsize = rsize_req;
if (rsize%2 != 0) rsize--;
Int_t csize = csize_req;
if (csize%2 != 0) csize--;
if (rsize%2 == 1 || csize%2 == 1) {
std::cout << "rsize: " << rsize <<std::endl;
std::cout << "csize: " << csize <<std::endl;
Fatal("spstress_norms","Sorry, size of the matrix to test must be even for this test\n");
}
TMatrixD m_d(rsize,csize); m_d = 1;
TMatrixDSparse m(rsize,csize); m.SetSparseIndex(m_d);
if (gVerbose)
std::cout << "\nAssign " << pattern << " to all the elements and check norms" << std::endl;
m = pattern;
if (gVerbose)
std::cout << " 1. (col) norm should be pattern*nrows" << std::endl;
ok &= ( m.Norm1() == pattern*m.GetNrows() ) ? kTRUE : kFALSE;
ok &= ( m.Norm1() == m.ColNorm() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Inf (row) norm should be pattern*ncols" << std::endl;
ok &= ( m.NormInf() == pattern*m.GetNcols() ) ? kTRUE : kFALSE;
ok &= ( m.NormInf() == m.RowNorm() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Square of the Eucl norm has got to be pattern^2 * no_elems" << std::endl;
ok &= ( m.E2Norm() == (pattern*pattern)*m.GetNoElements() ) ? kTRUE : kFALSE;
TMatrixDSparse m1(m); m1 = 1;
ok &= ( m.E2Norm() == E2Norm(m+1.,m1) ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(6,"Matrix Norms",ok);
}
//
//------------------------------------------------------------------------
// Verify matrix multiplications
//
void spstress_mm_multiplications()
{
Bool_t ok = kTRUE;
Int_t i;
Int_t iloop = 0;
Int_t nr = 0;
while (iloop <= gNrLoop) {
const Int_t msize = gSizeA[iloop];
const Double_t epsilon = EPSILON*msize/100;
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
if (msize <= 5) {
iloop++;
continue; // some references to m(3,..)
}
if (verbose)
std::cout << "\n---> Verify matrix multiplications "
"for matrices of the characteristic size " << msize << std::endl;
{
if (verbose)
std::cout << "\nTest inline multiplications of the UnitMatrix" << std::endl;
if (ok)
{
TMatrixDSparse m = TMatrixD(THilbertMatrixD(-1,msize,-1,msize));
TMatrixDSparse ur(TMatrixDSparse::kUnit,m);
TMatrixDSparse ul(TMatrixDSparse::kUnit,m);
m(3,1) = TMath::Pi();
ul *= m;
m *= ur;
ok &= VerifyMatrixIdentity(ul,m,verbose,epsilon);
}
if (ok)
{
TMatrixD m_d = THilbertMatrixD(-1,msize,-1,msize);
TMatrixDSparse ur(TMatrixDSparse::kUnit,m_d);
TMatrixDSparse ul(TMatrixDSparse::kUnit,m_d);
m_d(3,1) = TMath::Pi();
ul *= m_d;
m_d *= TMatrixD(ur);
ok &= VerifyMatrixIdentity(ul,m_d,verbose,epsilon);
}
}
if (ok)
{
if (verbose)
std::cout << "Test XPP = X where P is a permutation matrix" << std::endl;
TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize-1,msize));
m(2,3) = TMath::Pi();
TMatrixDSparse eth = m;
TMatrixDSparse p(msize,msize);
{
Int_t *irow = new Int_t[msize];
Int_t *icol = new Int_t[msize];
Double_t *val = new Double_t[msize];
Int_t n = 0;
for (i = p.GetRowLwb(); i <= p.GetRowUpb(); i++) {
irow[n] = p.GetRowUpb()+p.GetRowLwb()-i;
icol[n] = i;
val[n] = 1;
n++;
}
p.SetMatrixArray(msize,irow,icol,val);
delete [] irow;
delete [] icol;
delete [] val;
}
m *= p;
m *= p;
ok &= VerifyMatrixIdentity(m,eth,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Test general matrix multiplication through inline mult" << std::endl;
TMatrixDSparse m = TMatrixD(THilbertMatrixD(msize-2,msize));
m(3,3) = TMath::Pi();
TMatrixD p_d = THilbertMatrixD(msize,msize);
TMatrixDDiag(p_d) += 1;
TMatrixDSparse mp(m,TMatrixDSparse::kMult,p_d);
TMatrixDSparse m1 = m;
m *= p_d;
ok &= VerifyMatrixIdentity(m,mp,verbose,epsilon);
TMatrixDSparse pt_d(TMatrixDSparse::kTransposed,p_d);
TMatrixDSparse mp1(m1,TMatrixDSparse::kMultTranspose,pt_d);
VerifyMatrixIdentity(m,mp1,verbose,epsilon);
ok &= ( !(m1 == m) );
TMatrixDSparse mp2(TMatrixDSparse::kZero,m1);
ok &= ( mp2 == 0 );
mp2.SetSparseIndex(m1);
mp2.Mult(m1,p_d);
ok &= VerifyMatrixIdentity(m,mp2,verbose,epsilon);
if (verbose)
std::cout << "Test XP=X*P vs XP=X;XP*=P" << std::endl;
TMatrixDSparse mp3 = m1*p_d;
ok &= VerifyMatrixIdentity(m,mp3,verbose,epsilon);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop++;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "General Matrix Multiplications failed for size " << msize << std::endl;
break;
}
}
if (ok)
{
if (gVerbose)
std::cout << "Check to see UU' = U'U = E when U is the Haar matrix" << std::endl;
const Int_t order = 5;
const Int_t no_sub_cols = (1<<order)-order;
TMatrixDSparse haar_sub = TMatrixD(THaarMatrixD(order,no_sub_cols));
TMatrixDSparse haar_sub_t(TMatrixDSparse::kTransposed,haar_sub);
TMatrixDSparse hsths(haar_sub_t,TMatrixDSparse::kMult,haar_sub);
TMatrixDSparse square(no_sub_cols,no_sub_cols);
TMatrixDSparse units(TMatrixDSparse::kUnit,square);
ok &= ( hsths.GetNrows() == no_sub_cols && hsths.GetNcols() == no_sub_cols );
ok &= VerifyMatrixIdentity(hsths,units,gVerbose,EPSILON);
TMatrixDSparse haar = TMatrixD(THaarMatrixD(order));
TMatrixDSparse haar_t(TMatrixDSparse::kTransposed,haar);
TMatrixDSparse unit(TMatrixDSparse::kUnit,haar);
TMatrixDSparse hth(haar_t,TMatrixDSparse::kMult,haar);
ok &= VerifyMatrixIdentity(hth,unit,gVerbose,EPSILON);
TMatrixDSparse hht(haar,TMatrixDSparse::kMultTranspose,haar);
ok &= VerifyMatrixIdentity(hht,unit,gVerbose,EPSILON);
TMatrixDSparse hht1 = haar; hht1 *= haar_t;
ok &= VerifyMatrixIdentity(hht1,unit,gVerbose,EPSILON);
TMatrixDSparse hht2(TMatrixDSparse::kZero,haar);
hht2.SetSparseIndex(hht1);
hht2.Mult(haar,haar_t);
ok &= VerifyMatrixIdentity(hht2,unit,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(7,"General Matrix Multiplications",ok);
}
//
//------------------------------------------------------------------------
// Verify vector-matrix multiplications
//
void spstress_vm_multiplications()
{
Bool_t ok = kTRUE;
Int_t iloop = gNrLoop;
Int_t nr = 0;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Double_t epsilon = EPSILON*msize;
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
if (verbose)
std::cout << "\n---> Verify vector-matrix multiplications "
"for matrices of the characteristic size " << msize << std::endl;
{
if (verbose)
std::cout << "\nCheck shrinking a vector by multiplying by a non-sq unit matrix" << std::endl;
TVectorD vb(-2,msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()-i;
ok &= ( vb != 0 ) ? kTRUE : kFALSE;
TMatrixDSparse mc(1,msize-2,-2,msize); // contracting matrix
mc.UnitMatrix();
TVectorD v1 = vb;
TVectorD v2 = vb;
v1 *= mc;
v2.ResizeTo(1,msize-2);
ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check expanding a vector by multiplying by a non-sq unit matrix" << std::endl;
TVectorD vb(msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()+i;
ok &= ( vb != 0 ) ? kTRUE : kFALSE;
TMatrixDSparse me(2,msize+5,0,msize-1); // expanding matrix
me.UnitMatrix();
TVectorD v1 = vb;
TVectorD v2 = vb;
v1 *= me;
v2.ResizeTo(v1);
ok &= VerifyVectorIdentity(v1,v2,verbose,epsilon);
}
if (ok)
{
if (verbose)
std::cout << "Check general matrix-vector multiplication" << std::endl;
TVectorD vb(msize);
for (Int_t i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = TMath::Pi()+i;
TMatrixD vm(msize,1);
TMatrixDColumn(vm,0) = vb;
TMatrixD hilbert_with_zeros = THilbertMatrixD(0,msize,0,msize-1);
TMatrixDRow (hilbert_with_zeros,3) = 0.0;
TMatrixDColumn(hilbert_with_zeros,3) = 0.0;
const TMatrixDSparse m = hilbert_with_zeros;
vb *= m;
ok &= ( vb.GetLwb() == 0 ) ? kTRUE : kFALSE;
TMatrixDSparse mvm(m,TMatrixDSparse::kMult,vm);
TMatrixD mvb(msize+1,1);
TMatrixDColumn(mvb,0) = vb;
ok &= VerifyMatrixIdentity(mvb,mvm,verbose,epsilon);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop--;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "Vector Matrix Multiplications failed for size " << msize << std::endl;
break;
}
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(8,"Matrix Vector Multiplications",ok);
}
//
//------------------------------------------------------------------------
// Test operations with vectors and sparse matrix slices
//
void spstress_matrix_slices(Int_t vsize)
{
if (gVerbose)
std::cout << "\n---> Test operations with vectors and sparse matrix slices" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 8.625;
TVectorD vc(0,vsize);
TVectorD vr(0,vsize+1);
TMatrixD m_d(0,vsize,0,vsize+1); m_d = pattern;
TMatrixDSparse m(m_d);
Int_t i,j;
if (gVerbose)
std::cout << "\nCheck modifying the matrix row-by-row" << std::endl;
m = pattern;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
TMatrixDSparseRow(m,i) = pattern+2;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
vr = TMatrixDSparseRow(m,i);
ok &= VerifyVectorValue(vr,pattern+2,gVerbose,EPSILON);
vr = TMatrixDSparseRow(m,i+1 > m.GetRowUpb() ? m.GetRowLwb() : i+1);
ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
TMatrixDSparseRow(m,i) += -2;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
vr = pattern-2;
TMatrixDSparseRow(m,i) = vr;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
{
TMatrixDSparseRow mr(m,i);
for (j = m.GetColLwb(); j <= m.GetColUpb(); j++)
mr(j) *= 8;
}
vr = TMatrixDSparseRow(m,i);
ok &= VerifyVectorValue(vr,8*(pattern-2),gVerbose,EPSILON);
TMatrixDSparseRow(m,i) *= 1./8;
TMatrixDSparseRow(m,i) += 2;
vr = TMatrixDSparseRow(m,i);
ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nCheck modifying the matrix diagonal" << std::endl;
m = pattern;
TMatrixDSparseDiag td = m;
td = pattern-3;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
vc = TMatrixDSparseDiag(m);
ok &= VerifyVectorValue(vc,pattern-3,gVerbose,EPSILON);
td += 3;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
vc = pattern+3;
td = vc;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
{
TMatrixDSparseDiag md(m);
for (j = 0; j < md.GetNdiags(); j++)
md(j) /= 1.5;
}
vc = TMatrixDSparseDiag(m);
ok &= VerifyVectorValue(vc,(pattern+3)/1.5,gVerbose,EPSILON);
TMatrixDSparseDiag(m) *= 1.5;
TMatrixDSparseDiag(m) += -3;
vc = TMatrixDSparseDiag(m);
ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
ok &= ( m == pattern ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(9,"Matrix Slices to Vectors",ok);
}
//
//------------------------------------------------------------------------
// Test matrix I/O
//
void spstress_matrix_io()
{
if (gVerbose)
std::cout << "\n---> Test matrix I/O" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = TMath::Pi();
TFile *f = new TFile("stress-vmatrix.root", "RECREATE");
Char_t name[80];
Int_t iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TMatrixD m_d(msize,msize); m_d = pattern;
TMatrixDSparse m(m_d);
TMatrixDSparse ma;
ma.Use(m);
if (verbose)
std::cout << "\nWrite matrix m to database" << std::endl;
snprintf(name,80,"m_%d",msize);
m.Write(name);
if (verbose)
std::cout << "\nWrite matrix ma which adopts to database" << std::endl;
snprintf(name,80,"ma_%d",msize);
ma.Write(name);
iloop--;
}
if (gVerbose)
std::cout << "\nClose database" << std::endl;
delete f;
if (gVerbose)
std::cout << "\nOpen database in read-only mode and read matrix" << std::endl;
TFile *f1 = new TFile("stress-vmatrix.root");
iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TMatrixD m_d(msize,msize); m_d = pattern;
TMatrixDSparse m(m_d);
snprintf(name,80,"m_%d",msize);
TMatrixDSparse *mr = (TMatrixDSparse*) f1->Get(name);
snprintf(name,80,"ma_%d",msize);
TMatrixDSparse *mar = (TMatrixDSparse*) f1->Get(name);
snprintf(name,80,"ms_%d",msize);
if (verbose)
std::cout << "\nRead matrix should be same as original" << std::endl;
ok &= ((*mr) == m) ? kTRUE : kFALSE;
ok &= ((*mar) == m) ? kTRUE : kFALSE;
iloop--;
}
delete f1;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(10,"Matrix Persistence",ok);
}
//------------------------------------------------------------------------
// Test allocation functions and compatibility check
//
void vstress_allocation(Int_t msize)
{
if (gVerbose)
std::cout << "\n\n---> Test allocation and compatibility check" << std::endl;
Int_t i;
Bool_t ok = kTRUE;
TVectorD v1(msize);
TVectorD v2(0,msize-1);
TVectorD v3(1,msize);
TVectorD v4(v1);
if (gVerbose) {
std::cout << "\nStatus information reported for vector v3:" << std::endl;
std::cout << " Lower bound ... " << v3.GetLwb() << std::endl;
std::cout << " Upper bound ... " << v3.GetUpb() << std::endl;
std::cout << " No. of elements " << v3.GetNoElements() << std::endl;
}
if (gVerbose)
std::cout << "\nCheck vectors 1 & 2 for compatibility" << std::endl;
ok &= AreCompatible(v1,v2,gVerbose);
if (gVerbose)
std::cout << "Check vectors 1 & 4 for compatibility" << std::endl;
ok &= AreCompatible(v1,v4,gVerbose);
if (gVerbose)
std::cout << "v2 has to be compatible with v3 after resizing to v3" << std::endl;
v2.ResizeTo(v3);
ok &= AreCompatible(v2,v3,gVerbose);
TVectorD v5(v1.GetUpb()+5);
if (gVerbose)
std::cout << "v1 has to be compatible with v5 after resizing to v5.upb" << std::endl;
v1.ResizeTo(v5.GetNoElements());
ok &= AreCompatible(v1,v5,gVerbose);
{
if (gVerbose)
std::cout << "Check that shrinking does not change remaining elements" << std::endl;
TVectorD vb(-1,msize);
for (i = vb.GetLwb(); i <= vb.GetUpb(); i++)
vb(i) = i+TMath::Pi();
TVectorD v = vb;
ok &= ( v == vb ) ? kTRUE : kFALSE;
ok &= ( v != 0 ) ? kTRUE : kFALSE;
v.ResizeTo(0,msize/2);
for (i = v.GetLwb(); i <= v.GetUpb(); i++)
ok &= ( v(i) == vb(i) ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "Check that expansion expands by zeros" << std::endl;
const Int_t old_nelems = v.GetNoElements();
const Int_t old_lwb = v.GetLwb();
v.ResizeTo(vb);
ok &= ( !(v == vb) ) ? kTRUE : kFALSE;
for (i = old_lwb; i < old_lwb+old_nelems; i++)
ok &= ( v(i) == vb(i) ) ? kTRUE : kFALSE;
for (i = v.GetLwb(); i < old_lwb; i++)
ok &= ( v(i) == 0 ) ? kTRUE : kFALSE;
for (i = old_lwb+old_nelems; i <= v.GetUpb(); i++)
ok &= ( v(i) == 0 ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(1,"Allocation, Filling, Resizing",ok);
}
//
//------------------------------------------------------------------------
// Test uniform element operations
//
class SinAction : public TElementActionD {
void Operation(Double_t &element) const { element = TMath::Sin(element); }
public:
SinAction() { }
};
class CosAction : public TElementPosActionD {
Double_t factor;
void Operation(Double_t &element) const { element = TMath::Cos(factor*fI); }
public:
CosAction(Int_t no_elems): factor(2*TMath::Pi()/no_elems) { }
};
void vstress_element_op(Int_t vsize)
{
if (gVerbose)
std::cout << "\n---> Test operations that treat each element uniformly" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = TMath::Pi();
TVectorD v(-1,vsize-2);
TVectorD v1(v);
if (gVerbose)
std::cout << "\nWriting zeros to v..." << std::endl;
for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = 0;
ok &= VerifyVectorValue(v,0.0,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Clearing v1 ..." << std::endl;
v1.Zero();
ok &= VerifyVectorValue(v1,0.0,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing v1 with 0 ..." << std::endl;
ok &= (v1 == 0) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "Writing a pattern " << pattern << " by assigning to v(i)..." << std::endl;
{
for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = pattern;
ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "Writing the pattern by assigning to v1 as a whole ..." << std::endl;
v1 = pattern;
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "Comparing v and v1 ..." << std::endl;
ok &= (v == v1) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "Comparing (v=0) and v1 ..." << std::endl;
ok &= (!(v.Zero() == v1)) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nClear v and add the pattern" << std::endl;
v.Zero();
v += pattern;
ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " add the doubled pattern with the negative sign" << std::endl;
v += -2*pattern;
ok &= VerifyVectorValue(v,-pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtract the trippled pattern with the negative sign" << std::endl;
v -= -3*pattern;
ok &= VerifyVectorValue(v,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify comparison operations" << std::endl;
v = pattern;
ok &= ( v == pattern && !(v != pattern) && v >= pattern && v <= pattern ) ? kTRUE : kFALSE;
ok &= ( v > 0 && v >= 0 ) ? kTRUE : kFALSE;
ok &= ( v > -pattern && v >= -pattern ) ? kTRUE : kFALSE;
ok &= ( v < pattern+1 && v <= pattern+1 ) ? kTRUE : kFALSE;
v(v.GetUpb()) += 1;
ok &= ( !(v==pattern) && !(v != pattern) && v != pattern-1 ) ? kTRUE : kFALSE;
ok &= ( v >= pattern && !(v > pattern) && !(v >= pattern+1) ) ? kTRUE : kFALSE;
ok &= ( v <= pattern+1.001 && !(v < pattern+1) && !(v <= pattern) ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nAssign 2*pattern to v by repeating additions" << std::endl;
v = 0; v += pattern; v += pattern;
if (gVerbose)
std::cout << "Assign 2*pattern to v1 by multiplying by two" << std::endl;
v1 = pattern; v1 *= 2;
ok &= VerifyVectorValue(v1,2*pattern,gVerbose,EPSILON);
ok &= ( v == v1 ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "Multiply v1 by one half returning it to the 1*pattern" << std::endl;
v1 *= 1/2.;
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nAssign -pattern to v and v1" << std::endl;
v.Zero(); v -= pattern; v1 = -pattern;
ok &= VerifyVectorValue(v,-pattern,gVerbose,EPSILON);
ok &= ( v == v1 ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same" << std::endl;
v.Sqr();
ok &= VerifyVectorValue(v,pattern*pattern,gVerbose,EPSILON);
v.Sqrt();
ok &= VerifyVectorValue(v,pattern,gVerbose,EPSILON);
v1.Abs();
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
ok &= ( v == v1 ) ? kTRUE : kFALSE;
{
if (gVerbose)
std::cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1" << std::endl;
for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = 2*TMath::Pi()/v.GetNoElements()*i;
#ifndef __CINT__
SinAction s;
v.Apply(s);
CosAction c(v.GetNoElements());
v1.Apply(c);
#else
for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++) {
v(i) = TMath::Sin(v(i));
v1(i) = TMath::Cos(2*TMath::Pi()/v1.GetNrows()*i);
}
#endif
TVectorD v2 = v;
TVectorD v3 = v1;
v.Sqr();
v1.Sqr();
v += v1;
ok &= VerifyVectorValue(v,1.,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nVerify constructor with initialization" << std::endl;
#ifndef __CINT__
TVectorD vi(0,4,0.0,1.0,2.0,3.0,4.0,"END");
#else
Double_t vval[] = {0.0,1.0,2.0,3.0,4.0};
TVectorD vi(5,vval);
#endif
TVectorD vit(5);
{
for (Int_t i = vit.GetLwb(); i <= vit.GetUpb(); i++)
vit(i) = Double_t(i);
ok &= VerifyVectorIdentity(vi,vit,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(2,"Uniform vector operations",ok);
}
//
//------------------------------------------------------------------------
// Test binary vector operations
//
void vstress_binary_op(Int_t vsize)
{
if (gVerbose)
std::cout << "\n---> Test Binary Vector operations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = TMath::Pi();
const Double_t epsilon = EPSILON*vsize/10;
TVectorD v(2,vsize+1);
TVectorD v1(v);
if (gVerbose)
std::cout << "\nVerify assignment of a vector to the vector" << std::endl;
v = pattern;
v1.Zero();
v1 = v;
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
ok &= ( v1 == v ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nAdding one vector to itself, uniform pattern " << pattern << std::endl;
v.Zero(); v = pattern;
v1 = v; v1 += v1;
ok &= VerifyVectorValue(v1,2*pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting two vectors ..." << std::endl;
v1 -= v;
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
if (gVerbose)
std::cout << " subtracting the vector from itself" << std::endl;
v1 -= v1;
ok &= VerifyVectorValue(v1,0.,gVerbose,EPSILON);
if (gVerbose)
std::cout << " adding two vectors together" << std::endl;
v1 += v;
ok &= VerifyVectorValue(v1,pattern,gVerbose,EPSILON);
TVectorD vp(2,vsize+1);
{
for (Int_t i = vp.GetLwb(); i <= vp.GetUpb(); i++)
vp(i) = (i-vp.GetNoElements()/2.)*pattern;
}
if (gVerbose) {
std::cout << "\nArithmetic operations on vectors with not the same elements" << std::endl;
std::cout << " adding vp to the zero vector..." << std::endl;
}
v.Zero();
ok &= ( v == 0.0 ) ? kTRUE : kFALSE;
v += vp;
ok &= VerifyVectorIdentity(v,vp,gVerbose,epsilon);
v1 = v;
if (gVerbose)
std::cout << " making v = 3*vp and v1 = 3*vp, via add() and succesive mult" << std::endl;
Add(v,2.,vp);
v1 += v1; v1 += vp;
ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
if (gVerbose)
std::cout << " clear both v and v1, by subtracting from itself and via add()" << std::endl;
v1 -= v1;
Add(v,-3.,vp);
ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
if (gVerbose) {
std::cout << "\nTesting element-by-element multiplications and divisions" << std::endl;
std::cout << " squaring each element with sqr() and via multiplication" << std::endl;
}
v = vp; v1 = vp;
v.Sqr();
ElementMult(v1,v1);
ok &= VerifyVectorIdentity(v,v1,gVerbose,epsilon);
if (gVerbose)
std::cout << " compare (v = pattern^2)/pattern with pattern" << std::endl;
v = pattern; v1 = pattern;
v.Sqr();
ElementDiv(v,v1);
ok &= VerifyVectorValue(v,pattern,gVerbose,epsilon);
if (gVerbose)
Compare(v1,v);
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(3,"Binary vector element-by-element operations",ok);
}
//
//------------------------------------------------------------------------
// Verify the norm calculation
//
void vstress_norms(Int_t vsize)
{
if (gVerbose)
std::cout << "\n---> Verify norm calculations" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 10.25;
if ( vsize % 2 == 1 )
Fatal("vstress_norms", "size of the vector to test must be even for this test\n");
TVectorD v(vsize);
TVectorD v1(v);
if (gVerbose)
std::cout << "\nAssign " << pattern << " to all the elements and check norms" << std::endl;
v = pattern;
if (gVerbose)
std::cout << " 1. norm should be pattern*no_elems" << std::endl;
ok &= ( v.Norm1() == pattern*v.GetNoElements() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Square of the 2. norm has got to be pattern^2 * no_elems" << std::endl;
ok &= ( v.Norm2Sqr() == (pattern*pattern)*v.GetNoElements() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Inf norm should be pattern itself" << std::endl;
ok &= ( v.NormInf() == pattern ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Scalar product of vector by itself is the sqr(2. vector norm)" << std::endl;
ok &= ( v.Norm2Sqr() == v*v ) ? kTRUE : kFALSE;
Double_t ap_step = 1;
Double_t ap_a0 = -pattern;
Int_t n = v.GetNoElements();
if (gVerbose) {
std::cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
"\nand the difference " << ap_step << std::endl;
}
{
for (Int_t i = v.GetLwb(); i <= v.GetUpb(); i++)
v(i) = (i-v.GetLwb())*ap_step + ap_a0;
}
Int_t l = TMath::Min(TMath::Max((int)TMath::Ceil(-ap_a0/ap_step),0),n);
Double_t norm = (2*ap_a0+(l+n-1)*ap_step)/2*(n-l) +
(-2*ap_a0-(l-1)*ap_step)/2*l;
if (gVerbose)
std::cout << " 1. norm should be " << norm << std::endl;
ok &= ( v.Norm1() == norm ) ? kTRUE : kFALSE;
norm = n*( (ap_a0*ap_a0)+ap_a0*ap_step*(n-1)+(ap_step*ap_step)*(n-1)*(2*n-1)/6);
if (gVerbose) {
std::cout << " Square of the 2. norm has got to be "
"n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << std::endl;
}
ok &= ( TMath::Abs( (v.Norm2Sqr()-norm)/norm ) < EPSILON ) ? kTRUE : kFALSE;
norm = TMath::Max(TMath::Abs(v(v.GetLwb())),TMath::Abs(v(v.GetUpb())));
if (gVerbose)
std::cout << " Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm << std::endl;
ok &= ( v.NormInf() == norm ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << " Scalar product of vector by itself is the sqr(2. vector norm)" << std::endl;
ok &= ( v.Norm2Sqr() == v*v ) ? kTRUE : kFALSE;
#if 0
v1.Zero();
Compare(v,v1); // they are not equal (of course)
#endif
if (gVerbose)
std::cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)..." << std::endl;
{
for (Int_t i = 0; i < v1.GetNoElements(); i++)
v1(i+v1.GetLwb()) = v(v.GetUpb()-i) * ( i % 2 == 1 ? -1 : 1 );
}
if (gVerbose)
std::cout << "||v1|| has got to be equal ||v|| regardless of the norm def" << std::endl;
ok &= ( v1.Norm1() == v.Norm1() ) ? kTRUE : kFALSE;
ok &= ( v1.Norm2Sqr() == v.Norm2Sqr() ) ? kTRUE : kFALSE;
ok &= ( v1.NormInf() == v.NormInf() ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "But the scalar product has to be zero" << std::endl;
ok &= ( v1 * v == 0 ) ? kTRUE : kFALSE;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(4,"Vector Norms",ok);
}
//
//------------------------------------------------------------------------
// Test operations with vectors and matrix slices
//
void vstress_matrix_slices(Int_t vsize)
{
if (gVerbose)
std::cout << "\n---> Test operations with vectors and matrix slices" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = 8.625;
TVectorD vc(0,vsize);
TVectorD vr(0,vsize+1);
TMatrixD m(0,vsize,0,vsize+1);
Int_t i,j;
if (gVerbose)
std::cout << "\nCheck modifying the matrix column-by-column" << std::endl;
m = pattern;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
TMatrixDColumn(m,i) = pattern-1;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
TMatrixDColumn(m,i) *= 2;
vc = TMatrixDColumn(m,i);
ok &= VerifyVectorValue(vc,2*(pattern-1),gVerbose);
vc = TMatrixDColumn(m,i+1 > m.GetColUpb() ? m.GetColLwb() : i+1);
ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
TMatrixDColumn(m,i) *= 0.5;
TMatrixDColumn(m,i) += 1;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
vc = pattern+1;
TMatrixDColumn(m,i) = vc;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
{
TMatrixDColumn mc(m,i);
for (j = m.GetRowLwb(); j <= m.GetRowUpb(); j++)
mc(j) *= 4;
}
vc = TMatrixDColumn(m,i);
ok &= VerifyVectorValue(vc,4*(pattern+1),gVerbose,EPSILON);
TMatrixDColumn(m,i) *= 0.25;
TMatrixDColumn(m,i) += -1;
vc = TMatrixDColumn(m,i);
ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nCheck modifying the matrix row-by-row" << std::endl;
m = pattern;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
TMatrixDRow(m,i) = pattern+2;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
vr = TMatrixDRow(m,i);
ok &= VerifyVectorValue(vr,pattern+2,gVerbose,EPSILON);
vr = TMatrixDRow(m,i+1 > m.GetRowUpb() ? m.GetRowLwb() : i+1);
ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
TMatrixDRow(m,i) += -2;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
ok &= ( m == pattern ) ? kTRUE : kFALSE;
for (i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
vr = pattern-2;
TMatrixDRow(m,i) = vr;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
{
TMatrixDRow mr(m,i);
for (j = m.GetColLwb(); j <= m.GetColUpb(); j++)
mr(j) *= 8;
}
vr = TMatrixDRow(m,i);
ok &= VerifyVectorValue(vr,8*(pattern-2),gVerbose,EPSILON);
TMatrixDRow(m,i) *= 1./8;
TMatrixDRow(m,i) += 2;
vr = TMatrixDRow(m,i);
ok &= VerifyVectorValue(vr,pattern,gVerbose,EPSILON);
ok &= ( m == pattern ) ? kTRUE : kFALSE;
}
if (gVerbose)
std::cout << "\nCheck modifying the matrix diagonal" << std::endl;
m = pattern;
TMatrixDDiag td = m;
td = pattern-3;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
vc = TMatrixDDiag(m);
ok &= VerifyVectorValue(vc,pattern-3,gVerbose,EPSILON);
td += 3;
ok &= ( m == pattern ) ? kTRUE : kFALSE;
vc = pattern+3;
td = vc;
ok &= ( !( m == pattern ) && !( m != pattern ) ) ? kTRUE : kFALSE;
{
TMatrixDDiag md(m);
for (j = 0; j < md.GetNdiags(); j++)
md(j) /= 1.5;
}
vc = TMatrixDDiag(m);
ok &= VerifyVectorValue(vc,(pattern+3)/1.5,gVerbose,EPSILON);
TMatrixDDiag(m) *= 1.5;
TMatrixDDiag(m) += -3;
vc = TMatrixDDiag(m);
ok &= VerifyVectorValue(vc,pattern,gVerbose,EPSILON);
ok &= ( m == pattern ) ? kTRUE : kFALSE;
if (gVerbose) {
std::cout << "\nCheck out to see that multiplying by diagonal is column-wise"
"\nmatrix multiplication" << std::endl;
}
TMatrixD mm(m);
TMatrixD m1(m.GetRowLwb(),TMath::Max(m.GetRowUpb(),m.GetColUpb()),
m.GetColLwb(),TMath::Max(m.GetRowUpb(),m.GetColUpb()));
TVectorD vc1(vc),vc2(vc);
for (i = m.GetRowLwb(); i < m.GetRowUpb(); i++)
TMatrixDRow(m,i) = pattern+i; // Make a multiplicand
mm = m; // Save it
m1 = pattern+10;
for (i = vr.GetLwb(); i <= vr.GetUpb(); i++)
vr(i) = i+2;
TMatrixDDiag td2 = m1;
td2 = vr;
ok &= ( !(m1 == pattern+10) ) ? kTRUE : kFALSE;
m *= TMatrixDDiag(m1);
for (i = m.GetColLwb(); i <= m.GetColUpb(); i++) {
vc1 = TMatrixDColumn(mm,i);
vc1 *= vr(i); // Do a column-wise multiplication
vc2 = TMatrixDColumn(m,i);
ok &= VerifyVectorIdentity(vc1,vc2,gVerbose,EPSILON);
}
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(5,"Matrix Slices to Vectors",ok);
}
//
//------------------------------------------------------------------------
// Test vector I/O
//
void vstress_vector_io()
{
if (gVerbose)
std::cout << "\n---> Test vector I/O" << std::endl;
Bool_t ok = kTRUE;
const Double_t pattern = TMath::Pi();
TFile *f = new TFile("stress-vvector.root","RECREATE");
Char_t name[80];
Int_t iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TVectorD v(msize);
v = pattern;
Double_t *pattern_array = new Double_t[msize];
for (Int_t i = 0; i < msize; i++)
pattern_array[i] = pattern;
TVectorD va;
va.Use(msize,pattern_array);
if (verbose)
std::cout << "\nWrite vector v to database" << std::endl;
snprintf(name,80,"v_%d",msize);
v.Write(name);
snprintf(name,80,"va_%d",msize);
va.Write(name);
delete [] pattern_array;
iloop--;
}
if (gVerbose)
std::cout << "\nClose database" << std::endl;
delete f;
if (gVerbose)
std::cout << "\nOpen database in read-only mode and read vector" << std::endl;
TFile *f1 = new TFile("stress-vvector.root");
iloop = gNrLoop;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && iloop==gNrLoop);
TVectorD v(msize);
v = pattern;
snprintf(name,80,"v_%d",msize);
TVectorD *vr = (TVectorD*) f1->Get(name);
snprintf(name,80,"va_%d",msize);
TVectorD *var = (TVectorD*) f1->Get(name);
if (verbose)
std::cout << "\nRead vector should be same as original still in memory" << std::endl;
ok &= ((*vr) == v) ? kTRUE : kFALSE;
ok &= ((*var) == v) ? kTRUE : kFALSE;
iloop--;
}
delete f1;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(6,"Vector Persistence",ok);
}
Bool_t test_svd_expansion(const TMatrixD &A)
{
if (gVerbose)
std::cout << "\n\nSVD-decompose matrix A and check if we can compose it back\n" <<std::endl;
Bool_t ok = kTRUE;
TDecompSVD svd(A);
if (gVerbose) {
std::cout << "left factor U" <<std::endl;
svd.GetU().Print();
std::cout << "Vector of Singular values" <<std::endl;
svd.GetSig().Print();
std::cout << "right factor V" <<std::endl;
svd.GetV().Print();
}
{
if (gVerbose)
std::cout << "\tchecking that U is orthogonal indeed, i.e., U'U=E and UU'=E" <<std::endl;
const Int_t nRows = svd.GetU().GetNrows();
const Int_t nCols = svd.GetU().GetNcols();
TMatrixD E1(nRows,nRows); E1.UnitMatrix();
TMatrixD E2(nCols,nCols); E2.UnitMatrix();
TMatrixD ut(TMatrixD::kTransposed,svd.GetU());
ok &= VerifyMatrixIdentity(ut * svd.GetU(),E2,gVerbose,100*EPSILON);
ok &= VerifyMatrixIdentity(svd.GetU() * ut,E1,gVerbose,100*EPSILON);
}
{
if (gVerbose)
std::cout << "\tchecking that V is orthogonal indeed, i.e., V'V=E and VV'=E" <<std::endl;
const Int_t nRows = svd.GetV().GetNrows();
const Int_t nCols = svd.GetV().GetNcols();
TMatrixD E1(nRows,nRows); E1.UnitMatrix();
TMatrixD E2(nCols,nCols); E2.UnitMatrix();
TMatrixD vt(TMatrixD::kTransposed,svd.GetV());
ok &= VerifyMatrixIdentity(vt * svd.GetV(),E2,gVerbose,100*EPSILON);
ok &= VerifyMatrixIdentity(svd.GetV() * vt,E1,gVerbose,100*EPSILON);
}
{
if (gVerbose)
std::cout << "\tchecking that U*Sig*V' is indeed A" <<std::endl;
const Int_t nRows = svd.GetU().GetNrows();
const Int_t nCols = svd.GetV().GetNcols();
TMatrixD s(nRows,nCols);
TMatrixDDiag diag(s); diag = svd.GetSig();
TMatrixD vt(TMatrixD::kTransposed,svd.GetV());
TMatrixD tmp = s * vt;
ok &= VerifyMatrixIdentity(A,svd.GetU() * tmp,gVerbose,100*EPSILON);
if (gVerbose) {
std::cout << "U*Sig*V'" <<std::endl;
(svd.GetU()*tmp).Print();
}
}
return ok;
}
#ifndef __CINT__
// Make a matrix from an array (read it row-by-row)
class MakeMatrix : public TMatrixDLazy {
const Double_t *array;
Int_t no_elems;
void FillIn(TMatrixD& m) const {
R__ASSERT( m.GetNrows() * m.GetNcols() == no_elems );
const Double_t *ap = array;
Double_t *mp = m.GetMatrixArray();
for (Int_t i = 0; i < no_elems; i++)
*mp++ = *ap++;
}
public:
MakeMatrix(Int_t nrows,Int_t ncols,
const Double_t *_array,Int_t _no_elems)
:TMatrixDLazy(nrows,ncols), array(_array), no_elems(_no_elems) {}
MakeMatrix(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
const Double_t *_array,Int_t _no_elems)
: TMatrixDLazy(row_lwb,row_upb,col_lwb,col_upb),
array(_array), no_elems(_no_elems) {}
};
#else
TMatrixD MakeMatrix(Int_t nrows,Int_t ncols,
const Double_t *_array,Int_t _no_elems)
{
TMatrixD m(nrows,ncols,_array);
return m;
}
#endif
void astress_decomp()
{
Bool_t ok = kTRUE;
{
if (gVerbose)
std::cout << "\nBrandt example page 503\n" <<std::endl;
Double_t array0[] = { -2,1,0,0, 2,1,0,0, 0,0,0,0, 0,0,0,0 };
ok &= test_svd_expansion(MakeMatrix(4,4,array0,sizeof(array0)/sizeof(array0[0])));
}
{
if (gVerbose)
std::cout << "\nRotated by PI/2 Matrix Diag(1,4,9)\n" <<std::endl;
Double_t array1[] = {0,-4,0, 1,0,0, 0,0,9 };
ok &= test_svd_expansion(MakeMatrix(3,3,array1,sizeof(array1)/sizeof(array1[0])));
}
{
if (gVerbose)
std::cout << "\nExample from the Forsythe, Malcolm, Moler's book\n" <<std::endl;
Double_t array2[] =
{ 1,6,11, 2,7,12, 3,8,13, 4,9,14, 5,10,15};
ok &= test_svd_expansion(MakeMatrix(5,3,array2,sizeof(array2)/sizeof(array2[0])));
}
{
if (gVerbose)
std::cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
"Singular numbers are 0, 19.5959, 20, 0, 35.3270\n" <<std::endl;
Double_t array3[] =
{ 22, 10, 2, 3, 7, 14, 7, 10, 0, 8,
-1, 13, -1, -11, 3, -3, -2, 13, -2, 4,
9, 8, 1, -2, 4, 9, 1, -7, 5, -1,
2, -6, 6, 5, 1, 4, 5, 0, -2, 2 };
ok &= test_svd_expansion(MakeMatrix(8,5,array3,sizeof(array3)/sizeof(array3[0])));
}
{
if (gVerbose)
std::cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
"Ordered singular numbers are Sig[21-k] = sqrt(k*(k-1))\n" <<std::endl;
TMatrixD A(21,20);
for (Int_t irow = A.GetRowLwb(); irow <= A.GetRowUpb(); irow++)
for (Int_t icol = A.GetColLwb(); icol <= A.GetColUpb(); icol++)
A(irow,icol) = (irow>icol ? 0 : ((irow==icol) ? 20-icol : -1));
ok &= test_svd_expansion(A);
}
if (0)
{
if (gVerbose) {
std::cout << "\nTest by W. Meier <wmeier@manu.com> to catch an obscure "
<< "bug in QR\n" <<std::endl;
std::cout << "expect singular values to be\n"
<< "1.4666e-024 1.828427 3.828427 4.366725 7.932951\n" <<std::endl;
}
Double_t array4[] =
{ 1, 2, 0, 0, 0,
0, 2, 3, 0, 0,
0, 0, 0, 4, 0,
0, 0, 0, 4, 5,
0, 0, 0, 0, 5 };
ok &= test_svd_expansion(MakeMatrix(5,5,array4,sizeof(array4)/sizeof(array4[0])));
}
{
const TMatrixD m1 = THilbertMatrixD(5,5);
TDecompLU lu(m1);
ok &= VerifyMatrixIdentity(lu.GetMatrix(),m1,gVerbose,100*EPSILON);
}
{
const TMatrixD m2 = THilbertMatrixD(5,5);
const TMatrixDSym mtm(TMatrixDSym::kAtA,m2);
TDecompChol chol(mtm);
ok &= VerifyMatrixIdentity(chol.GetMatrix(),mtm,gVerbose,100*EPSILON);
}
if (gVerbose)
std::cout << "\nDone" <<std::endl;
StatusPrint(1,"Decomposition / Reconstruction",ok);
}
void astress_lineqn()
{
if (gVerbose)
std::cout << "\nSolve Ax=b where A is a Hilbert matrix and b(i) = sum_j Aij\n" <<std::endl;
Bool_t ok = kTRUE;
Int_t iloop = gNrLoop;
Int_t nr = 0;
while (iloop >= 0) {
const Int_t msize = gSizeA[iloop];
const Int_t verbose = (gVerbose && nr==0 && iloop==gNrLoop);
if (verbose)
std::cout << "\nSolve Ax=b for size = " << msize <<std::endl;
// Since The Hilbert matrix is accuracy "challenged", I will use a diagonaly
// dominant one fore sizes > 100, otherwise the verification might fail
TMatrixDSym m1 = THilbertMatrixDSym(-1,msize-2);
TMatrixDDiag diag1(m1);
diag1 += 1.;
TVectorD rowsum1(-1,msize-2); rowsum1.Zero();
TVectorD colsum1(-1,msize-2); colsum1.Zero();
for (Int_t irow = m1.GetRowLwb(); irow <= m1.GetColUpb(); irow++) {
for (Int_t icol = m1.GetColLwb(); icol <= m1.GetColUpb(); icol++) {
rowsum1(irow) += m1(irow,icol);
colsum1(icol) += m1(irow,icol);
}
}
TMatrixDSym m2 = THilbertMatrixDSym(msize);
TMatrixDDiag diag2(m2);
diag2 += 1.;
TVectorD rowsum2(msize); rowsum2.Zero();
TVectorD colsum2(msize); colsum2.Zero();
for (Int_t irow = m2.GetRowLwb(); irow <= m2.GetColUpb(); irow++) {
for (Int_t icol = m2.GetColLwb(); icol <= m2.GetColUpb(); icol++) {
rowsum2(irow) += m2(irow,icol);
colsum2(icol) += m2(irow,icol);
}
}
TVectorD b1(-1,msize-2);
TVectorD b2(msize);
{
TDecompLU lu(m1,1.0e-20);
b1 = rowsum1;
lu.Solve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
b1 = colsum1;
lu.TransSolve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
}
{
TDecompChol chol(m1,1.0e-20);
b1 = rowsum1;
chol.Solve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
b1 = colsum1;
chol.TransSolve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
}
{
TDecompQRH qrh1(m1,1.0e-20);
b1 = rowsum1;
qrh1.Solve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
TDecompQRH qrh2(m2,1.0e-20);
b2 = colsum2;
qrh2.TransSolve(b2);
if (msize < 10)
ok &= VerifyVectorValue(b2,1.0,verbose,msize*EPSILON);
}
{
TDecompSVD svd1(m1);
b1 = rowsum1;
svd1.Solve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
TDecompSVD svd2(m2);
b2 = colsum2;
svd2.TransSolve(b2);
if (msize < 10)
ok &= VerifyVectorValue(b2,1.0,verbose,msize*EPSILON);
}
{
TDecompBK bk(m1,1.0e-20);
b1 = rowsum1;
bk.Solve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
b1 = colsum1;
bk.TransSolve(b1);
if (msize < 10)
ok &= VerifyVectorValue(b1,1.0,verbose,msize*EPSILON);
}
#ifndef __CINT__
if (nr >= Int_t(1.0e+5/msize/msize)) {
#else
if (nr >= Int_t(1.0e+3/msize/msize)) {
#endif
nr = 0;
iloop--;
} else
nr++;
if (!ok) {
if (gVerbose)
std::cout << "Linear Equations failed for size " << msize << std::endl;
break;
}
}
if (gVerbose)
std::cout << "\nDone" <<std::endl;
StatusPrint(2,"Linear Equations",ok);
}
void astress_pseudo()
{
// The pseudo-inverse of A is found by "inverting" the SVD of A.
// To be more precise, we use SVD to solve the equation
// AX=B where B is a unit matrix.
//
// After we compute the pseudo-inverse, we verify the Moore-Penrose
// conditions: Matrix X is a pseudo-inverse of A iff
// AXA = A
// XAX = X
// XA = (XA)' (i.e., XA is symmetric)
// AX = (AX)' (i.e., AX is symmetric)
Bool_t ok = kTRUE;
// Allocate and fill matrix A
enum {nrows = 4, ncols = 3};
#ifndef __CINT__
const Double_t A_data [nrows*ncols] =
#else
const Double_t A_data [12] =
#endif
{0, 0, 0,
0, 0, 0,
1, 1, 0,
4, 3, 0};
TMatrixD A(nrows,ncols,A_data);
// Perform the SVD decomposition of the transposed matrix.
TDecompSVD svd(A);
if (gVerbose)
std::cout << "\ncondition number " << svd.Condition() << std::endl;
// Compute the inverse as a solution of A*Ainv = E
TMatrixD Ainv(nrows,nrows); Ainv.UnitMatrix();
svd.MultiSolve(Ainv);
Ainv.ResizeTo(ncols,nrows);
TMatrixD Ainv2(nrows,nrows);
svd.Invert(Ainv2);
Ainv2.ResizeTo(ncols,nrows);
ok &= VerifyMatrixIdentity(Ainv,Ainv2,gVerbose,EPSILON);
if (gVerbose) {
std::cout << "\nChecking the Moore-Penrose conditions for the Ainv" << std::endl;
std::cout << "\nVerify that Ainv * A * Ainv is Ainv" << std::endl;
}
ok &= VerifyMatrixIdentity(Ainv * A * Ainv, Ainv,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify that A * Ainv * A is A" << std::endl;
ok &= VerifyMatrixIdentity(A * Ainv * A, A,gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify that Ainv * A is symmetric" << std::endl;
ok &= VerifyMatrixIdentity(Ainv * A, (Ainv*A).T(),gVerbose,EPSILON);
if (gVerbose)
std::cout << "\nVerify that A * Ainv is symmetric" << std::endl;
ok &= VerifyMatrixIdentity(A * Ainv, (A*Ainv).T(),gVerbose,EPSILON);
StatusPrint(3,"Pseudo-Inverse, Moore-Penrose",ok);
}
void astress_eigen(Int_t msize)
{
Bool_t ok = kTRUE;
// Check :
// 1. the eigen values of M' * M are the same as the singular values
// squared of the SVD of M .
// 2. M' * M x_i = lambda_i x_i where x_i and lambda_i are the ith
// eigen - vector/value .
const TMatrixD m = THilbertMatrixD(msize,msize);
TDecompSVD svd(m);
TVectorD sig = svd.GetSig(); sig.Sqr();
// Symmetric matrix EigenVector algorithm
TMatrixDSym mtm1(TMatrixDSym::kAtA,m);
const TMatrixDSymEigen eigen1(mtm1);
const TVectorD eigenVal1 = eigen1.GetEigenValues();
ok &= VerifyVectorIdentity(sig,eigenVal1,gVerbose,EPSILON);
// Check that the eigen vectors by comparing M'M x to lambda x
const TMatrixD mtm1x = mtm1 * eigen1.GetEigenVectors();
TMatrixD lam_x1 = eigen1.GetEigenVectors();
lam_x1.NormByRow(eigen1.GetEigenValues(),"");
ok &= VerifyMatrixIdentity(mtm1x,lam_x1,gVerbose,EPSILON);
// General matrix EigenVector algorithm tested on symmetric matrix
TMatrixD mtm2(m,TMatrixD::kTransposeMult,m);
const TMatrixDEigen eigen2(mtm2);
const TVectorD eigenVal2 = eigen2.GetEigenValuesRe();
ok &= VerifyVectorIdentity(sig,eigenVal2,gVerbose,EPSILON);
// Check that the eigen vectors by comparing M'M x to lambda x
const TMatrixD mtm2x = mtm2 * eigen2.GetEigenVectors();
TMatrixD lam_x2 = eigen2.GetEigenVectors();
lam_x2.NormByRow(eigen2.GetEigenValuesRe(),"");
ok &= VerifyMatrixIdentity(mtm2x,lam_x2,gVerbose,EPSILON);
// The Imaginary part of the eigenvalues should be zero
const TVectorD eigenValIm = eigen2.GetEigenValuesIm();
TVectorD epsVec(msize); epsVec = EPSILON;
ok &= VerifyVectorIdentity(epsVec,eigenValIm,gVerbose,EPSILON);
StatusPrint(4,"Eigen - Values/Vectors",ok);
}
void astress_decomp_io(Int_t msize)
{
//
//------------------------------------------------------------------------
// Test decomposition I/O
//
if (gVerbose)
std::cout << "\n---> Test decomp I/O" << std::endl;
Bool_t ok = kTRUE;
const TMatrixDSym m = THilbertMatrixDSym(msize);
TVectorD rowsum(msize); rowsum.Zero();
TVectorD colsum(msize); colsum.Zero();
for (Int_t irow = 0; irow < m.GetNrows(); irow++) {
for (Int_t icol = 0; icol < m.GetNcols(); icol++) {
rowsum(irow) += m(irow,icol);
colsum(icol) += m(irow,icol);
}
}
if (gVerbose)
std::cout << "\nWrite decomp m to database" << std::endl;
TFile *f = new TFile("stress-vdecomp.root", "RECREATE");
TDecompLU lu(m,1.0e-20);
TDecompQRH qrh(m,1.0e-20);
TDecompChol chol(m,1.0e-20);
TDecompSVD svd(m);
TDecompBK bk(m,1.0e-20);
lu.Write("lu");
qrh.Write("qrh");
chol.Write("chol");
svd.Write("svd");
bk.Write("bk");
if (gVerbose)
std::cout << "\nClose database" << std::endl;
delete f;
if (gVerbose)
std::cout << "\nOpen database in read-only mode and read matrix" << std::endl;
TFile *f1 = new TFile("stress-vdecomp.root");
if (gVerbose)
std::cout << "\nRead decompositions should create same solutions" << std::endl;
{
TDecompLU *rlu = (TDecompLU*) f1->Get("lu");
TVectorD b1(rowsum);
lu.Solve(b1);
TVectorD b2(rowsum);
rlu->Solve(b2);
ok &= (b1 == b2);
b1 = colsum;
lu.TransSolve(b1);
b2 = colsum;
rlu->TransSolve(b2);
ok &= (b1 == b2);
}
{
TDecompChol *rchol = (TDecompChol*) f1->Get("chol");
TVectorD b1(rowsum);
chol.Solve(b1);
TVectorD b2(rowsum);
rchol->Solve(b2);
ok &= (b1 == b2);
b1 = colsum;
chol.TransSolve(b1);
b2 = colsum;
rchol->TransSolve(b2);
ok &= (b1 == b2);
}
{
TDecompQRH *rqrh = (TDecompQRH*) f1->Get("qrh");
TVectorD b1(rowsum);
qrh.Solve(b1);
TVectorD b2(rowsum);
rqrh->Solve(b2);
ok &= (b1 == b2);
b1 = colsum;
qrh.TransSolve(b1);
b2 = colsum;
rqrh->TransSolve(b2);
ok &= (b1 == b2);
}
{
TDecompSVD *rsvd = (TDecompSVD*) f1->Get("svd");
TVectorD b1(rowsum);
svd.Solve(b1);
TVectorD b2(rowsum);
rsvd->Solve(b2);
ok &= (b1 == b2);
b1 = colsum;
svd.TransSolve(b1);
b2 = colsum;
rsvd->TransSolve(b2);
ok &= (b1 == b2);
}
{
TDecompBK *rbk = (TDecompBK*) f1->Get("bk");
TVectorD b1(rowsum);
bk.Solve(b1);
TVectorD b2(rowsum);
rbk->Solve(b2);
ok &= (b1 == b2);
b1 = colsum;
bk.TransSolve(b1);
b2 = colsum;
rbk->TransSolve(b2);
ok &= (b1 == b2);
}
delete f1;
if (gVerbose)
std::cout << "\nDone\n" << std::endl;
StatusPrint(5,"Decomposition Persistence",ok);
}
void stress_backward_io()
{
TFile::SetCacheFileDir(".");
TFile *f = TFile::Open("http://root.cern.ch/files/linearIO.root","CACHEREAD");
TMatrixF mf1 = THilbertMatrixF(-5,5,-5,5);
mf1[1][2] = TMath::Pi();
TMatrixFSym mf2 = THilbertMatrixFSym(-5,5);
TVectorF vf_row(mf1.GetRowLwb(),mf1.GetRowUpb()); vf_row = TMatrixFRow(mf1,3);
TMatrixF *mf1_r = (TMatrixF*) f->Get("mf1");
TMatrixFSym *mf2_r = (TMatrixFSym*) f->Get("mf2");
TVectorF *vf_row_r = (TVectorF*) f->Get("vf_row");
Bool_t ok = kTRUE;
ok &= ((*mf1_r) == mf1) ? kTRUE : kFALSE;
ok &= ((*mf2_r) == mf2) ? kTRUE : kFALSE;
ok &= ((*vf_row_r) == vf_row) ? kTRUE : kFALSE;
TMatrixD md1 = THilbertMatrixD(-5,5,-5,5);
md1[1][2] = TMath::Pi();
TMatrixDSym md2 = THilbertMatrixDSym(-5,5);
TMatrixDSparse md3 = md1;
TVectorD vd_row(md1.GetRowLwb(),md1.GetRowUpb()); vd_row = TMatrixDRow(md1,3);
TMatrixD *md1_r = (TMatrixD*) f->Get("md1");
TMatrixDSym *md2_r = (TMatrixDSym*) f->Get("md2");
TMatrixDSparse *md3_r = (TMatrixDSparse*) f->Get("md3");
TVectorD *vd_row_r = (TVectorD*) f->Get("vd_row");
ok &= ((*md1_r) == md1) ? kTRUE : kFALSE;
ok &= ((*md2_r) == md2) ? kTRUE : kFALSE;
ok &= ((*md3_r) == md3) ? kTRUE : kFALSE;
ok &= ((*vd_row_r) == vd_row) ? kTRUE : kFALSE;
StatusPrint(1,"Streamers",ok);
}
void cleanup()
{
gSystem->Unlink("stress-vmatrix.root");
gSystem->Unlink("stress-vvector.root");
gSystem->Unlink("stress-vdecomp.root");
}
Computing file changes ...
