Hi,
I’m trying to use TDecompSparse to solve Ax=b where A is a matrix of the form below which I’m using a TMatrixDSparse to model.
1 0 0 0 0
1 -2 1 0 0
0 1 -2 1 0
0 0 1 -2 1
0 0 0 0 1
and b equals to:
1
0
0
0
10
Of course these are just a prototype matrix and vector that in the future are meant to have very higher number of rows and possibly two or four extra diagonals. Hence, the use of Sparse functionalities. But for some reason I get the wrong answer when I get TDecompSparse to solve the equation. Attached is a part of my code and its output.
TMatrixDSparse A(0,nmax-1,0,nmax-1,nel,row,col,dat);
Int_t *rIndex = A.GetRowIndexArray();
Int_t *cIndex = A.GetColIndexArray();
Double_t *pData = A.GetMatrixArray();
for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
Int_t sIndex = rIndex[irow];
Int_t eIndex = rIndex[irow+1];
for (Int_t index = sIndex; index < eIndex; index++) {
Int_t icol = cIndex[index];
Double_t data = pData[index];
printf("data(%d,%d) = %.4e\n",irow+A.GetRowLwb(),icol+A.GetColLwb(),data);
}
}
TVectorD b(0,nmax-1);
TVectorD r(0,nmax-1);
b(0)=0.0;
b(nmax-1)=10.0;
for(int i=0; i<nmax; i++)cout <<"b("<< i << "): " << b(i) << endl;
//TDecompSVD ls(A);
TDecompSparse ls(A,0);
bool fine = ls.Decompose();
r = ls.Solve(b,fine);
for(int i=0; i<nmax; i++)cout <<"r("<< i << "): " << r(i) << endl;
The output I get is:
data(0,0) = 1.0000e+00
data(1,0) = 1.0000e+00
data(1,1) = -2.0000e+00
data(1,2) = 1.0000e+00
data(2,1) = 1.0000e+00
data(2,2) = -2.0000e+00
data(2,3) = 1.0000e+00
data(3,2) = 1.0000e+00
data(3,3) = -2.0000e+00
data(3,4) = 1.0000e+00
data(4,4) = 1.0000e+00
b(0): 0
b(1): 0
b(2): 0
b(3): 0
b(4): 10
r(0): 0
r(1): 0
r(2): 0
r(3): 0
r(4): 10
Where r vector gives me a clear incorrect answer. Using TDecompSVD it works just fine for this case, but I cannot use it for a high enough number of rows. So… does anybody can help me how to get TDecompSparse to work properly?
Best regards,
fmarinho