The LAPACK forum has moved to https://github.com/Reference-LAPACK/lapack/discussions.

Matrix inversion with clapack

Open discussion regarding features, bugs, issues, vendors, etc.

Matrix inversion with clapack

Postby WolframW » Tue Sep 25, 2007 2:42 pm

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.
WolframW
 
Posts: 2
Joined: Tue Sep 25, 2007 2:32 pm

Postby Julien Langou » Tue Sep 25, 2007 7:23 pm

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.
Julien Langou
 
Posts: 835
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby WolframW » Wed Sep 26, 2007 4:18 am

Arghh, I'm really sorry for that. I did not see the mistake - now everything works fine.
WolframW
 
Posts: 2
Joined: Tue Sep 25, 2007 2:32 pm


Return to User Discussion

Who is online

Users browsing this forum: No registered users and 4 guests