I need to apply the following matrix calculations:
trp(M) * M //multiply transposed matrix with matrix
vector * M //multiply matrix with vector
const * M //multiply matrix with constant
diagonale(M) //get the matrix diagonale
M1 * M2 //multiply 2 matrices
M1 + M2 //add 2 matrices
The matrices that I need to handle are Double_t and can have sizes of up to 500 x 600 (or even larger),
and I want to do the calculation on a computer with only 1GB RAM.
My question is:
Can I use class TMatrixT for this purpose or does ROOT supply a better alternative?
The reason for my question is that I want to convert an algorithm implemented in R using matrix notation to C++.
Meanwhile I am able to aks my question in a more specific way, since I have converted most of the R code using matrices in the usual way.
Since I need to pass the original matrix X as vector X[i*ncol+j] I have implemented e.g. “mm = tr(X)*X” and “XX = (mm + tr(mm))/2” as follows:
Double_t *mm = new Double_t[nrow*nrow];
Double_t *XX = new Double_t[nrow*nrow];
// mm = tr(X)*X/ncol
Double_t sum = 0.0;
for (Int_t i=0; i<nrow; i++) {
for (Int_t k=0; k<nrow; k++) {
sum = 0.0;
for (Int_t j=0; j<ncol; j++) {
sum += X[i*ncol + j]*X[k*ncol + j];
}//for_j
mm[i*nrow + k] = sum;
}//for_k
}//for_i
// XX = (mm + tr(mm))/2 but positive definit
for (Int_t i=0; i<nrow; i++) {
for (Int_t k=0; k<nrow; k++) {
sum = (mm[i*nrow + k] + mm[k*nrow + i])/2.0;
XX[i*nrow + k] = (sum > 0.0) ? sum : 0.0;
}//for_k
}//for_i
delete [] XX;
delete [] mm;
Thus my questions are:
What is the advantage of using TMatrixT or TMatrixD?
Is using TMatrixT faster than my current implementation?
Does using TMatrixT need less RAM than my current implementation?
TMatrixD is just a “Double_t” representation of the templated TMatrixT …
I see no advantage of using TMatrixT insteade of your code as far as speed
or memory usage is concerned .
However, I think using TMatrixT results in compacter code and less
errors because the algorithms are already contained in this class .
It also allows for easy futher manipulation of the matrix .
For instance “mm = tr(X)*X” would be:
TMatrixD *X = new TMatrixD(nrow,nrow);
TMatrixD mm(TMatrixD::kAtA,*X);
and XX = (mm + tr(mm))/2":
TMatrixD XX = 0.5*(mm+TMatrixD(TMatrixD::kTransposed,mm);
trp(M) * M //multiply transposed matrix with matrix: see above
vector * M //multiply matrix with vector : M * vector
const * M //multiply matrix with constant: see above
diagonale(M) //get the matrix diagonale : TVectorD v = TMatrixDDiag(M)
M1 * M2 //multiply 2 matrices : TMatrixD(M1,TMatrixD::kMult,M2) or M1*m2
M1 + M2 //add 2 matrices : TMatrixD(M1,TMatrixD::kPlus,M2) or M1 + M2
Thank you for your reply and especially for the demonstration code on how to use TMatrixD.
I agree with you that using TMatrixT results in compacter code and is less error prone, and
I intend to use it in the future, too.