-
Notifications
You must be signed in to change notification settings - Fork 0
/
inverse.f90
executable file
·113 lines (76 loc) · 3.07 KB
/
inverse.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
subroutine inv(ndim,Amat)
implicit none
integer,parameter :: dp=8
integer :: i
integer :: info
integer,intent(in):: ndim
! IPIV : INTEGER. Array, DIMENSION at least max(1, n). The pivot indices that define
! the permutation matrix P; row i of the matrix was interchanged with
! row ipiv(i). Corresponds to the single precision factorization (if info=
! 0 and iter ≥ 0) or the double precision factorization (if info= 0 and
! iter < 0).
integer,allocatable :: ipiv(:)
complex(dp),parameter :: zone=(1.0d0,0.0d0)
! Amat :
! Overwritten by the factors L and U from the factorization of A = P*L*U;
! the unit diagonal elements of L are not stored.
! If iterative refinement has been successfully used (info= 0 and iter≥
! 0), then A is unchanged.
! If double precision factorization has been used (info= 0 and iter <
! 0), then the array A contains the factors L and U from the factorization
! A = P*L*U; the unit diagonal elements of L are not stored.
complex(dp),intent(inout):: Amat(ndim,ndim)
! Bmat :
! Overwritten by the solution matrix X for dgesv, sgesv,zgesv,zgesv.
complex(dp),allocatable :: Bmat(:,:)
allocate(ipiv(ndim))
allocate(Bmat(ndim,ndim))
ipiv=0
! unit matrix
Bmat= (0d0, 0d0)
do i=1,ndim
Bmat(i,i)= zone
enddo
call zgesv(ndim,ndim,Amat,ndim,ipiv,Bmat,ndim,info)
if(info.ne.0)print *,'something wrong with zgesv'
Amat=Bmat
return
end subroutine inv
!> get the inverse of a real matrix
subroutine inv_r(ndim,Amat)
implicit none
integer,parameter :: dp=8
integer :: i
integer :: info
integer,intent(in):: ndim
! IPIV : INTEGER. Array, DIMENSION at least max(1, n). The pivot indices that define
! the permutation matrix P; row i of the matrix was interchanged with
! row ipiv(i). Corresponds to the single precision factorization (if info=
! 0 and iter ≥ 0) or the double precision factorization (if info= 0 and
! iter < 0).
integer,allocatable :: ipiv(:)
! Amat :
! Overwritten by the factors L and U from the factorization of A = P*L*U;
! the unit diagonal elements of L are not stored.
! If iterative refinement has been successfully used (info= 0 and iter≥
! 0), then A is unchanged.
! If double precision factorization has been used (info= 0 and iter <
! 0), then the array A contains the factors L and U from the factorization
! A = P*L*U; the unit diagonal elements of L are not stored.
real(dp),intent(inout):: Amat(ndim,ndim)
! Bmat :
! Overwritten by the solution matrix X for dgesv, sgesv,zgesv,zgesv.
real(dp),allocatable :: Bmat(:,:)
allocate(ipiv(ndim))
allocate(Bmat(ndim,ndim))
ipiv=0
! unit matrix
Bmat= 0d0
do i=1,ndim
Bmat(i,i)= 1d0
enddo
call dgesv(ndim,ndim,Amat,ndim,ipiv,Bmat,ndim,info)
if(info.ne.0)stop 'something wrong with dgesv in inverse.f90'
Amat=Bmat
return
end subroutine inv_r