// @(#)root/matrix:$Name: $:$Id: TMatrixD.cxx,v 1.80 2005/04/05 12:47:11 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Nov 2003
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
//////////////////////////////////////////////////////////////////////////
// //
// TMatrixD //
// //
// Implementation of a general matrix in the linear algebra package //
// //
//////////////////////////////////////////////////////////////////////////
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixF.h"
#include "TMatrixDLazy.h"
#include "TMatrixDCramerInv.h"
#include "TDecompLU.h"
#include "TMatrixDEigen.h"
ClassImp(TMatrixD)
//______________________________________________________________________________
TMatrixD::TMatrixD(Int_t no_rows,Int_t no_cols)
{
Allocate(no_rows,no_cols,0,0,1);
}
//______________________________________________________________________________
TMatrixD::TMatrixD(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
{
Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1);
}
//______________________________________________________________________________
TMatrixD::TMatrixD(Int_t no_rows,Int_t no_cols,const Double_t *elements,Option_t *option)
{
// option="F": array elements contains the matrix stored column-wise
// like in Fortran, so a[i,j] = elements[i+no_rows*j],
// else it is supposed that array elements are stored row-wise
// a[i,j] = elements[i*no_cols+j]
//
// array elements are copied
Allocate(no_rows,no_cols);
SetMatrixArray(elements,option);
}
//______________________________________________________________________________
TMatrixD::TMatrixD(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
const Double_t *elements,Option_t *option)
{
// array elements are copied
Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb);
SetMatrixArray(elements,option);
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixD &another) : TMatrixDBase(another)
{
Assert(another.IsValid());
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
*this = another;
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixF &another)
{
Assert(another.IsValid());
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
*this = another;
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixDSym &another)
{
Assert(another.IsValid());
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
*this = another;
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixDSparse &another)
{
Assert(another.IsValid());
Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
*this = another;
}
//______________________________________________________________________________
TMatrixD::TMatrixD(EMatrixCreatorsOp1 op,const TMatrixD &prototype)
{
// Create a matrix applying a specific operation to the prototype.
// Example: TMatrixD a(10,12); ...; TMatrixD b(TMatrixDBase::kTransposed, a);
// Supported operations are: kZero, kUnit, kTransposed, kInverted and kAtA.
Assert(this != &prototype);
Invalidate();
Assert(prototype.IsValid());
switch(op) {
case kZero:
Allocate(prototype.GetNrows(),prototype.GetNcols(),
prototype.GetRowLwb(),prototype.GetColLwb(),1);
break;
case kUnit:
Allocate(prototype.GetNrows(),prototype.GetNcols(),
prototype.GetRowLwb(),prototype.GetColLwb(),1);
UnitMatrix();
break;
case kTransposed:
Allocate(prototype.GetNcols(), prototype.GetNrows(),
prototype.GetColLwb(),prototype.GetRowLwb());
Transpose(prototype);
break;
case kInverted:
{
Allocate(prototype.GetNrows(),prototype.GetNcols(),
prototype.GetRowLwb(),prototype.GetColLwb(),1);
*this = prototype;
// Since the user can not control the tolerance of this newly created matrix
// we put it to the smallest possible number
const Double_t oldTol = this->SetTol(DBL_MIN);
this->Invert();
this->SetTol(oldTol);
break;
}
case kAtA:
AtMultB(prototype,prototype);
break;
default:
Error("TMatrixD(EMatrixCreatorOp1)", "operation %d not yet implemented", op);
}
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixD &a,EMatrixCreatorsOp2 op,const TMatrixD &b)
{
// Create a matrix applying a specific operation to two prototypes.
// Example: TMatrixD a(10,12), b(12,5); ...; TMatrixD c(a, TMatrixDBase::kMult, b);
// Supported operations are: kMult (a*b), kTransposeMult (a'*b), kInvMult (a^(-1)*b)
Invalidate();
Assert(a.IsValid());
Assert(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b);
break;
case kTransposeMult:
AtMultB(a,b);
break;
case kMultTranspose:
AMultBt(a,b);
break;
case kInvMult:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
const Double_t oldTol = this->SetTol(DBL_MIN);
this->Invert();
this->SetTol(oldTol);
*this *= b;
break;
}
case kPlus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this += b;
break;
}
case kMinus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this -= b;
break;
}
default:
Error("TMatrixD(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
}
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixD &a,EMatrixCreatorsOp2 op,const TMatrixDSym &b)
{
Invalidate();
Assert(a.IsValid());
Assert(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b);
break;
case kTransposeMult:
AtMultB(a,b);
break;
case kMultTranspose:
AMultBt(a,b);
break;
case kInvMult:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
const Double_t oldTol = this->SetTol(DBL_MIN);
this->Invert();
this->SetTol(oldTol);
*this *= b;
break;
}
case kPlus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this += b;
break;
}
case kMinus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this -= b;
break;
}
default:
Error("TMatrixD(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
}
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixDSym &a,EMatrixCreatorsOp2 op,const TMatrixD &b)
{
Invalidate();
Assert(a.IsValid());
Assert(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b);
break;
case kTransposeMult:
AtMultB(a,b);
break;
case kMultTranspose:
AMultBt(a,b);
break;
case kInvMult:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
const Double_t oldTol = this->SetTol(DBL_MIN);
this->Invert();
this->SetTol(oldTol);
*this *= b;
break;
}
case kPlus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this += b;
break;
}
case kMinus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this -= b;
break;
}
default:
Error("TMatrixD(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
}
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixDSym &a,EMatrixCreatorsOp2 op,const TMatrixDSym &b)
{
Invalidate();
Assert(a.IsValid());
Assert(b.IsValid());
switch(op) {
case kMult:
AMultB(a,b);
break;
case kTransposeMult:
AtMultB(a,b);
break;
case kMultTranspose:
AMultBt(a,b);
break;
case kInvMult:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
const Double_t oldTol = this->SetTol(DBL_MIN);
this->Invert();
this->SetTol(oldTol);
*this *= b;
break;
}
case kPlus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this += b;
break;
}
case kMinus:
{
Allocate(a.GetNrows(),a.GetNcols(),
a.GetRowLwb(),a.GetColLwb(),1);
*this = a;
*this -= b;
break;
}
default:
Error("TMatrixD(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
}
}
//______________________________________________________________________________
TMatrixD::TMatrixD(const TMatrixDLazy &lazy_constructor)
{
Allocate(lazy_constructor.GetRowUpb()-lazy_constructor.GetRowLwb()+1,
lazy_constructor.GetColUpb()-lazy_constructor.GetColLwb()+1,
lazy_constructor.GetRowLwb(),lazy_constructor.GetColLwb(),1);
lazy_constructor.FillIn(*this);
}
//______________________________________________________________________________
void TMatrixD::Delete_m(Int_t size,Double_t *&m)
{
// delete data pointer m, if it was assigned on the heap
if (m) {
if (size > kSizeMax)
delete [] m;
m = 0;
}
}
//______________________________________________________________________________
Double_t* TMatrixD::New_m(Int_t size)
{
// return data pointer . if requested size <= kSizeMax, assign pointer
// to the stack space
if (size == 0) return 0;
else {
if ( size <= kSizeMax )
return fDataStack;
else {
Double_t *heap = new Double_t[size];
return heap;
}
}
}
//______________________________________________________________________________
Int_t TMatrixD::Memcpy_m(Double_t *newp,const Double_t *oldp,Int_t copySize,
Int_t newSize,Int_t oldSize)
{
// copy copySize doubles from *oldp to *newp . However take care of the
// situation where both pointers are assigned to the same stack space
if (copySize == 0 || oldp == newp)
return 0;
else {
if ( newSize <= kSizeMax && oldSize <= kSizeMax ) {
// both pointers are inside fDataStack, be careful with copy direction !
if (newp > oldp) {
for (Int_t i = copySize-1; i >= 0; i--)
newp[i] = oldp[i];
} else {
for (Int_t i = 0; i < copySize; i++)
newp[i] = oldp[i];
}
}
else
memcpy(newp,oldp,copySize*sizeof(Double_t));
}
return 0;
}
//______________________________________________________________________________
void TMatrixD::Allocate(Int_t no_rows,Int_t no_cols,Int_t row_lwb,Int_t col_lwb,
Int_t init,Int_t /*nr_nonzeros*/)
{
// Allocate new matrix. Arguments are number of rows, columns, row
// lowerbound (0 default) and column lowerbound (0 default).
if (no_rows < 0 || no_cols < 0)
{
Error("Allocate","no_rows=%d no_cols=%d",no_rows,no_cols);
Invalidate();
return;
}
MakeValid();
fNrows = no_rows;
fNcols = no_cols;
fRowLwb = row_lwb;
fColLwb = col_lwb;
fNelems = fNrows*fNcols;
fIsOwner = kTRUE;
fTol = DBL_EPSILON;
if (fNelems > 0) {
fElements = New_m(fNelems);
if (init)
memset(fElements,0,fNelems*sizeof(Double_t));
} else
fElements = 0;
}
//______________________________________________________________________________
void TMatrixD::AMultB(const TMatrixD &a,const TMatrixD &b,Int_t constr)
{
// General matrix multiplication. Create a matrix C such that C = A * B.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
Error("AMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == &a) {
Error("AMultB","this == &a");
Invalidate();
return;
}
if (this == &b) {
Error("AMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNcols(),a.GetRowLwb(),b.GetColLwb(),1);
#ifdef CBLAS
const Double_t *ap = a.GetMatrixArray();
const Double_t *bp = b.GetMatrixArray();
Double_t *cp = this->GetMatrixArray();
cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,fNrows,fNcols,a.GetNcols(),
1.0,ap,a.GetNcols(),bp,b.GetNcols(),1.0,cp,fNcols);
#else
const Int_t na = a.GetNoElements();
const Int_t nb = b.GetNoElements();
const Int_t ncolsa = a.GetNcols();
const Int_t ncolsb = b.GetNcols();
const Double_t * const ap = a.GetMatrixArray();
const Double_t * const bp = b.GetMatrixArray();
Double_t * cp = this->GetMatrixArray();
const Double_t *arp0 = ap; // Pointer to A[i,0];
while (arp0 < ap+na) {
for (const Double_t *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of B, Start bcp = B[0,0]
const Double_t *arp = arp0; // Pointer to the i-th row of A, reset to A[i,0]
Double_t cij = 0;
while (bcp < bp+nb) { // Scan the i-th row of A and
cij += *arp++ * *bcp; // the j-th col of B
bcp += ncolsb;
}
*cp++ = cij;
bcp -= nb-1; // Set bcp to the (j+1)-th col
}
arp0 += ncolsa; // Set ap to the (i+1)-th row
}
Assert(cp == this->GetMatrixArray()+fNelems && arp0 == ap+na);
#endif
}
//______________________________________________________________________________
void TMatrixD::AMultB(const TMatrixDSym &a,const TMatrixD &b,Int_t constr)
{
// Matrix multiplication, with A symmetric and B general.
// Create a matrix C such that C = A * B.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
Error("AMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&a)) {
Error("AMultB","this == &a");
Invalidate();
return;
}
if (this == &b) {
Error("AMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNcols(),a.GetRowLwb(),b.GetColLwb(),1);
const Double_t *ap1 = a.GetMatrixArray();
const Double_t *bp1 = b.GetMatrixArray();
Double_t *cp1 = this->GetMatrixArray();
#ifdef CBLAS
cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
ap1,a.GetNcols(),bp1,b.GetNcols(),0.0,cp1,fNcols);
#else
const Double_t *bp2 = b.GetMatrixArray();
Double_t *cp2 = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
for (Int_t j = 0; j < fNcols; j++) {
const Double_t b_ij = *bp1++;
*cp1 += b_ij*(*ap1);
Double_t tmp = 0.0;
const Double_t *ap2 = ap1+1;
for (Int_t k = i+1; k < fNrows; k++) {
const Int_t index_kj = k*fNcols+j;
const Double_t a_ik = *ap2++;
const Double_t b_kj = bp2[index_kj];
cp2[index_kj] += a_ik*b_ij;
tmp += a_ik*b_kj;
}
*cp1++ += tmp;
}
ap1 += fNrows+1;
}
#endif
}
//______________________________________________________________________________
void TMatrixD::AMultB(const TMatrixD &a,const TMatrixDSym &b,Int_t constr)
{
// Matrix multiplication, with A general and B symmetric.
// Create a matrix C such that C = A * B.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
Error("AMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == &a) {
Error("AMultB","this == &a");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&b)) {
Error("AMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNcols(),a.GetRowLwb(),b.GetColLwb(),1);
const Double_t *ap1 = a.GetMatrixArray();
Double_t *cp1 = this->GetMatrixArray();
#ifdef CBLAS
const Double_t *bp1 = b.GetMatrixArray();
cblas_dsymm (CblasRowMajor,CblasRight,CblasUpper,fNrows,fNcols,1.0,
bp1,b.GetNcols(),ap1,a.GetNcols(),0.0,cp1,fNcols);
#else
for (Int_t i = 0; i < fNrows; i++) {
const Double_t *bp1 = b.GetMatrixArray();
for (Int_t j = 0; j < fNcols; j++) {
const Double_t a_ij = *ap1++;
*cp1 += a_ij*(*bp1);
Double_t tmp = 0.0;
const Double_t *ap2 = ap1;
const Double_t *bp2 = bp1+1;
Double_t *cp2 = cp1+1;
for (Int_t k = j+1; k < fNcols; k++) {
const Double_t a_ik = *ap2++;
const Double_t b_jk = *bp2++;
*cp2++ += a_ij*b_jk;
tmp += a_ik*b_jk;
}
*cp1++ += tmp;
bp1 += fNcols+1;
}
}
#endif
}
//______________________________________________________________________________
void TMatrixD::AMultB(const TMatrixDSym &a,const TMatrixDSym &b,Int_t constr)
{
// Matrix multiplication, with A symmetric and B symmetric.
// (Actually copied for the moment routine for B general)
// Create a matrix C such that C = A * B.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
Error("AMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&a)) {
Error("AMultB","this == &a");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&b)) {
Error("AMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNcols(),a.GetRowLwb(),b.GetColLwb(),1);
const Double_t *ap1 = a.GetMatrixArray();
const Double_t *bp1 = b.GetMatrixArray();
Double_t *cp1 = this->GetMatrixArray();
#ifdef CBLAS
cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
ap1,a.GetNcols(),bp1,b.GetNcols(),0.0,cp1,fNcols);
#else
const Double_t *bp2 = b.GetMatrixArray();
Double_t *cp2 = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
for (Int_t j = 0; j < fNcols; j++) {
const Double_t b_ij = *bp1++;
*cp1 += b_ij*(*ap1);
Double_t tmp = 0.0;
const Double_t *ap2 = ap1+1;
for (Int_t k = i+1; k < fNrows; k++) {
const Int_t index_kj = k*fNcols+j;
const Double_t a_ik = *ap2++;
const Double_t b_kj = bp2[index_kj];
cp2[index_kj] += a_ik*b_ij;
tmp += a_ik*b_kj;
}
*cp1++ += tmp;
}
ap1 += fNrows+1;
}
#endif
}
//______________________________________________________________________________
void TMatrixD::AtMultB(const TMatrixD &a,const TMatrixD &b,Int_t constr)
{
// Create a matrix C such that C = A' * B. In other words,
// c[i,j] = SUM{ a[k,i] * b[k,j] }. Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetRowLwb() != b.GetRowLwb()) {
Error("AtMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == &a) {
Error("AtMultB","this == &a");
Invalidate();
return;
}
if (this == &b) {
Error("AtMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNcols(),b.GetNcols(),a.GetColLwb(),b.GetColLwb(),1);
#ifdef CBLAS
const Double_t *ap = a.GetMatrixArray();
const Double_t *bp = b.GetMatrixArray();
Double_t *cp = this->GetMatrixArray();
cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,a.GetNrows(),
1.0,ap,a.GetNcols(),bp,b.GetNcols(),1.0,cp,fNcols);
#else
const Int_t nb = b.GetNoElements();
const Int_t ncolsa = a.GetNcols();
const Int_t ncolsb = b.GetNcols();
const Double_t * const ap = a.GetMatrixArray();
const Double_t * const bp = b.GetMatrixArray();
Double_t * cp = this->GetMatrixArray();
const Double_t *acp0 = ap; // Pointer to A[i,0];
while (acp0 < ap+ncolsa) {
for (const Double_t *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of B, Start bcp = B[0,0]
const Double_t *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
Double_t cij = 0;
while (bcp < bp+nb) { // Scan the i-th column of A and
cij += *acp * *bcp; // the j-th col of B
acp += ncolsa;
bcp += ncolsb;
}
*cp++ = cij;
bcp -= nb-1; // Set bcp to the (j+1)-th col
}
acp0++; // Set acp0 to the (i+1)-th col
}
Assert(cp == this->GetMatrixArray()+fNelems && acp0 == ap+ncolsa);
#endif
}
//______________________________________________________________________________
void TMatrixD::AtMultB(const TMatrixD &a,const TMatrixDSym &b,Int_t constr)
{
// Create a matrix C such that C = A' * B. In other words,
// c[i,j] = SUM{ a[k,i] * b[k,j] }. Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNrows() != b.GetNrows() || a.GetRowLwb() != b.GetRowLwb()) {
Error("AtMultB","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == &a) {
Error("AtMultB","this == &a");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&b)) {
Error("AtMultB","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNcols(),b.GetNcols(),a.GetColLwb(),b.GetColLwb(),1);
#ifdef CBLAS
const Double_t *ap = a.GetMatrixArray();
const Double_t *bp = b.GetMatrixArray();
Double_t *cp = this->GetMatrixArray();
cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,a.GetNrows(),
1.0,ap,a.GetNcols(),bp,b.GetNcols(),1.0,cp,fNcols);
#else
Double_t *cp1 = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
const Double_t *ap1 = a.GetMatrixArray()+i; // i-column of a
const Double_t *bp1 = b.GetMatrixArray();
for (Int_t j = 0; j < fNcols; j++) {
const Double_t a_ji = *ap1;
*cp1 += a_ji*(*bp1);
Double_t tmp = 0.0;
const Double_t *ap2 = ap1+fNrows;
const Double_t *bp2 = bp1+1;
Double_t *cp2 = cp1+1;
for (Int_t k = j+1; k < fNcols; k++) {
const Double_t b_jk = *bp2++;
*cp2++ += a_ji*b_jk;
tmp += (*ap2) * b_jk;
ap2 += fNrows;
}
*cp1++ += tmp;
ap1 += fNrows;
bp1 += fNcols+1;
}
}
#endif
}
//______________________________________________________________________________
void TMatrixD::AMultBt(const TMatrixD &a,const TMatrixD &b,Int_t constr)
{
// General matrix multiplication. Create a matrix C such that C = A * B^T.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
Error("AMultBt","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == &a) {
Error("AMultBt","this == &a");
Invalidate();
return;
}
if (this == &b) {
Error("AMultBt","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1);
#ifdef CBLAS
const Double_t *ap = a.GetMatrixArray();
const Double_t *bp = b.GetMatrixArray();
Double_t *cp = this->GetMatrixArray();
cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,fNrows,fNcols,a.GetNcols(),
1.0,ap,a.GetNcols(),bp,b.GetNcols(),1.0,cp,fNcols);
#else
const Int_t na = a.GetNoElements();
const Int_t nb = b.GetNoElements();
const Int_t ncolsa = a.GetNcols();
const Int_t ncolsb = b.GetNcols();
const Double_t * const ap = a.GetMatrixArray();
const Double_t * const bp = b.GetMatrixArray();
Double_t * cp = this->GetMatrixArray();
const Double_t *arp0 = ap; // Pointer to A[i,0];
while (arp0 < ap+na) {
const Double_t *brp0 = bp; // Pointer to B[j,0];
while (brp0 < bp+nb) {
const Double_t *arp = arp0; // Pointer to the i-th row of A, reset to A[i,0]
const Double_t *brp = brp0; // Pointer to the j-th row of B, reset to B[j,0]
Double_t cij = 0;
while (brp < brp0+ncolsb) // Scan the i-th row of A and
cij += *arp++ * *brp++; // the j-th row of B
*cp++ = cij;
brp0 += ncolsb; // Set brp0 to the (j+1)-th row
}
arp0 += ncolsa; // Set arp0 to the (i+1)-th row
}
Assert(cp == this->GetMatrixArray()+fNelems && arp0 == ap+na);
#endif
}
//______________________________________________________________________________
void TMatrixD::AMultBt(const TMatrixDSym &a,const TMatrixD &b,Int_t constr)
{
// Matrix multiplication, with A symmetric and B general.
// Create a matrix C such that C = A * B^T.
// Note, matrix C is allocated for constr=1.
Assert(a.IsValid());
Assert(b.IsValid());
if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
Error("AMultBt","A rows and B columns incompatible");
Invalidate();
return;
}
if (this == dynamic_cast<const TMatrixD *>(&a)) {
Error("AMultBt","this == &a");
Invalidate();
return;
}
if (this == &b) {
Error("AMultBt","this == &b");
Invalidate();
return;
}
if (constr)
Allocate(a.GetNrows(),b.GetNrows(),a.GetRowLwb(),b.GetRowLwb(),1);
const Double_t *ap1 = a.GetMatrixArray();
const Double_t *bp1 = b.GetMatrixArray();
Double_t *cp1 = this->GetMatrixArray();
#ifdef CBLAS
cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,fNrows,fNcols,a.GetNcols(),
1.0,ap,a.GetNcols(),bp,b.GetNcols(),1.0,cp,fNcols);
#else
const Int_t nb = b.GetNoElements();
const Int_t ncolsb = b.GetNcols();
Double_t *cp2 = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
for (Int_t j = 0; j < fNcols; j++) {
const Double_t b_ji = *bp1;
*cp1 += b_ji*(*ap1);
Double_t tmp = 0.0;
const Double_t *ap2 = ap1+1;
const Double_t *bp2 = bp1+1;
for (Int_t k = i+1; k < fNrows; k++) {
const Int_t index_kj = k*fNcols+j;
const Double_t a_ik = *ap2++;
const Double_t b_jk = *bp2++;
cp2[index_kj] += a_ik*b_ji;
tmp += a_ik*b_jk;
}
*cp1++ += tmp;
bp1 += ncolsb;
}
ap1 += fNrows+1;
bp1 -= nb-1;
}
#endif
}
//______________________________________________________________________________
TMatrixD &TMatrixD::Use(Int_t row_lwb,Int_t row_upb,
Int_t col_lwb,Int_t col_upb,Double_t *data)
{
if (row_upb < row_lwb)
{
Error("Use","row_upb=%d < row_lwb=%d",row_upb,row_lwb);
Invalidate();
return *this;
}
if (col_upb < col_lwb)
{
Error("Use","col_upb=%d < col_lwb=%d",col_upb,col_lwb);
Invalidate();
return *this;
}
Clear();
fNrows = row_upb-row_lwb+1;
fNcols = col_upb-col_lwb+1;
fRowLwb = row_lwb;
fColLwb = col_lwb;
fNelems = fNrows*fNcols;
fElements = data;
fIsOwner = kFALSE;
return *this;
}
//______________________________________________________________________________
TMatrixDBase &TMatrixD::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
TMatrixDBase &target,Option_t *option) const
{
// Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the
// returned matrix depends on the argument option:
//
// option == "S" : return [0..row_upb-row_lwb+1][0..col_upb-col_lwb+1] (default)
// else : return [row_lwb..row_upb][col_lwb..col_upb]
Assert(IsValid());
if (row_lwb < fRowLwb || row_lwb > fRowLwb+fNrows-1) {
Error("GetSub","row_lwb out of bounds");
target.Invalidate();
return target;
}
if (col_lwb < fColLwb || col_lwb > fColLwb+fNcols-1) {
Error("GetSub","col_lwb out of bounds");
target.Invalidate();
return target;
}
if (row_upb < fRowLwb || row_upb > fRowLwb+fNrows-1) {
Error("GetSub","row_upb out of bounds");
target.Invalidate();
return target;
}
if (col_upb < fColLwb || col_upb > fColLwb+fNcols-1) {
Error("GetSub","col_upb out of bounds");
target.Invalidate();
return target;
}
if (row_upb < row_lwb || col_upb < col_lwb) {
Error("GetSub","row_upb < row_lwb || col_upb < col_lwb");
target.Invalidate();
return target;
}
TString opt(option);
opt.ToUpper();
const Int_t shift = (opt.Contains("S")) ? 1 : 0;
const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
target.ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
if (target.GetRowIndexArray() && target.GetColIndexArray()) {
for (Int_t irow = 0; irow < nrows_sub; irow++) {
for (Int_t icol = 0; icol < ncols_sub; icol++) {
target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
}
}
} else {
const Double_t *ap = this->GetMatrixArray()+(row_lwb-fRowLwb)*fNcols+(col_lwb-fColLwb);
Double_t *bp = target.GetMatrixArray();
for (Int_t irow = 0; irow < nrows_sub; irow++) {
const Double_t *ap_sub = ap;
for (Int_t icol = 0; icol < ncols_sub; icol++) {
*bp++ = *ap_sub++;
}
ap += fNcols;
}
}
return target;
}
//______________________________________________________________________________
TMatrixDBase &TMatrixD::SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixDBase &source)
{
// Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part
// [row_lwb..row_lwb+nrows_source][col_lwb..col_lwb+ncols_source];
Assert(IsValid());
Assert(source.IsValid());
if (row_lwb < fRowLwb || row_lwb > fRowLwb+fNrows-1) {
Error("SetSub","row_lwb outof bounds");
Invalidate();
return *this;
}
if (col_lwb < fColLwb || col_lwb > fColLwb+fNcols-1) {
Error("SetSub","col_lwb outof bounds");
Invalidate();
return *this;
}
const Int_t nRows_source = source.GetNrows();
const Int_t nCols_source = source.GetNcols();
if (row_lwb+nRows_source > fRowLwb+fNrows || col_lwb+nCols_source > fColLwb+fNcols) {
Error("SetSub","source matrix too large");
Invalidate();
return *this;
}
if (source.GetRowIndexArray() && source.GetColIndexArray()) {
const Int_t rowlwb_s = source.GetRowLwb();
const Int_t collwb_s = source.GetColLwb();
for (Int_t irow = 0; irow < nRows_source; irow++) {
for (Int_t icol = 0; icol < nCols_source; icol++) {
(*this)(row_lwb+irow,col_lwb+icol) = source(rowlwb_s+irow,collwb_s+icol);
}
}
} else {
const Double_t *bp = source.GetMatrixArray();
Double_t *ap = this->GetMatrixArray()+(row_lwb-fRowLwb)*fNcols+(col_lwb-fColLwb);
for (Int_t irow = 0; irow < nRows_source; irow++) {
Double_t *ap_sub = ap;
for (Int_t icol = 0; icol < nCols_source; icol++) {
*ap_sub++ = *bp++;
}
ap += fNcols;
}
}
return *this;
}
//______________________________________________________________________________
TMatrixDBase &TMatrixD::ResizeTo(Int_t nrows,Int_t ncols,Int_t /*nr_nonzeros*/)
{
// Set size of the matrix to nrows x ncols
// New dynamic elements are created, the overlapping part of the old ones are
// copied to the new structures, then the old elements are deleted.
Assert(IsValid());
if (!fIsOwner) {
Error("ResizeTo(Int_t,Int_t)","Not owner of data array,cannot resize");
Invalidate();
return *this;
}
if (fNelems > 0) {
if (fNrows == nrows && fNcols == ncols)
return *this;
else if (nrows == 0 || ncols == 0) {
fNrows = nrows; fNcols = ncols;
Clear();
return *this;
}
Double_t *elements_old = GetMatrixArray();
const Int_t nelems_old = fNelems;
const Int_t nrows_old = fNrows;
const Int_t ncols_old = fNcols;
Allocate(nrows,ncols);
Assert(IsValid());
Double_t *elements_new = GetMatrixArray();
// new memory should be initialized but be careful ot to wipe out the stack
// storage. Initialize all when old or new storage was on the heap
if (fNelems > kSizeMax || nelems_old > kSizeMax)
memset(elements_new,0,fNelems*sizeof(Double_t));
else if (fNelems > nelems_old)
memset(elements_new+nelems_old,0,(fNelems-nelems_old)*sizeof(Double_t));
// Copy overlap
const Int_t ncols_copy = TMath::Min(fNcols,ncols_old);
const Int_t nrows_copy = TMath::Min(fNrows,nrows_old);
const Int_t nelems_new = fNelems;
if (ncols_old < fNcols) {
for (Int_t i = nrows_copy-1; i >= 0; i--)
Memcpy_m(elements_new+i*fNcols,elements_old+i*ncols_old,ncols_copy,
nelems_new,nelems_old);
} else {
for (Int_t i = 0; i < nrows_copy; i++)
Memcpy_m(elements_new+i*fNcols,elements_old+i*ncols_old,ncols_copy,
nelems_new,nelems_old);
}
Delete_m(nelems_old,elements_old);
} else {
Allocate(nrows,ncols,0,0,1);
}
return *this;
}
//______________________________________________________________________________
TMatrixDBase &TMatrixD::ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
Int_t /*nr_nonzeros*/)
{
// Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb]
// New dynamic elemenst are created, the overlapping part of the old ones are
// copied to the new structures, then the old elements are deleted.
Assert(IsValid());
if (!fIsOwner) {
Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
Invalidate();
return *this;
}
const Int_t new_nrows = row_upb-row_lwb+1;
const Int_t new_ncols = col_upb-col_lwb+1;
if (fNelems > 0) {
if (fNrows == new_nrows && fNcols == new_ncols &&
fRowLwb == row_lwb && fColLwb == col_lwb)
return *this;
else if (new_nrows == 0 || new_ncols == 0) {
fNrows = new_nrows; fNcols = new_ncols;
fRowLwb = row_lwb; fColLwb = col_lwb;
Clear();
return *this;
}
Double_t *elements_old = GetMatrixArray();
const Int_t nelems_old = fNelems;
const Int_t nrows_old = fNrows;
const Int_t ncols_old = fNcols;
const Int_t rowLwb_old = fRowLwb;
const Int_t colLwb_old = fColLwb;
Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
Assert(IsValid());
Double_t *elements_new = GetMatrixArray();
// new memory should be initialized but be careful ot to wipe out the stack
// storage. Initialize all when old or new storag ewas on the heap
if (fNelems > kSizeMax || nelems_old > kSizeMax)
memset(elements_new,0,fNelems*sizeof(Double_t));
else if (fNelems > nelems_old)
memset(elements_new+nelems_old,0,(fNelems-nelems_old)*sizeof(Double_t));
// Copy overlap
const Int_t rowLwb_copy = TMath::Max(fRowLwb,rowLwb_old);
const Int_t colLwb_copy = TMath::Max(fColLwb,colLwb_old);
const Int_t rowUpb_copy = TMath::Min(fRowLwb+fNrows-1,rowLwb_old+nrows_old-1);
const Int_t colUpb_copy = TMath::Min(fColLwb+fNcols-1,colLwb_old+ncols_old-1);
const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
if (nrows_copy > 0 && ncols_copy > 0) {
const Int_t colOldOff = colLwb_copy-colLwb_old;
const Int_t colNewOff = colLwb_copy-fColLwb;
if (ncols_old < fNcols) {
for (Int_t i = nrows_copy-1; i >= 0; i--) {
const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
const Int_t iRowNew = rowLwb_copy+i-fRowLwb;
Memcpy_m(elements_new+iRowNew*fNcols+colNewOff,
elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,fNelems,nelems_old);
}
} else {
for (Int_t i = 0; i < nrows_copy; i++) {
const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
const Int_t iRowNew = rowLwb_copy+i-fRowLwb;
Memcpy_m(elements_new+iRowNew*fNcols+colNewOff,
elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,fNelems,nelems_old);
}
}
}
Delete_m(nelems_old,elements_old);
} else {
Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
}
return *this;
}
//______________________________________________________________________________
Double_t TMatrixD::Determinant() const
{
const TMatrixD &tmp = *this;
TDecompLU lu(tmp,fTol);
Double_t d1,d2;
lu.Det(d1,d2);
return d1*TMath::Power(2.0,d2);
}
//______________________________________________________________________________
void TMatrixD::Determinant(Double_t &d1,Double_t &d2) const
{
const TMatrixD &tmp = *this;
TDecompLU lu(tmp,fTol);
lu.Det(d1,d2);
}
//______________________________________________________________________________
TMatrixD &TMatrixD::Invert(Double_t *det)
{
// Invert the matrix and calculate its determinant
Assert(IsValid());
TDecompLU::InvertLU(*this,fTol,det);
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::InvertFast(Double_t *det)
{
// Invert the matrix and calculate its determinant
Assert(IsValid());
const Char_t nRows = Char_t(GetNrows());
switch (nRows) {
case 1:
{
if (GetNrows() != GetNcols() || GetRowLwb() != GetColLwb()) {
Error("Invert()","matrix should be square");
Invalidate();
} else {
Double_t *pM = this->GetMatrixArray();
if (*pM == 0.) {
Error("InvertFast","matrix is singular");
Invalidate();
}
else
*pM = 1.0/(*pM);
}
return *this;
}
case 2:
{
TMatrixDCramerInv::Inv2x2(*this,det);
return *this;
}
case 3:
{
TMatrixDCramerInv::Inv3x3(*this,det);
return *this;
}
case 4:
{
TMatrixDCramerInv::Inv4x4(*this,det);
return *this;
}
case 5:
{
TMatrixDCramerInv::Inv5x5(*this,det);
return *this;
}
case 6:
{
TMatrixDCramerInv::Inv6x6(*this,det);
return *this;
}
default:
{
TDecompLU::InvertLU(*this,fTol,det);
return *this;
}
}
}
//______________________________________________________________________________
TMatrixD &TMatrixD::Transpose(const TMatrixD &source)
{
// Transpose a matrix.
Assert(IsValid());
Assert(source.IsValid());
if (this == &source) {
Double_t *ap = this->GetMatrixArray();
if (fNrows == fNcols && fRowLwb == fColLwb) {
for (Int_t i = 0; i < fNrows; i++) {
const Int_t off_i = i*fNrows;
for (Int_t j = i+1; j < fNcols; j++) {
const Int_t off_j = j*fNcols;
const Double_t tmp = ap[off_i+j];
ap[off_i+j] = ap[off_j+i];
ap[off_j+i] = tmp;
}
}
} else {
Double_t *oldElems = new Double_t[source.GetNoElements()];
memcpy(oldElems,source.GetMatrixArray(),source.GetNoElements()*sizeof(Double_t));
const Int_t nrows_old = fNrows;
const Int_t ncols_old = fNcols;
const Int_t rowlwb_old = fRowLwb;
const Int_t collwb_old = fColLwb;
fNrows = ncols_old; fNcols = nrows_old;
fRowLwb = collwb_old; fColLwb = rowlwb_old;
for (Int_t irow = fRowLwb; irow < fRowLwb+fNrows; irow++) {
for (Int_t icol = fColLwb; icol < fColLwb+fNcols; icol++) {
const Int_t off = (icol-collwb_old)*ncols_old;
(*this)(irow,icol) = oldElems[off+irow-rowlwb_old];
}
}
delete [] oldElems;
}
} else {
if (fNrows != source.GetNcols() || fNcols != source.GetNrows() ||
fRowLwb != source.GetColLwb() || fColLwb != source.GetRowLwb())
{
Error("Transpose","matrix has wrong shape");
Invalidate();
return *this;
}
const Double_t *sp1 = source.GetMatrixArray();
const Double_t *scp = sp1; // Row source pointer
Double_t *tp = this->GetMatrixArray();
const Double_t * const tp_last = this->GetMatrixArray()+fNelems;
// (This: target) matrix is traversed row-wise way,
// whilst the source matrix is scanned column-wise
while (tp < tp_last) {
const Double_t *sp2 = scp++;
// Move tp to the next elem in the row and sp to the next elem in the curr col
while (sp2 < sp1+fNelems) {
*tp++ = *sp2;
sp2 += fNrows;
}
}
Assert(tp == tp_last && scp == sp1+fNrows);
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::Rank1Update(const TVectorD &v,Double_t alpha)
{
// Perform a rank 1 operation on the matrix:
// A += alpha * v * v^T
Assert(IsValid());
Assert(v.IsValid());
if (v.GetNoElements() < TMath::Max(fNrows,fNcols)) {
Error("Rank1Update","vector too short");
Invalidate();
return *this;
}
const Double_t * const pv = v.GetMatrixArray();
Double_t *mp = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
const Double_t tmp = alpha*pv[i];
for (Int_t j = 0; j < fNcols; j++)
*mp++ += tmp*pv[j];
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::Rank1Update(const TVectorD &v1,const TVectorD &v2,Double_t alpha)
{
// Perform a rank 1 operation on the matrix:
// A += alpha * v1 * v2^T
Assert(IsValid());
Assert(v1.IsValid());
Assert(v2.IsValid());
if (v1.GetNoElements() < fNrows) {
Error("Rank1Update","vector v1 too short");
Invalidate();
return *this;
}
if (v2.GetNoElements() < fNcols) {
Error("Rank1Update","vector v2 too short");
Invalidate();
return *this;
}
const Double_t * const pv1 = v1.GetMatrixArray();
const Double_t * const pv2 = v2.GetMatrixArray();
Double_t *mp = this->GetMatrixArray();
for (Int_t i = 0; i < fNrows; i++) {
const Double_t tmp = alpha*pv1[i];
for (Int_t j = 0; j < fNcols; j++)
*mp++ += tmp*pv2[j];
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::NormByColumn(const TVectorD &v,Option_t *option)
{
// Multiply/divide matrix columns by a vector:
// option:
// "D" : b(i,j) = a(i,j)/v(i) i = 0,fNrows-1 (default)
// else : b(i,j) = a(i,j)*v(i)
Assert(IsValid());
Assert(v.IsValid());
if (v.GetNoElements() < fNrows) {
Error("NormByColumn","vector shorter than matrix column");
Invalidate();
return *this;
}
TString opt(option);
opt.ToUpper();
const Int_t divide = (opt.Contains("D")) ? 1 : 0;
const Double_t* pv = v.GetMatrixArray();
Double_t *mp = this->GetMatrixArray();
const Double_t * const mp_last = mp+fNelems;
if (divide) {
for ( ; mp < mp_last; pv++) {
for (Int_t j = 0; j < fNcols; j++)
{
Assert(*pv != 0.0);
*mp++ /= *pv;
}
}
} else {
for ( ; mp < mp_last; pv++)
for (Int_t j = 0; j < fNcols; j++)
*mp++ *= *pv;
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::NormByRow(const TVectorD &v,Option_t *option)
{
// Multiply/divide matrix rows with a vector:
// option:
// "D" : b(i,j) = a(i,j)/v(j) i = 0,fNcols-1 (default)
// else : b(i,j) = a(i,j)*v(j)
Assert(IsValid());
Assert(v.IsValid());
if (v.GetNoElements() < fNcols) {
Error("NormByRow","vector shorter than matrix column");
Invalidate();
return *this;
}
TString opt(option);
opt.ToUpper();
const Int_t divide = (opt.Contains("D")) ? 1 : 0;
const Double_t *pv0 = v.GetMatrixArray();
const Double_t *pv = pv0;
Double_t *mp = this->GetMatrixArray();
const Double_t * const mp_last = mp+fNelems;
if (divide) {
for ( ; mp < mp_last; pv = pv0 )
for (Int_t j = 0; j < fNcols; j++) {
Assert(*pv != 0.0);
*mp++ /= *pv++;
}
} else {
for ( ; mp < mp_last; pv = pv0 )
for (Int_t j = 0; j < fNcols; j++)
*mp++ *= *pv++;
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(const TMatrixD &source)
{
if (!AreCompatible(*this,source)) {
Error("operator=(const TMatrixD &)","matrices not compatible");
Invalidate();
return *this;
}
if (this != &source) {
TObject::operator=(source);
memcpy(fElements,source.GetMatrixArray(),fNelems*sizeof(Double_t));
fTol = source.GetTol();
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(const TMatrixF &source)
{
if (!AreCompatible(*this,source)) {
Error("operator=(const TMatrixF &)","matrices not compatible");
Invalidate();
return *this;
}
if (dynamic_cast<TMatrixF *>(this) != &source) {
TObject::operator=(source);
const Float_t * const ps = source.GetMatrixArray();
Double_t * const pt = GetMatrixArray();
for (Int_t i = 0; i < fNelems; i++)
pt[i] = (Double_t) ps[i];
fTol = (Double_t)source.GetTol();
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(const TMatrixDSym &source)
{
if (!AreCompatible(*this,source)) {
Error("operator=(const TMatrixDSym &)","matrices not compatible");
Invalidate();
return *this;
}
if ((TMatrixDBase *)this != (TMatrixDBase *)&source) {
TObject::operator=(source);
memcpy(fElements,source.GetMatrixArray(),fNelems*sizeof(Double_t));
fTol = source.GetTol();
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(const TMatrixDSparse &source)
{
if (GetNrows() != source.GetNrows() || GetNcols() != source.GetNcols() ||
GetRowLwb() != source.GetRowLwb() || GetColLwb() != source.GetColLwb()) {
Error("operator=(const TMatrixDSparse &","matrices not compatible");
Invalidate();
return *this;
}
if ((TMatrixDBase *)this != (TMatrixDBase *)&source) {
TObject::operator=(source);
memset(fElements,0,fNelems*sizeof(Double_t));
const Double_t * const sp = source.GetMatrixArray();
Double_t * tp = this->GetMatrixArray();
const Int_t * const pRowIndex = source.GetRowIndexArray();
const Int_t * const pColIndex = source.GetColIndexArray();
for (Int_t irow = 0; irow < fNrows; irow++ ) {
const Int_t off = irow*fNcols;
const Int_t sIndex = pRowIndex[irow];
const Int_t eIndex = pRowIndex[irow+1];
for (Int_t index = sIndex; index < eIndex; index++)
tp[off+pColIndex[index]] = sp[index];
}
fTol = source.GetTol();
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(const TMatrixDLazy &lazy_constructor)
{
Assert(IsValid());
if (lazy_constructor.GetRowUpb() != GetRowUpb() ||
lazy_constructor.GetColUpb() != GetColUpb() ||
lazy_constructor.GetRowLwb() != GetRowLwb() ||
lazy_constructor.GetColLwb() != GetColLwb()) {
Error("operator=(const TMatrixDLazy&)", "matrix is incompatible with "
"the assigned Lazy matrix");
Invalidate();
return *this;
}
lazy_constructor.FillIn(*this);
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator=(Double_t val)
{
// Assign val to every element of the matrix.
Assert(IsValid());
Double_t *ep = this->GetMatrixArray();
const Double_t * const ep_last = ep+fNelems;
while (ep < ep_last)
*ep++ = val;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator+=(Double_t val)
{
// Add val to every element of the matrix.
Assert(IsValid());
Double_t *ep = this->GetMatrixArray();
const Double_t * const ep_last = ep+fNelems;
while (ep < ep_last)
*ep++ += val;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator-=(Double_t val)
{
// Subtract val from every element of the matrix.
Assert(IsValid());
Double_t *ep = this->GetMatrixArray();
const Double_t * const ep_last = ep+fNelems;
while (ep < ep_last)
*ep++ -= val;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(Double_t val)
{
// Multiply every element of the matrix with val.
Assert(IsValid());
Double_t *ep = this->GetMatrixArray();
const Double_t * const ep_last = ep+fNelems;
while (ep < ep_last)
*ep++ *= val;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator+=(const TMatrixD &source)
{
// Add the source matrix.
if (!AreCompatible(*this,source)) {
Error("operator+=(const TMatrixD &)","matrices not compatible");
Invalidate();
return *this;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *tp = this->GetMatrixArray();
const Double_t * const tp_last = tp+fNelems;
while (tp < tp_last)
*tp++ += *sp++;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator+=(const TMatrixDSym &source)
{
// Add the source matrix.
if (!AreCompatible(*this,source)) {
Error("operator+=(const TMatrixDSym &)","matrices not compatible");
Invalidate();
return *this;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < fNrows; i++) {
sp += i;
trp += i; // point to [i,i]
tcp += i*fNcols; // point to [i,i]
for (Int_t j = i; j < fNcols; j++) {
if (j > i) *tcp += *sp;
*trp++ += *sp++;
tcp += fNcols;
}
tcp -= fNelems-1; // point to [0,i]
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator-=(const TMatrixD &source)
{
// Subtract the source matrix.
if (!AreCompatible(*this,source)) {
Error("operator=-(const TMatrixD &)","matrices not compatible");
Invalidate();
return *this;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *tp = this->GetMatrixArray();
const Double_t * const tp_last = tp+fNelems;
while (tp < tp_last)
*tp++ -= *sp++;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator-=(const TMatrixDSym &source)
{
// Subtract the source matrix.
if (!AreCompatible(*this,source)) {
Error("operator=-(const TMatrixDSym &)","matrices not compatible");
Invalidate();
return *this;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < fNrows; i++) {
sp += i;
trp += i; // point to [i,i]
tcp += i*fNcols; // point to [i,i]
for (Int_t j = i; j < fNcols; j++) {
if (j > i) *tcp -= *sp;
*trp++ -= *sp++;
tcp += fNcols;
}
tcp -= fNelems-1; // point to [0,i]
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(const TMatrixD &source)
{
// Compute target = target * source inplace. Strictly speaking, it can't be
// done inplace, though only the row of the target matrix needs to be saved.
// "Inplace" multiplication is only allowed when the 'source' matrix is square.
Assert(IsValid());
Assert(source.IsValid());
if (fNcols != source.GetNrows() || fColLwb != source.GetRowLwb() ||
fNcols != source.GetNcols() || fColLwb != source.GetColLwb()) {
Error("operator*=(const TMatrixD &)","source matrix has wrong shape");
Invalidate();
return *this;
}
// Check for A *= A;
const Double_t *sp;
TMatrixD tmp;
if (this == &source) {
tmp.ResizeTo(source);
tmp = source;
sp = tmp.GetMatrixArray();
}
else
sp = source.GetMatrixArray();
// One row of the old_target matrix
Double_t work[kWorkMax];
Bool_t isAllocated = kFALSE;
Double_t *trp = work;
if (fNcols > kWorkMax) {
isAllocated = kTRUE;
trp = new Double_t[fNcols];
}
Double_t *cp = this->GetMatrixArray();
const Double_t *trp0 = cp; // Pointer to target[i,0];
const Double_t * const trp0_last = trp0+fNelems;
while (trp0 < trp0_last) {
memcpy(trp,trp0,fNcols*sizeof(Double_t)); // copy the i-th row of target, Start at target[i,0]
for (const Double_t *scp = sp; scp < sp+fNcols; ) { // Pointer to the j-th column of source,
// Start scp = source[0,0]
Double_t cij = 0;
for (Int_t j = 0; j < fNcols; j++) {
cij += trp[j] * *scp; // the j-th col of source
scp += fNcols;
}
*cp++ = cij;
scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
}
trp0 += fNcols; // Set trp0 to the (i+1)-th row
Assert(trp0 == cp);
}
Assert(cp == trp0_last && trp0 == trp0_last);
if (isAllocated)
delete [] trp;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(const TMatrixDSym &source)
{
// Compute target = target * source inplace. Strictly speaking, it can't be
// done inplace, though only the row of the target matrix needs to be saved.
Assert(IsValid());
Assert(source.IsValid());
if (fNcols != source.GetNrows() || fColLwb != source.GetRowLwb()) {
Error("operator*=(const TMatrixDSym &)","source matrix has wrong shape");
Invalidate();
return *this;
}
// Check for A *= A;
const Double_t *sp;
TMatrixD tmp;
if ((TMatrixDBase *)this == (TMatrixDBase *)&source) {
tmp.ResizeTo(source);
tmp = source;
sp = tmp.GetMatrixArray();
}
else
sp = source.GetMatrixArray();
// One row of the old_target matrix
Double_t work[kWorkMax];
Bool_t isAllocated = kFALSE;
Double_t *trp = work;
if (fNcols > kWorkMax) {
isAllocated = kTRUE;
trp = new Double_t[fNcols];
}
Double_t *cp = this->GetMatrixArray();
const Double_t *trp0 = cp; // Pointer to target[i,0];
const Double_t * const trp0_last = trp0+fNelems;
while (trp0 < trp0_last) {
memcpy(trp,trp0,fNcols*sizeof(Double_t)); // copy the i-th row of target, Start at target[i,0]
for (const Double_t *scp = sp; scp < sp+fNcols; ) { // Pointer to the j-th column of source,
//Start scp = source[0,0]
Double_t cij = 0;
for (Int_t j = 0; j < fNcols; j++) {
cij += trp[j] * *scp; // the j-th col of source
scp += fNcols;
}
*cp++ = cij;
scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
}
trp0 += fNcols; // Set trp0 to the (i+1)-th row
Assert(trp0 == cp);
}
Assert(cp == trp0_last && trp0 == trp0_last);
if (isAllocated)
delete [] trp;
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(const TMatrixDDiag_const &diag)
{
// Multiply a matrix row by the diagonal of another matrix
// matrix(i,j) *= diag(j), j=1,fNcols
Assert(IsValid());
Assert(diag.GetMatrix()->IsValid());
Assert(fNcols == diag.GetNdiags());
if (fNcols != diag.GetNdiags()) {
Error("operator*=(const TMatrixDDiag_const &)","wrong diagonal length");
Invalidate();
return *this;
}
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Int_t inc = diag.GetInc();
while (mp < mp_last) {
const Double_t *dp = diag.GetPtr();
for (Int_t j = 0; j < fNcols; j++) {
*mp++ *= *dp;
dp += inc;
}
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator/=(const TMatrixDDiag_const &diag)
{
// Divide a matrix row by the diagonal of another matrix
// matrix(i,j) /= diag(j)
Assert(IsValid());
Assert(diag.GetMatrix()->IsValid());
if (fNcols != diag.GetNdiags()) {
Error("operator/=(const TMatrixDDiag_const &)","wrong diagonal length");
Invalidate();
return *this;
}
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Int_t inc = diag.GetInc();
while (mp < mp_last) {
const Double_t *dp = diag.GetPtr();
for (Int_t j = 0; j < fNcols; j++) {
Assert(*dp != 0.0);
*mp++ /= *dp;
dp += inc;
}
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(const TMatrixDColumn_const &col)
{
// Multiply a matrix by the column of another matrix
// matrix(i,j) *= another(i,k) for fixed k
const TMatrixDBase *mt = col.GetMatrix();
Assert(IsValid());
Assert(mt->IsValid());
if (fNrows != mt->GetNrows()) {
Error("operator*=(const TMatrixDColumn_const &)","wrong column length");
Invalidate();
return *this;
}
const Double_t * const endp = col.GetPtr()+mt->GetNoElements();
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Double_t *cp = col.GetPtr(); // ptr
const Int_t inc = col.GetInc();
while (mp < mp_last) {
Assert(cp < endp);
for (Int_t j = 0; j < fNcols; j++)
*mp++ *= *cp;
cp += inc;
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator/=(const TMatrixDColumn_const &col)
{
// Divide a matrix by the column of another matrix
// matrix(i,j) /= another(i,k) for fixed k
const TMatrixDBase *mt = col.GetMatrix();
Assert(IsValid());
Assert(mt->IsValid());
if (fNrows != mt->GetNrows()) {
Error("operator/=(const TMatrixDColumn_const &)","wrong column matrix");
Invalidate();
return *this;
}
const Double_t * const endp = col.GetPtr()+mt->GetNoElements();
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Double_t *cp = col.GetPtr(); // ptr
const Int_t inc = col.GetInc();
while (mp < mp_last) {
Assert(cp < endp);
Assert(*cp != 0.0);
for (Int_t j = 0; j < fNcols; j++)
*mp++ /= *cp;
cp += inc;
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator*=(const TMatrixDRow_const &row)
{
// Multiply a matrix by the row of another matrix
// matrix(i,j) *= another(k,j) for fixed k
const TMatrixDBase *mt = row.GetMatrix();
Assert(IsValid());
Assert(mt->IsValid());
if (fNcols != mt->GetNcols()) {
Error("operator*=(const TMatrixDRow_const &)","wrong row length");
Invalidate();
return *this;
}
const Double_t * const endp = row.GetPtr()+mt->GetNoElements();
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Int_t inc = row.GetInc();
while (mp < mp_last) {
const Double_t *rp = row.GetPtr(); // Row ptr
for (Int_t j = 0; j < fNcols; j++) {
Assert(rp < endp);
*mp++ *= *rp;
rp += inc;
}
}
return *this;
}
//______________________________________________________________________________
TMatrixD &TMatrixD::operator/=(const TMatrixDRow_const &row)
{
// Divide a matrix by the row of another matrix
// matrix(i,j) /= another(k,j) for fixed k
const TMatrixDBase *mt = row.GetMatrix();
Assert(IsValid());
Assert(mt->IsValid());
if (fNcols != mt->GetNcols()) {
Error("operator/=(const TMatrixDRow_const &)","wrong row length");
Invalidate();
return *this;
}
const Double_t * const endp = row.GetPtr()+mt->GetNoElements();
Double_t *mp = this->GetMatrixArray(); // Matrix ptr
const Double_t * const mp_last = mp+fNelems;
const Int_t inc = row.GetInc();
while (mp < mp_last) {
const Double_t *rp = row.GetPtr(); // Row ptr
for (Int_t j = 0; j < fNcols; j++) {
Assert(rp < endp);
Assert(*rp != 0.0);
*mp++ /= *rp;
rp += inc;
}
}
return *this;
}
//______________________________________________________________________________
const TMatrixD TMatrixD::EigenVectors(TVectorD &eigenValues) const
{
// Return a matrix containing the eigen-vectors ordered by descending values
// of Re^2+Im^2 of the complex eigen-values .
// If the matrix is asymmetric, only the real part of the eigen-values is
// returned . For full functionality use TMatrixDEigen .
if (!IsSymmetric())
Warning("EigenVectors(TVectorD &)","Only real part of eigen-values will be returned");
TMatrixDEigen eigen(*this);
eigenValues.ResizeTo(fNrows);
eigenValues = eigen.GetEigenValuesRe();
return eigen.GetEigenVectors();
}
//______________________________________________________________________________
TMatrixD operator+(const TMatrixD &source1,const TMatrixD &source2)
{
TMatrixD target(source1);
target += source2;
return target;
}
//______________________________________________________________________________
TMatrixD operator+(const TMatrixD &source1,const TMatrixDSym &source2)
{
TMatrixD target(source1);
target += source2;
return target;
}
//______________________________________________________________________________
TMatrixD operator+(const TMatrixDSym &source1,const TMatrixD &source2)
{
return operator+(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator+(const TMatrixD &source,Double_t val)
{
TMatrixD target(source);
target += val;
return target;
}
//______________________________________________________________________________
TMatrixD operator+(Double_t val,const TMatrixD &source)
{
return operator+(source,val);
}
//______________________________________________________________________________
TMatrixD operator-(const TMatrixD &source1,const TMatrixD &source2)
{
TMatrixD target(source1);
target -= source2;
return target;
}
//______________________________________________________________________________
TMatrixD operator-(const TMatrixD &source1,const TMatrixDSym &source2)
{
TMatrixD target(source1);
target -= source2;
return target;
}
//______________________________________________________________________________
TMatrixD operator-(const TMatrixDSym &source1,const TMatrixD &source2)
{
return -1.*(operator-(source2,source1));
}
//______________________________________________________________________________
TMatrixD operator-(const TMatrixD &source,Double_t val)
{
TMatrixD target(source);
target -= val;
return target;
}
//______________________________________________________________________________
TMatrixD operator-(Double_t val,const TMatrixD &source)
{
return -1.0*operator-(source,val);
}
//______________________________________________________________________________
TMatrixD operator*(Double_t val,const TMatrixD &source)
{
TMatrixD target(source);
target *= val;
return target;
}
//______________________________________________________________________________
TMatrixD operator*(const TMatrixD &source,Double_t val)
{
return operator*(val,source);
}
//______________________________________________________________________________
TMatrixD operator*(const TMatrixD &source1,const TMatrixD &source2)
{
TMatrixD target(source1,TMatrixD::kMult,source2);
return target;
}
//______________________________________________________________________________
TMatrixD operator*(const TMatrixD &source1,const TMatrixDSym &source2)
{
TMatrixD target(source1,TMatrixD::kMult,source2);
return target;
}
//______________________________________________________________________________
TMatrixD operator*(const TMatrixDSym &source1,const TMatrixD &source2)
{
TMatrixD target(source1,TMatrixD::kMult,source2);
return target;
}
//______________________________________________________________________________
TMatrixD operator*(const TMatrixDSym &source1,const TMatrixDSym &source2)
{
TMatrixD target(source1,TMatrixD::kMult,source2);
return target;
}
//______________________________________________________________________________
TMatrixD operator&&(const TMatrixD &source1,const TMatrixD &source2)
{
// Logical AND
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator&&(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last)
*tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
return target;
}
//______________________________________________________________________________
TMatrixD operator&&(const TMatrixD &source1,const TMatrixDSym &source2)
{
// Logical AND
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator&&(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j > i) *tcp = (*scp1 != 0.0) & (*sp2 != 0.0);
*trp++ = (*srp1++ != 0.0 && *sp2++ != 0.0);
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator&&(const TMatrixDSym &source1,const TMatrixD &source2)
{
// Logical AND
return operator&&(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator||(const TMatrixD &source1,const TMatrixD &source2)
{
// Logical OR
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator||(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last)
*tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
return target;
}
//______________________________________________________________________________
TMatrixD operator||(const TMatrixD &source1,const TMatrixDSym &source2)
{
// Logical OR
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator||(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j > i) *tcp = (*scp1 != 0.0) || (*sp2 != 0.0);
*trp++ = (*srp1++ != 0.0 || *sp2++ != 0.0);
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator||(const TMatrixDSym &source1,const TMatrixD &source2)
{
// Logical OR
return operator||(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator>(const TMatrixD &source1,const TMatrixD &source2)
{
// source1 > source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator|(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last) {
*tp++ = (*sp1) > (*sp2); sp1++; sp2++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator>(const TMatrixD &source1,const TMatrixDSym &source2)
{
// source1 > source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator>(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j > i) *tcp = (*scp1) > (*sp2);
*trp++ = (*srp1) > (*sp2); srp1++; sp2++;
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator>(const TMatrixDSym &source1,const TMatrixD &source2)
{
// source1 > source2
return operator<=(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator>=(const TMatrixD &source1,const TMatrixD &source2)
{
// source1 >= source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator>=(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last) {
*tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator>=(const TMatrixD &source1,const TMatrixDSym &source2)
{
// source1 >= source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator>=(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j > i) *tcp = (*scp1) >= (*sp2);
*trp++ = (*srp1) >= (*sp2); srp1++; sp2++;
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator>=(const TMatrixDSym &source1,const TMatrixD &source2)
{
// source1 >= source2
return operator<(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator<=(const TMatrixD &source1,const TMatrixD &source2)
{
// source1 <= source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator<=(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last) {
*tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator<=(const TMatrixD &source1,const TMatrixDSym &source2)
{
// source1 <= source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator<=(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j < i) *tcp = (*scp1) <= (*sp2);
*trp++ = (*srp1) <= (*sp2); srp1++; sp2++;
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator<=(const TMatrixDSym &source1,const TMatrixD &source2)
{
// source1 <= source2
return operator>(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator<(const TMatrixD &source1,const TMatrixD &source2)
{
// source1 < source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator<(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp < tp_last) {
*tp++ = (*sp1) < (*sp2); sp1++; sp2++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator<(const TMatrixD &source1,const TMatrixDSym &source2)
{
// source1 < source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator<(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j < i) *tcp = (*scp1) < (*sp2);
*trp++ = (*srp1) < (*sp2); srp1++; sp2++;
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator<(const TMatrixDSym &source1,const TMatrixD &source2)
{
// source1 < source2
return operator>=(source2,source1);
}
//______________________________________________________________________________
TMatrixD operator!=(const TMatrixD &source1,const TMatrixD &source2)
{
// source1 != source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator!=(const TMatrixD&,const TMatrixD&)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp1 = source1.GetMatrixArray();
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp != tp_last) {
*tp++ = (*sp1) != (*sp2); sp1++; sp2++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator!=(const TMatrixD &source1,const TMatrixDSym &source2)
{
// source1 != source2
TMatrixD target;
if (!AreCompatible(source1,source2)) {
Error("operator!=(const TMatrixD&,const TMatrixDSym&)","matrices not compatible");
target.Invalidate();
return target;
}
target.ResizeTo(source1);
const Int_t nelems = target.GetNoElements();
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Double_t *srp1 = source1.GetMatrixArray();
const Double_t *scp1 = srp1;
const Double_t *sp2 = source2.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = trp; // pointer to LL part, traverse col-wise
for (Int_t i = 0; i != nrows; i++) {
sp2 += i;
srp1 += i; trp += i; // point to [i,i]
scp1 += i*ncols; tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j != ncols; j++) {
if (j != i) *tcp = (*scp1) != (*sp2);
*trp++ = (*srp1) != (*sp2); srp1++; sp2++;
scp1 += ncols;
tcp += ncols;
}
scp1 -= nelems-1; // point to [0,i]
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD operator!=(const TMatrixDSym &source1,const TMatrixD &source2)
{
// source1 != source2
return operator!=(source2,source1);
}
/*
//______________________________________________________________________________
TMatrixD operator!=(const TMatrixD &source1,Double_t val)
{
// source1 != val
TMatrixD target; target.ResizeTo(source1);
const Double_t *sp = source1.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t * const tp_last = tp+target.GetNoElements();
while (tp != tp_last) {
*tp++ = (*sp != val); sp++;
}
return target;
}
//______________________________________________________________________________
TMatrixD operator!=(Double_t val,const TMatrixD &source1)
{
// val != source1
return operator!=(source1,val);
}
*/
//______________________________________________________________________________
TMatrixD &Add(TMatrixD &target,Double_t scalar,const TMatrixD &source)
{
// Modify addition: target += scalar * source.
if (!AreCompatible(target,source)) {
::Error("Add(TMatrixD &,Double_t,const TMatrixD &)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t *ftp = tp+target.GetNoElements();
while ( tp < ftp )
*tp++ += scalar * (*sp++);
return target;
}
//______________________________________________________________________________
TMatrixD &Add(TMatrixD &target,Double_t scalar,const TMatrixDSym &source)
{
// Modify addition: target += scalar * source.
if (!AreCompatible(target,source)) {
::Error("Add(TMatrixD &,Double_t,const TMatrixDSym &)","matrices not compatible");
target.Invalidate();
return target;
}
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Int_t nelems = target.GetNoElements();
const Double_t *sp = source.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp += i;
trp += i; // point to [i,i]
tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
const Double_t tmp = scalar * *sp++;
if (j > i) *tcp += tmp;
tcp += ncols;
*trp++ += tmp;
}
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD &ElementMult(TMatrixD &target,const TMatrixD &source)
{
// Multiply target by the source, element-by-element.
if (!AreCompatible(target,source)) {
::Error("ElementMult(TMatrixD &,const TMatrixD &)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t *ftp = tp+target.GetNoElements();
while ( tp < ftp )
*tp++ *= *sp++;
return target;
}
//______________________________________________________________________________
TMatrixD &ElementMult(TMatrixD &target,const TMatrixDSym &source)
{
// Multiply target by the source, element-by-element.
if (!AreCompatible(target,source)) {
::Error("ElementMult(TMatrixD &,const TMatrixDSym &)","matrices not compatible");
target.Invalidate();
return target;
}
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Int_t nelems = target.GetNoElements();
const Double_t *sp = source.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp += i;
trp += i; // point to [i,i]
tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
if (j > i) *tcp *= *sp;
*trp++ *= *sp++;
tcp += ncols;
}
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
TMatrixD &ElementDiv(TMatrixD &target,const TMatrixD &source)
{
// Divide target by the source, element-by-element.
if (!AreCompatible(target,source)) {
::Error("ElementDiv(TMatrixD &,const TMatrixD &)","matrices not compatible");
target.Invalidate();
return target;
}
const Double_t *sp = source.GetMatrixArray();
Double_t *tp = target.GetMatrixArray();
const Double_t *ftp = tp+target.GetNoElements();
while ( tp < ftp ) {
Assert(*sp != 0.0);
*tp++ /= *sp++;
}
return target;
}
//______________________________________________________________________________
TMatrixD &ElementDiv(TMatrixD &target,const TMatrixDSym &source)
{
// Multiply target by the source, element-by-element.
if (!AreCompatible(target,source)) {
::Error("ElementDiv(TMatrixD &,const TMatrixDSym &)","matrices not compatible");
target.Invalidate();
return target;
}
const Int_t nrows = target.GetNrows();
const Int_t ncols = target.GetNcols();
const Int_t nelems = target.GetNoElements();
const Double_t *sp = source.GetMatrixArray();
Double_t *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
Double_t *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
for (Int_t i = 0; i < nrows; i++) {
sp += i;
trp += i; // point to [i,i]
tcp += i*ncols; // point to [i,i]
for (Int_t j = i; j < ncols; j++) {
Assert(*sp != 0.0);
if (j > i) *tcp /= *sp;
*trp++ /= *sp++;
tcp += ncols;
}
tcp -= nelems-1; // point to [0,i]
}
return target;
}
//______________________________________________________________________________
void TMatrixD::Streamer(TBuffer &R__b)
{
// Stream an object of class TMatrixD.
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v > 2) {
Clear();
TMatrixD::Class()->ReadBuffer(R__b,this,R__v,R__s,R__c);
} else if (R__v == 2) { //process old version 2
Clear();
TObject::Streamer(R__b);
MakeValid();
R__b >> fNrows;
R__b >> fNcols;
R__b >> fNelems;
R__b >> fRowLwb;
R__b >> fColLwb;
Char_t isArray;
R__b >> isArray;
if (isArray) {
if (fNelems > 0) {
fElements = new Double_t[fNelems];
R__b.ReadFastArray(fElements,fNelems);
} else
fElements = 0;
}
R__b.CheckByteCount(R__s,R__c,TMatrixD::IsA());
} else { //====process old versions before automatic schema evolution
TObject::Streamer(R__b);
MakeValid();
R__b >> fNrows;
R__b >> fNcols;
R__b >> fRowLwb;
R__b >> fColLwb;
fNelems = R__b.ReadArray(fElements);
R__b.CheckByteCount(R__s,R__c,TMatrixD::IsA());
}
// in version <=2 , the matrix was stored column-wise
if (R__v <= 2) {
for (Int_t i = 0; i < fNrows; i++) {
const Int_t off_i = i*fNcols;
for (Int_t j = i; j < fNcols; j++) {
const Int_t off_j = j*fNrows;
const Double_t tmp = fElements[off_i+j];
fElements[off_i+j] = fElements[off_j+i];
fElements[off_j+i] = tmp;
}
}
}
if (fNelems > 0 && fNelems <= kSizeMax) {
memcpy(fDataStack,fElements,fNelems*sizeof(Double_t));
delete [] fElements;
fElements = fDataStack;
} else if (fNelems < 0)
Invalidate();
} else {
TMatrixD::Class()->WriteBuffer(R__b,this);
}
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.