Clicky

Fortran Wiki
Matrix inversion (changes)

Skip the Navigation Links | Home Page | All Pages | Recently Revised | Authors | Feeds | Export |

Showing changes from revision #1 to #2: (追記) Added (追記ここまで) | (削除) Removed (削除ここまで) | (削除) Chan (削除ここまで)(追記) ged (追記ここまで)

(追記) (追記ここまで)(追記) (追記ここまで)(追記)

Using LAPACK (large matrices)

(追記ここまで)
(追記) (追記ここまで)
! Returns the inverse of a matrix calculated by finding the LU
! decomposition. Depends on LAPACK.
function inv(A) result(Ainv)
 real(dp), dimension(:,:), intent(in) :: A
 real(dp), dimension(size(A,1),size(A,2)) :: Ainv
 real(dp), dimension(size(A,1)) :: work ! work array for LAPACK
 integer, dimension(size(A,1)) :: ipiv ! pivot indices
 integer :: n, info
 ! External procedures defined in LAPACK
 external DGETRF
 external DGETRI
 ! Store A in Ainv to prevent it from being overwritten by LAPACK
 Ainv = A
 n = size(A,1)
 ! DGETRF computes an LU factorization of a general M-by-N matrix A
 ! using partial pivoting with row interchanges.
 call DGETRF(n, n, Ainv, n, ipiv, info)
 if (info /= 0) then
 stop 'Matrix is numerically singular!'
 end if
 ! DGETRI computes the inverse of a matrix using the LU factorization
 ! computed by DGETRF.
 call DGETRI(n, Ainv, n, ipiv, work, n, info)
 if (info /= 0) then
 stop 'Matrix inversion failed!'
 end if
end function inv
(追記) (追記ここまで)(追記)

Direct calculation (small matrices)

(追記ここまで)
(追記) (追記ここまで)(追記)

In my experience, LAPACK is great when you wish to invert huge ×ばつN matrices, but it can be really slow for inverting smaller ×ばつ2, ×ばつ3, and ×ばつ4 matrices. For my use case, where I need to invert billions of ×ばつ2 and ×ばつ4 matrices instead of a few large ×ばつN matrices, I got a 30% speedup of my program replacing the LAPACK calls by direct calculations of the matrix inversions. I have attached the code that I’ve used for the ×ばつ2, ×ばつ3, and ×ばつ4 cases below. The ×ばつ2 version is quite easy to derive analytically. The ×ばつ3 and ×ばつ4 versions are based on the subroutines M33INV and M44INV by David G. Simpson; I just converted them from subroutines to pure functions.

