determinant by LU-decomposition with ?getrf

determinant by LU-decomposition with ?getrf

Bild des Benutzers dierksen
dierksen
Mi, 01/12/2004 - 17:54

Hi,

I want to calculate the determinant of a matrix by LU-decomposition with ?getrf. Since ?getrf does not compute the decomposition of the original matrix but of a matrix obtained by row permutations of the original one, the sign of the determinant depends on the number of those permutations. The MKL reference tells me that this information is contained in the array 'ipiv' where the pivot indices mean: row i was interchanged with row ipiv(i). So when ipiv(i)==j with i/=j this means that the rows i and j have been exchanged. As far as I understand this also means that ipiv(j)==i. However I checked this and this is not the case. So what's the meaning of this 'ipiv' array?

Best regards,
Marc

RSS Oben
5 Beiträge / 0 neu
Letzter Beitrag
Nähere Informationen zur Compiler-Optimierung finden Sie in unserem Optimierungshinweis.
Bild des Benutzers Community Admin
  • Zum Verfassen von Kommentaren bitte anmelden.
Community Admin
Do, 02/12/2004 - 19:32

A common error is to think ob ipiv as a permutation vector - it isn't. Rather it is a sequential record of the swaps that have taken place during factoization. An acceptable ipiv vector could have all the vales the same, namely the number of eqautions. o the first row is swapped with the last. Next the second rpow is swapped with the last. then the third and so on. Further, ipiv(i)/=j, j

Oben
Bild des Benutzers dierksen
  • Zum Verfassen von Kommentaren bitte anmelden.
dierksen
Fr, 03/12/2004 - 18:26

Thank you very much for your answer. It seems to be that your message got tructated somehow, but I think I got the point.
I also did some searching on the web and found the following two pieces of code which both claim to calculate the determinant of a square matrix.

============================================

call dgetrf(n,n,a,n,piv,info)
det = 0d0
if (info.ne.0) then
return
endif
det = 1d0
do 10,i=1,n
if (piv(i).ne.i) then
det = -det * a(i,i)
else
det = det * a(i,i)
endif
10 continue
end

============================================

call dgetrf(n, n, a, n, ipiv, info)
if(info < 0) then
write(message(1),'(a,i4)') 'In states_mpdotp, LAPACK dgetrf returned error code ',info
call write_fatal(1)
endif
ddet = M_ONE
do i = 1, n
ddet = ddet*a(i, i)
if(ipiv(i).ne.i) then
ipiv(ipiv(i)) = ipiv(i)
ddet = -ddet
endif
enddo

============================================

So according to your explanation above the author of the first code did it right, while the other one falsely assumed that ipiv corresponds to a permutation vector as I did. Do I get this right?

Oben
Bild des Benutzers Community Admin
  • Zum Verfassen von Kommentaren bitte anmelden.
Community Admin
Sa, 04/12/2004 - 20:08

You are correct. The second case treats ipiv as a permutation vector. I have not verified that the first loop is correct, either, but at least it is treating the ipiv vector as a pivot vector and not a permutation vector. I don't think the second case has the effect the user desires. Assume that for N = 3 the first two rows are exchanged with the last row. The pivot vector is (3,3,3). ipiv(ipiv(i)) = ipiv(i) for ipiv(i)/=i does not change the pivot vector. For instance, for i = 1, ipiv(1) /= 1 so the code says make ipiv(ipiv(1)) = ipiv(1). So ipiv(3) = ipiv(1) = 3, which is unchanged. The only reason the pivot vector needs to come into play for determining the determinant is to get the signs right. One of the things I had typed but which had gotten truncated in the previous message was that ipiv(i) >= i for all i since you never pivot with a row prior to row you are currently working on.

Oben
Bild des Benutzers dierksen
  • Zum Verfassen von Kommentaren bitte anmelden.
dierksen
So, 05/12/2004 - 19:31

The first one seems to be correct since the total number of permutations P is equal to the number of elements of ipiv with ipiv(i)/=i and the determinant D is changed to D*(-1)**P. I checked this for a couple of random matrices with Mathematica and it seems to work. Thanks again for your help.

Best wishes,
Marc

Oben

Melden Sie sich an, um einen Kommentar zu hinterlassen.