首页 > Uncategorized > Computing Determinant in Fortran

Computing Determinant in Fortran

2012年01月6日 留下评论 Go to comments

It may seem that one should go to LAPACK and find a routine that does just this job. However, it turns out that LAPACK does not really have such a routine. So a little bit of detour is needed. If you want to deal with a general, complex matrix, not symmetric, not Hermitian, then the way to go is using Gauss elimination(GE). LAPACK provides routines to find the LU decomposition of a matrix, which is essentially GE and can be used to compute the determinant. Here is the code snippet:


 1 complex(kind=8) function det(N, mat)
 2     implicit none
 3     integer(kind=8), intent(in) :: N
 4     complex(kind=8), intent(inout), dimension(:,:) :: mat
 5     integer(kind=8) :: i, info
 6     integer, allocatable :: ipiv(:)
 7     real(kind=8) :: sgn
 8
 9     allocate(ipiv(N))
10
11     ipiv = 0
12     
13     call zgetrf(N, N, mat, N, ipiv, info)
14     
15     det = ONE
16     do i = 1, N
17         det = det*mat(i, i)
18     end do
19     sgn = ONE
20     do i = 1, N
21         if(ipiv(i) /= i) then
22             sgn = -sgn
23         end if
24     end do
25     det = sgn*det  
26 end function det
Advertisements
分类:Uncategorized
  1. 2013年07月23日 @ 3:19 上午

    This is very useful. Thanks.
    BTW, how do you highlight your code?

    • 2013年07月23日 @ 5:31 上午

      Vim can turn the highlighted code into html

  1. No trackbacks yet.

发表评论

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 更改 )

Twitter picture

You are commenting using your Twitter account. Log Out / 更改 )

Facebook photo

You are commenting using your Facebook account. Log Out / 更改 )

Google+ photo

You are commenting using your Google+ account. Log Out / 更改 )

Connecting to %s

%d 博主赞过: