Search:

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   ```

```