[Fortran] 纯文本查看 复制代码
subroutine gband(a,d,x,n,m,eps,ierr,ifrst)
include'param.f'
double precision a,d,x,eps,s
dimension a(14*idim-12),d(idim),x(idim)
integer n,m,ierr,ifrst,j,k,i,ie,ieux,ie1,mbig,j1,j2,jo,j1k,
& iaux,ieaux,ik,jk,np1,iinv
ierr =0
j=1
do 10 i=1,n
ie=m
if(i+m-n)21,21,22
22 ie=n-i
21 ieaux=m
if(i-m)23,23,24
23 ieaux=i
24 ie1=ie+ieaux
mbig=ie
j1=j+ie1
j2=j1
if(ifrst.gt.0)goto 27
if(dabs(a(j))-eps)25,25,27
25 ierr=ierr+1
27 if(mbig)10,10,26
26 do 20 jo=1,mbig
s=a(j1)/a(j)
if(ifrst.gt.0)goto 35
do 30 k=1,mbig
j1k=j1+k
jk=j+k
30 a(j1k)=a(j1k)-a(jk)*s
35 continue
iaux=jo+i
d(iaux)=d(iaux)-d(i)*s
ie=m
if(iaux+m-n)31,31,32
32 ie=n-iaux
31 ieaux=m
if(iaux-m)33,33,34
33 ieaux=iaux
34 ie1=ie+ieaux
20 j1=j1+ie1
10 j=j2+1
j=j-m-1
np1=n+1
c
do 40 iinv=1,n
i=np1-iinv
ie=m
if(i+m-n)41,41,42
42 ie=n-i
41 mbig=ie
x(i)=d(i)
if(mbig)44,44,43
c
43 do 50 k=1,mbig
ik=i+k
jk=j+k
50 x(i)=x(i)-x(ik)*a(ik)
c
44 x(i)=x(i)/a(j)
ie=m
if(i+m-np1)51,51,52
52 ie=np1-i
51 ieaux=m
if(i-1-m)53,53,54
53 ieaux=i-1
54 ie1=ie+ieaux
j=j-ie1-1
40 continue
return
end