1
2!***** Double precision linear algebra *******************************
3
4      subroutine loupd(x,n,nm,iperm,nperm)
5! Lower-upper decomposition of matrix x, logical size n x n,
6! stored as size nm x nm.  Returns matrix x(i,j),i<=j as upper
7! triangular matrix and x(i,j),i>j as lower diagonal matrix with 1 on
8! the diagonal.  The product of these two gives the original matrix x.
9! Nperm is the order of the permutations of rows, either +1 or -1,
10! unless the matrix is singular, in which case nperm is returned as 0
11! and evaluation of the subroutine is stopped.
12      implicit none
13      integer n,nm,nperm
14      integer iperm(nm)
15      integer i,ii,j,jj,k,kk,im
16      double precision x(nm,nm),temp
17      double precision cmax,scl(nm),check
18! scl(i) is scaling factor for row i, to normalize largest element to 1
19      nperm=1
20      do i=1,n
21        cmax=0.d0
22        do j=1,n
23          if(abs(x(i,j)).gt.cmax) cmax=abs(x(i,j))
24        enddo
25        if (cmax.eq.0.d0) then
26          nperm=0
27          return
28        endif
29        scl(i)=1.d0/cmax
30      enddo
31      do j=1,n
32        do i=1,j-1
33          temp=x(i,j)
34          do k=1,i-1
35            temp=temp-x(i,k)*x(k,j)
36          enddo
37          x(i,j)=temp
38        enddo
39        cmax=0.d0
40        do i=j,n
41          temp=x(i,j)
42          do k=1,j-1
43            temp=temp-x(i,k)*x(k,j)
44          enddo
45          x(i,j)=temp
46          check=scl(i)*abs(temp)
47          if (check.ge.cmax) then
48            im=i
49            cmax=check
50          endif
51        enddo
52        if (j.ne.im) then
53          do k=1,n
54            temp=x(im,k)
55            x(im,k)=x(j,k)
56            x(j,k)=temp
57          enddo
58          nperm=-nperm
59          check=scl(j)
60          scl(j)=scl(im)
61          scl(im)=check
62        endif
63        iperm(j)=im
64        if(j.ne.n) then
65          if (x(j,j).ne.0) then
66            temp=(1.d0,0.d0)/x(j,j)
67            do i=j+1,n
68              x(i,j)=x(i,j)*temp
69            enddo
70          else
71            nperm=0
72            return
73          endif
74        endif
75      enddo
76      return
77      end
78
79      subroutine baksublu(x,v,n,nm,iperm)
80! Solves a linear equation x'*w=v for w.  X is the lower-upper
81! decomposition of x' from subroutine loupd.  Vector v is replaced by
82! the solutions w as output.
83      implicit none
84      integer n,nm,iperm(nm),i,j,ii,jj
85      double precision x(nm,nm),v(nm),temp
86      ii=0
87      do i=1,n
88        temp=v(iperm(i))
89        v(iperm(i))=v(i)
90        if (ii.ne.0) then
91          do j=ii,i-1
92            temp=temp-x(i,j)*v(j)
93          enddo
94        elseif (temp.ne.0.d0) then
95          ii=i
96        endif
97        v(i)=temp
98      enddo
99      do i=n,1,-1
100        temp=v(i)
101        do j=i+1,n
102          temp=temp-x(i,j)*v(j)
103        enddo
104        v(i)=temp/x(i,i)
105      enddo
106      return
107      end
108
109
110      subroutine inverse(x,y,n,nm)
111! Returns matrix y as the inverse of matrix x, both
112! of dimension n x n.
113      implicit none
114      integer n,nm
115      double precision x(nm,nm),y(nm,nm),v(nm)
116      integer i,j,k,ii,jj,kk
117      integer iperm(nm),nperm
118      do i=1,n
119        do j=1,n
120          if (i.eq.j) then
121            y(i,j)=1.d0
122          else
123            y(i,j)=0.d0
124          endif
125        enddo
126      enddo
127      call loupd(x,n,nm,iperm,nperm)
128      do j=1,n
129        call baksublu(x,y(1,j),n,nm,iperm)
130      enddo
131      return
132      end
133
134!***** Complex linear algebra *********************************************
135
136      subroutine loupdc(x,n,nm,iperm,nperm)
137! Lower-upper decomposition of matrix x, logical size n x n,
138! stored as size nm x nm.  Returns matrix x(i,j),i<=j as upper
139! triangular matrix and x(i,j),i>j as lower diagonal matrix with 1 on
140! the diagonal.  The product of these two gives the original matrix x.
141! Nperm is the order of the permutations of rows, either +1 or -1,
142! unless the matrix is singular, in which case nperm is returned as 0
143! and evaluation of the subroutine is stopped.
144      implicit none
145      integer n,nm,nperm
146      integer iperm(nm)
147      integer i,ii,j,jj,k,kk,im
148      double complex x(nm,nm),temp
149      double precision cmax,scl(nm),check
150! scl(i) is scaling factor for row i, to normalize largest element to 1
151      nperm=1
152      do i=1,n
153        cmax=0.d0
154        do j=1,n
155          if(abs(x(i,j)).gt.cmax) cmax=abs(x(i,j))
156        enddo
157        if (cmax.eq.0.d0) then
158          nperm=0
159          return
160        endif
161        scl(i)=1.d0/cmax
162      enddo
163      do j=1,n
164        do i=1,j-1
165          temp=x(i,j)
166          do k=1,i-1
167            temp=temp-x(i,k)*x(k,j)
168          enddo
169          x(i,j)=temp
170        enddo
171        cmax=0.d0
172        do i=j,n
173          temp=x(i,j)
174          do k=1,j-1
175            temp=temp-x(i,k)*x(k,j)
176          enddo
177          x(i,j)=temp
178          check=scl(i)*abs(temp)
179          if (check.ge.cmax) then
180            im=i
181            cmax=check
182          endif
183        enddo
184        if (j.ne.im) then
185          do k=1,n
186            temp=x(im,k)
187            x(im,k)=x(j,k)
188            x(j,k)=temp
189          enddo
190          nperm=-nperm
191          check=scl(j)
192          scl(j)=scl(im)
193          scl(im)=check
194        endif
195        iperm(j)=im
196        if(j.ne.n) then
197          if (x(j,j).ne.0) then
198            temp=(1.d0,0.d0)/x(j,j)
199            do i=j+1,n
200              x(i,j)=x(i,j)*temp
201            enddo
202          else
203            nperm=0
204            return
205          endif
206        endif
207      enddo
208      return
209      end
210
211      subroutine baksubluc(x,v,n,nm,iperm)
212! Solves a linear equation x'*w=v for w.  X is the lower-upper
213! decomposition of x' from subroutine loupdc.  Vector v is replaced by
214! the solutions w as output.
215      implicit none
216      integer n,nm,iperm(nm),i,j,ii,jj
217      double complex x(nm,nm),v(nm),temp
218      ii=0
219      do i=1,n
220        temp=v(iperm(i))
221        v(iperm(i))=v(i)
222        if (ii.ne.0) then
223          do j=ii,i-1
224            temp=temp-x(i,j)*v(j)
225          enddo
226        elseif (temp.ne.0.d0) then
227          ii=i
228        endif
229        v(i)=temp
230      enddo
231      do i=n,1,-1
232        temp=v(i)
233        do j=i+1,n
234          temp=temp-x(i,j)*v(j)
235        enddo
236        v(i)=temp/x(i,i)
237      enddo
238      return
239      end
240
241
242      subroutine inversec(x,y,n,nm)
243! Returns matrix y as the inverse of matrix x, both
244! of dimension n x n.
245      implicit none
246      integer n,nm
247      double complex x(nm,nm),y(nm,nm),v(nm)
248      integer i,j,k,ii,jj,kk
249      integer iperm(nm),nperm
250!
251      double complex ll(nm,nm),uu(nm,nm),AA(nm,nm)
252      do i=1,n
253        do j=1,n
254          if (i.eq.j) then
255            y(i,j)=1.d0
256          else
257            y(i,j)=0.d0
258          endif
259        enddo
260      enddo
261      call loupdc(x,n,nm,iperm,nperm)
262      do j=1,n
263        call baksubluc(x,y(1,j),n,nm,iperm)
264      enddo
265      return
266      end
267
268
269