File: /nfs/home/0/users/jenkins/mfix.git/model/dqmom/gaussj.f
1 SUBROUTINE gaussj(a,n,np,b,m,mp)
2 INTEGER m,mp,n,np,NMAX
3 DOUBLE PRECISION a(np,np),b(np,mp)
4 PARAMETER (NMAX=50)
5 INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
6 DOUBLE PRECISION big,dum,pivinv
7
8 ipiv(:)=0
9
10 do i=1,n
11 big=0.
12 do j=1,n
13 if(ipiv(j)/=1)then
14 do k=1,n
15 if (ipiv(k)==0) then
16 if (abs(a(j,k))>=big)then
17 big=abs(a(j,k))
18 irow=j
19 icol=k
20 endif
21 else if (ipiv(k)>1) then
22 write(*,*) 'WARNING: singular matrix in gaussj'
23 endif
24 enddo
25 endif
26 enddo
27 ipiv(icol)=ipiv(icol)+1
28 if (irow/=icol) then
29 do l=1,n
30 dum=a(irow,l)
31 a(irow,l)=a(icol,l)
32 a(icol,l)=dum
33 enddo
34 do l=1,m
35 dum=b(irow,l)
36 b(irow,l)=b(icol,l)
37 b(icol,l)=dum
38 enddo
39 endif
40 indxr(i)=irow
41 indxc(i)=icol
42 if (a(icol,icol)==0.) write(*,*) 'WARNING: singular matrix in gaussj'
43 pivinv=1./a(icol,icol)
44 a(icol,icol)=1.
45 do l=1,n
46 a(icol,l)=a(icol,l)*pivinv
47 enddo
48 do l=1,m
49 b(icol,l)=b(icol,l)*pivinv
50 enddo
51 do ll=1,n
52 if(ll.ne.icol)then
53 dum=a(ll,icol)
54 a(ll,icol)=0.
55 do l=1,n
56 a(ll,l)=a(ll,l)-a(icol,l)*dum
57 enddo
58 do l=1,m
59 b(ll,l)=b(ll,l)-b(icol,l)*dum
60 enddo
61 endif
62 enddo
63 enddo
64 do l=n,1,-1
65 if(indxr(l)/=indxc(l))then
66 do k=1,n
67 dum=a(k,indxr(l))
68 a(k,indxr(l))=a(k,indxc(l))
69 a(k,indxc(l))=dum
70 enddo
71 endif
72 enddo
73 return
74 END
75