-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmatinv.f
130 lines (127 loc) · 3.66 KB
/
matinv.f
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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
c Truong Dang added 02-11-2016
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc ccc
ccc purpose: Inverse computation of a matrix ccc
ccc ccc
ccc Input and output varibles: ccc
ccc input n size of matrix ccc
ccc input a input matrix ccc
ccc output y output inverse matrix ccc
ccc ccc
ccc Use: ccc
ccc ludcmp LU docomposition ccc
ccc lubksb find inverse of matrix by columns ccc
ccc ccc
ccc Warning: After matinv(a,n,y) is called, the input matrix ccc
ccc 'a' would be changed (used as working space) ccc
ccc unless it is a diagnal matix. So if 'a' is ccc
ccc not a diagnal matrix and will be used after ccc
ccc calling matinv.f, 'a' should be saved before ccc
ccc calling matinv.f. ccc
ccc ccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine matinv(a,n,y)
implicit none
c
integer(4)::n,i,j
real(8):: a(n,n),y(n,n),d
integer(4)::indx(n)
y=0.
do i=1,n
y(i,i)=1.
end do
call ludcmp(a,n,n,indx,d)
do j=1,n
call lubksb(a,n,n,indx,y(1,j))
end do
c
return
end
c
SUBROUTINE ludcmp(a,n,np,indx,d)
implicit none
INTEGER(4):: n,np,indx(n),NMAX
REAL(8):: d,a(np,np),TINY
PARAMETER (NMAX=100,TINY=1.0e-20)
INTEGER(4):: i,imax,j,k
REAL(8):: aamax,dum,sum,vv(n)
d=1.0d0
do 12 i=1,n
aamax=0.
do 11 j=1,n
if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11 continue
if (aamax.eq.0.) print*, 'singular matrix in ludcmp'
vv(i)=1./aamax
12 continue
do 19 j=1,n
do 14 i=1,j-1
sum=a(i,j)
do 13 k=1,i-1
sum=sum-a(i,k)*a(k,j)
13 continue
a(i,j)=sum
14 continue
aamax=0.
do 16 i=j,n
sum=a(i,j)
do 15 k=1,j-1
sum=sum-a(i,k)*a(k,j)
15 continue
a(i,j)=sum
dum=vv(i)*abs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
16 continue
if (j.ne.imax)then
do 17 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
17 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.)a(j,j)=TINY
if(j.ne.n)then
dum=1./a(j,j)
do 18 i=j+1,n
a(i,j)=a(i,j)*dum
18 continue
endif
19 continue
return
END
c
SUBROUTINE lubksb(a,n,np,indx,b)
implicit none
INTEGER(4):: n,np,indx(n)
REAL(8):: a(np,np),b(n)
INTEGER(4):: i,ii,j,ll
REAL(8):: sum
ii=0
do 12 i=1,n
ll=indx(i)
sum=b(ll)
b(ll)=b(i)
if (ii.ne.0)then
do 11 j=ii,i-1
sum=sum-a(i,j)*b(j)
11 continue
else if (sum.ne.0.) then
ii=i
endif
b(i)=sum
12 continue
do 14 i=n,1,-1
sum=b(i)
do 13 j=i+1,n
sum=sum-a(i,j)*b(j)
13 continue
b(i)=sum/a(i,i)
14 continue
return
END