Page 1 of 1

Matrix inversion with clapack

PostPosted: Tue Sep 25, 2007 2:42 pm
by WolframW
Presumably I am the 100th person with that mistake, but I unfortunately haven't found the answer when skimming through this forum, so sorry in case I repeat an answered question.

I try to invert a general matrix in clapack. The code I use is the following one:

===
(...)
int *ipiv = new int[dimensionalityX];
for (int i = 0; i < dimensionalityX * dimensionalityX; i++) inverseHessian[i] = hessian[i];
std::cout << "return code of clapack_cgetrf: " << clapack_cgetrf (CblasRowMajor, dimensionalityX, dimensionalityX, inverseHessian, dimensionalityX, ipiv) << std::endl;
std::cout << "return code of clapack_dgetri: " << clapack_dgetri (CblasRowMajor, dimensionalityX, inverseHessian, dimensionalityX, ipiv) << std::endl;
(...)
===

(In fact, we should use Cholesky for this but I would like to start with this simple case.) The array hessian has the following content:

hessian = (118.31, -114.225, -114.225, 111.939),

and dimensionalityX = 2, the number of rows and columns in hessian.

The result should be (according to an external calculation):

inverseHessian = (0.57067342997703, 0.58232762968337, 0.58232762968337, 0.60315326651643),

but oddly, the result of my code is:

inverseHessian = (0.00845236, -5.27401e-19, 2.85921e+19, 2.96148e+19).

Does anybody have an explanation for this? I would be very grateful for the answer.

Best,

Wolfram W.

PostPosted: Tue Sep 25, 2007 7:23 pm
by Julien Langou
Hello,
1- where did you get your clapack from. This is not the interface for LAPACK fortran routine.
The interface you are describing is a C-frienly interface.
2- you should definetly try to call dgetrf and not cgetrf for the first call.
Julien.

PostPosted: Wed Sep 26, 2007 4:18 am
by WolframW
Arghh, I'm really sorry for that. I did not see the mistake - now everything works fine.