Contents of file 'zgaussj.f':
1 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
2 C File: zgaussj.f (Fortran 77 source for the zgaussj routine)
3 C Author: Fredrik Jonsson <fj@phys.soton.ac.uk>
4 C Date: December 26, 2005
5 C Last change: December 29, 2006
6 C Description: The zgaussj(a,n,np,b,m,mp) routine solves systems of
7 C linear equations by Gauss-Jordan elimination, in similar
8 C to as in Eq. (2.1.1) of W.H. Press et al, "Numerical
9 C Recipes in Fortran 77, The Art of Scientific Computing",
10 C 2nd Edn. (Cambridge University Press, Cambridge, 1992),
11 C but here modified so as to solve complex-valued systems.
12 C Here a(1:n,1:n) is an input matrix stored in an array
13 C of physical dimensions np by np. b(1:n,1:m) is an input
14 C matrix containing the m right-hand side vectors, stored
15 C in an array of physical dimensions np by mp. On output,
16 C a(1:n,1:n) is replaced by its matrix inverse, and
17 C b(1:n,1:m) is replaced by the corresponding set of solu-
18 C tion vectors. The parameter nmax contains the largest
19 C anticipated value of n.
20 C
21 C Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
22 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
23 subroutine zgaussj(a,n,np,b,m,mp)
24 integer m,mp,n,np,nmax
25 complex*8 a(np,np),b(np,mp)
26 parameter (nmax=50)
27 integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
28 complex*8 big,dum,pivinv
29 do 11 j=1,n
30 ipiv(j)=0
31 11 continue
32 do 22 i=1,n
33 big=0.
34 do 13 j=1,n
35 if(ipiv(j).ne.1)then
36 do 12 k=1,n
37 if (ipiv(k).eq.0) then
38 if (abs(a(j,k)).ge.abs(big))then
39 big=abs(a(j,k))
40 irow=j
41 icol=k
42 endif
43 else if (ipiv(k).gt.1) then
44 pause 'Singular matrix ecountered in zgaussj()!'
45 endif
46 12 continue
47 endif
48 13 continue
49 ipiv(icol)=ipiv(icol)+1
50 if (irow.ne.icol) then
51 do 14 l=1,n
52 dum=a(irow,l)
53 a(irow,l)=a(icol,l)
54 a(icol,l)=dum
55 14 continue
56 do 15 l=1,m
57 dum=b(irow,l)
58 b(irow,l)=b(icol,l)
59 b(icol,l)=dum
60 15 continue
61 endif
62 indxr(i)=irow
63 indxc(i)=icol
64 if (abs(a(icol,icol)).eq.(0.0))
65 pause 'Singular matrix ecountered in zgaussj()! [Case 2]'
66 pivinv=1.0/a(icol,icol)
67 a(icol,icol)=1.
68 do 16 l=1,n
69 a(icol,l)=a(icol,l)*pivinv
70 16 continue
71 do 17 l=1,m
72 b(icol,l)=b(icol,l)*pivinv
73 17 continue
74 do 21 ll=1,n
75 if(ll.ne.icol)then
76 dum=a(ll,icol)
77 a(ll,icol)=0.
78 do 18 l=1,n
79 a(ll,l)=a(ll,l)-a(icol,l)*dum
80 18 continue
81 do 19 l=1,m
82 b(ll,l)=b(ll,l)-b(icol,l)*dum
83 19 continue
84 endif
85 21 continue
86 22 continue
87 do 24 l=n,1,-1
88 if(indxr(l).ne.indxc(l))then
89 do 23 k=1,n
90 dum=a(k,indxr(l))
91 a(k,indxr(l))=a(k,indxc(l))
92 a(k,indxc(l))=dum
93 23 continue
94 endif
95 24 continue
96 return
97 end
98
Generated by ::viewsrc::