请教高斯解方程中gaussj(a,n,np,b,m,mp)中参数a,n,np,b,m,mp分别代表什么
请教高斯解方程中gaussj(a,n,np,b,m,mp)中参数a,n,np,b,m,mp分别代表什么?
SUBROUTINE gaussj(a,n,np,b,m,mp)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=4000)
C INTEGER m,mp,n,np,NMAX
c REAL a(np,np),b(np,mp)
DIMENSION a(np,np),b(np,mp)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
c REAL big,dum,pivinv
C DOUBLE PRECISION big,dum,pivinv
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
c WRITE(*,*) "I=",I
big=0.
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif
C WRITE(*,*) I,J
else if (ipiv(k).gt.1) then
WRITE(*,*) I,J
pause 'is singular matrix in gaussj'
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
pivinv=1./a(icol,icol)
a(icol,icol)=1.
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
END
SUBROUTINE gaussj(a,n,np,b,m,mp)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=4000)
C INTEGER m,mp,n,np,NMAX
c REAL a(np,np),b(np,mp)
DIMENSION a(np,np),b(np,mp)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
c REAL big,dum,pivinv
C DOUBLE PRECISION big,dum,pivinv
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
c WRITE(*,*) "I=",I
big=0.
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif
C WRITE(*,*) I,J
else if (ipiv(k).gt.1) then
WRITE(*,*) I,J
pause 'is singular matrix in gaussj'
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
pivinv=1./a(icol,icol)
a(icol,icol)=1.
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
END