# TDecompSparse failing to solve

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

the matrix has to be symmetric:

///////////////////////////////////////////////////////////////////////////
// //
// Sparse Symmetric Decomposition class //
// //
// Solve a sparse symmetric system of linear equations using a method //
// based on Gaussian elimination as discussed in Duff and Reid, //
// ACM Trans. Math. Software 9 (1983), 302-325. //
// //
///////////////////////////////////////////////////////////////////////////

Thanks! Just thought of having it all in ROOT code. Found myself a dedicated package as a solution (Eigen).

[quote=“Eddy Offermann”]the matrix has to be symmetric:

///////////////////////////////////////////////////////////////////////////
// //
// Sparse Symmetric Decomposition class //
// //
// Solve a sparse symmetric system of linear equations using a method //
// based on Gaussian elimination as discussed in Duff and Reid, //
// ACM Trans. Math. Software 9 (1983), 302-325. //
// //
///////////////////////////////////////////////////////////////////////////[/quote]