我编程用共轭梯度法解方程，结果运行时出现以下提示，请问如何改正？


Loaded 'ntdll.dll', no matching symbolic information found.
Loaded 'C:\WINDOWS\system32\kernel32.dll', no matching symbolic information found.
First-chance exception in gongetidufa.exe: 0xC0000005: Access Violation.
The thread 0x1C00 has exited with code -1073741667 (0xC000009D).
The program 'E:\R_K Method\jisuanliutizuoye\Debug\gongetidufa.exe' has exited with code -1073741667 (0xC000009D).


源程序：
subroutine sparse(b,n,asub,atsub,x,rsq)
parameter(nmax=500,eps=1.e-6)
dimension b(n),x(n),g(nmax),h(nmax),xi(nmax),xj(nmax)
logical done
eps2=n*eps**2
irst=0
do
done=-1
irst=irst+1
call asub(x,xi)
rp=0.0
bsq=0.0
do j=1,n
bsq=bsq+b(j)**2
xi(j)=xi(j)-b(j)
rp=rp+xi(j)**2
enddo
call atsub(xi,g)
do j=1,n
g(j)=-g(j)
h(j)=g(j)
enddo
do iter=1,10*n
call asub(h,xi)
anum=0.
aden=0.
do j=1,n
anum=anum+g(j)*h(j)
aden=aden+xi(j)**2
enddo
if(aden==0.0)pause'very singular matrix'
anum=anum/aden
do j=1,n
xi(j)=x(j)
x(j)=x(j)+anum*h(j)
enddo
call asub(x,xj)
rsq=0.
do j=1,n
xj(j)=xj(j)-b(j)
rsq=rsq+xj(j)**2
enddo
if(rsq==rp.or.rsq<=bsq*eps2)return
if(rsq>rp)then
do j=1,n
x(j)=xi(j)
enddo
if(irst>=3)return
done=0
endif
if(.not.done)exit
rp=rsq
call atsub(xj,xi)
gg=0.0
dgg=0.0
do j=1,n
gg=gg+g(j)**2
dgg=dgg+(xi(j)+g(j))*xi(j)
enddo
if(gg==0.)return
gam=dgg/gg
do j=1,n
g(j)=-xi(j)
h(j)=g(j)+gam*h(j)
enddo
enddo
if(.not.done)exit
enddo
pause 'too many iterations'
end subroutine sparse

program d1r9
!driver program for routine sparse
parameter(n=10000,nn=100,deltx=0.01,delty=0.01)
common m
dimension b(n),x(n),bcmp(n),a(n,n),bb(nn,nn)
external asub,atsub
do i=1,nn
   do j=1,nn
    a((i-1)*nn+j,(j-1)*nn+i)=-(2+2.0*deltx**2/delty**2)
    a((i-1)*nn+j,(j-2)*nn+i)=deltx**2/delty**2
    a((i-1)*nn+j,j*nn+i)=deltx**2/delty**2
    a((i-1)*nn+j,(j-1)*nn+i-1)=1.0
    a((i-1)*nn+j,(j-1)*nn+i+1)=1.0
    enddo
enddo
do i=1,nn
    do j=1,nn
    bb(i,j)=deltx**2*sin(i*deltx)*cos(j*delty)
    b((i-1)*nn+j)=bb(i,j)
    enddo
enddo
m=n
do i=1,n
x(i)=0.0
b(i)=b(i)
enddo
call sparse(b,n,asub,atsub,x,rsq)
write(*,'(/1x,a,e15.6)')'sum-squared residual:',rsq
write(*,'(/1x,a)')'solution vector:'
write(*,'(/1x,5f12.6)')(x(i),i=1,n)
call asub(x,bcmp)
write(*,'(/1x,a)')'press return to end do...'
read(*,*)
write(*,'(1x,a/t8,a,t22,a)')'test of solution vector:','a*x','b'
do i=1,n**2
write(*,'(1x,2f12.6)')bcmp(i),b(i)
enddo
end program
subroutine asub(xin,xout)
common n
dimension xin(n),xout(n)
xout(1)=xin(i+1)-(2+2.0*deltx**2/delty**2)*xin(i)+deltx**2/delty**2*xin(101)
xout(n)=xin(9999)-(2+2.0*deltx**2/delty**2)*xin(10000)+deltx**2/delty**2*xin(9900)
do i=2,nn-1
   do j=2,nn-1
xout(i)=xin((j-1)*nn+i-1)+xin((j-1)*nn+i+1)-(2+2.0*deltx**2/delty**2)*xin((j-1)*nn+i)&
+deltx**2/delty**2*xin((j-2)*nn+i)+deltx**2/delty**2*xin(j*nn+i)
enddo
enddo
end subroutine asub
subroutine atsub(xin,xout)
common n
dimension xin(n),xout(n)
xout(1)=xin(2)-(2+2.0*deltx**2/delty**2)*xin(1)+deltx**2/delty**2*xin(101)
xout(n)=xin(9999)-(2+2.0*deltx**2/delty**2)*xin(10000)+deltx**2/delty**2*xin(9900)
do i=2,nn-1
   do j=2,nn-1
xout(i)=xin((i-1)*nn+j-1)+xin((i-1)*nn+j+1)-(2+2.0*deltx**2/delty**2)*xin((i-1)*nn+j)&
+deltx**2/delty**2*xin((i-2)*nn+j)+deltx**2/delty**2*xin(i*nn+j)
enddo
enddo
end subroutine atsub