(追記ここまで)
(追記) (追記ここまで)(追記)
 pure function matinv2(A) result(B)
 !! Performs a direct calculation of the inverse of a ×ばつ2 matrix.
 complex(wp), intent(in) :: A(2,2) !! Matrix
 complex(wp) :: B(2,2) !! Inverse matrix
 complex(wp) :: detinv
 ! Calculate the inverse determinant of the matrix
 detinv = 1/(A(1,1)*A(2,2) - A(1,2)*A(2,1))
 ! Calculate the inverse of the matrix
 B(1,1) = +detinv * A(2,2)
 B(2,1) = -detinv * A(2,1)
 B(1,2) = -detinv * A(1,2)
 B(2,2) = +detinv * A(1,1)
 end function
 pure function matinv3(A) result(B)
 !! Performs a direct calculation of the inverse of a ×ばつ3 matrix.
 complex(wp), intent(in) :: A(3,3) !! Matrix
 complex(wp) :: B(3,3) !! Inverse matrix
 complex(wp) :: detinv
 ! Calculate the inverse determinant of the matrix
 detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
 - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
 + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))
 ! Calculate the inverse of the matrix
 B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
 B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
 B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
 B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
 B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
 B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
 B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
 B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
 B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1))
 end function
 pure function matinv4(A) result(B)
 !! Performs a direct calculation of the inverse of a ×ばつ4 matrix.
 complex(wp), intent(in) :: A(4,4) !! Matrix
 complex(wp) :: B(4,4) !! Inverse matrix
 complex(wp) :: detinv
 ! Calculate the inverse determinant of the matrix
 detinv = &
 1/(A(1,1)*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2)))&
 - A(1,2)*(A(2,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1)))&
 + A(1,3)*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))&
 - A(1,4)*(A(2,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(2,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(2,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))))
 ! Calculate the inverse of the matrix
 B(1,1) = detinv*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2)))
 B(2,1) = detinv*(A(2,1)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(2,3)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(2,4)*(A(3,3)*A(4,1)-A(3,1)*A(4,3)))
 B(3,1) = detinv*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))
 B(4,1) = detinv*(A(2,1)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))+A(2,2)*(A(3,1)*A(4,3)-A(3,3)*A(4,1))+A(2,3)*(A(3,2)*A(4,1)-A(3,1)*A(4,2)))
 B(1,2) = detinv*(A(1,2)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(1,3)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(1,4)*(A(3,3)*A(4,2)-A(3,2)*A(4,3)))
 B(2,2) = detinv*(A(1,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(1,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(1,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1)))
 B(3,2) = detinv*(A(1,1)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(1,2)*(A(3,1)*A(4,4)-A(3,4)*A(4,1))+A(1,4)*(A(3,2)*A(4,1)-A(3,1)*A(4,2)))
 B(4,2) = detinv*(A(1,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(1,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(1,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))
 B(1,3) = detinv*(A(1,2)*(A(2,3)*A(4,4)-A(2,4)*A(4,3))+A(1,3)*(A(2,4)*A(4,2)-A(2,2)*A(4,4))+A(1,4)*(A(2,2)*A(4,3)-A(2,3)*A(4,2)))
 B(2,3) = detinv*(A(1,1)*(A(2,4)*A(4,3)-A(2,3)*A(4,4))+A(1,3)*(A(2,1)*A(4,4)-A(2,4)*A(4,1))+A(1,4)*(A(2,3)*A(4,1)-A(2,1)*A(4,3)))
 B(3,3) = detinv*(A(1,1)*(A(2,2)*A(4,4)-A(2,4)*A(4,2))+A(1,2)*(A(2,4)*A(4,1)-A(2,1)*A(4,4))+A(1,4)*(A(2,1)*A(4,2)-A(2,2)*A(4,1)))
 B(4,3) = detinv*(A(1,1)*(A(2,3)*A(4,2)-A(2,2)*A(4,3))+A(1,2)*(A(2,1)*A(4,3)-A(2,3)*A(4,1))+A(1,3)*(A(2,2)*A(4,1)-A(2,1)*A(4,2)))
 B(1,4) = detinv*(A(1,2)*(A(2,4)*A(3,3)-A(2,3)*A(3,4))+A(1,3)*(A(2,2)*A(3,4)-A(2,4)*A(3,2))+A(1,4)*(A(2,3)*A(3,2)-A(2,2)*A(3,3)))
 B(2,4) = detinv*(A(1,1)*(A(2,3)*A(3,4)-A(2,4)*A(3,3))+A(1,3)*(A(2,4)*A(3,1)-A(2,1)*A(3,4))+A(1,4)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)))
 B(3,4) = detinv*(A(1,1)*(A(2,4)*A(3,2)-A(2,2)*A(3,4))+A(1,2)*(A(2,1)*A(3,4)-A(2,4)*A(3,1))+A(1,4)*(A(2,2)*A(3,1)-A(2,1)*A(3,2)))
 B(4,4) = detinv*(A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)))
 end function
(追記ここまで)

category: code

Revised on April 22, 2016 16:57:10 by jabirali (46.9.153.214) (6444 characters / 2.0 pages)
Edit | Back in time (1 revision) | Hide changes | History | Views: Print | TeX | Source

AltStyle によって変換されたページ (->オリジナル) /