您所在位置：主题：为何0乘以一个常数输的是NaN
siyuan
[专家分：0] 发布于 2010-05-25 19:03:00
下面的程序，提示没有错误，但运行没结果。我一步一步查看，发现第一次call asub(x,xi)(在第十行)就出了问题，此时x全部为零，xi等于部分x乘以一个常数，结果应该也全部为零，但是结果大部分是NaN，只有极少数为0，感觉好奇怪，看了半天也没找到原因。谁可以指出问题么
subroutine sparse(b,n,asub,atsub,x,rsq)
parameter(nmax=10000,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)
write(*,*) 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+1
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
real deltx,delty
parameter(n=99**2,deltx=0.01,delty=0.01)
common m
dimension b(n),x(n),bcmp(n)
external asub,atsub
m=n
do i=1,98
b(i)=sin(i*deltx)*cos(delty)-0.5*sin(i*deltx)/delty**2
enddo
b(99)=sin(99*deltx)*cos(delty)-0.5*sin(99*deltx)/delty**2+sin(1.0)*cos(delty)/deltx**2
do i=2,98
b((i-1)*99+99)=sin(99*deltx)*cos(i*delty)+(i*delty-0.2*sin(1.0)*cos(i*delty))/deltx**2
enddo
do i=2,98
do j=1,98
b((i-1)*99+j)=sin(j*deltx)*cos(i*delty)
enddo
enddo
do i=99*98+1,99*99-1
b(i)=sin((i-98*99)*deltx)*cos(99*delty)+((i-98*99)*deltx-0.5*sin((i-98*99)*deltx)*cos(1.0))/delty**2
enddo
b(99*99)=sin(99*deltx)*cos(99*delty)+(99*deltx-0.5*sin(99*deltx)*cos(1.0))/delty**2+(99*delty-0.5*sin(1.0)*cos(99*delty))/deltx**2
do i=1,99*99
x(i)=0.0
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
write(*,'(1x,2f12.6)') bcmp(i),b(i)
enddo
end program
subroutine asub(xin,xout)
common n
dimension xout(n),xin(n)
xout(1)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(1)-xin(2)/deltx**2-xin(100)/delty**2
xout(99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99)-xin(98)/deltx**2-xin(198)/delty**2
do i=2,98
xout(i)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(i)-xin(i-1)/deltx**2-xin(i+1)/deltx**2-xin(i+99)/delty**2
enddo
do i=2,98
xout((i-1)*99+1)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+1)-xin((i-1)*99+2)/deltx**2-xin((i-1)*99+100)/delty**2-xin((i-1)*99-98)/delty**2
xout((i-1)*99+99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+99)-xin((i-1)*99+98)/deltx**2-xin((i-1)*99+198)/delty**2-xin((i-1)*99)/delty**2
enddo
do i=2,98
do j=2,98
xout((i-1)*99+j)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+j)-xin((i-1)*99+j+1)/deltx**2-xin((i-1)*99+j-1)/deltx**2-xin((i-1)*99+j+99)/delty**2-xin((i-1)*99+j-99)/delty**2
enddo
enddo
xout(99*99-98)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99*99-98)-xin(99*99-97)/deltx**2-xin(99*99-197)/delty**2
do i=99*98+2,99*98+98
xout(i)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(i)-xin(i-1)/deltx**2-xin(i+1)/deltx**2-xin(i-99)/delty**2
enddo
xout(99*99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99*99)-xin(99*99-1)/deltx**2-xin(99*99-99)/delty**2
end subroutine asub
subroutine atsub(xin,xout)
common n
dimension xin(n),xout(n)
xout(1)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(1)-xin(2)/deltx**2-xin(100)/delty**2
xout(99*99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99*99)-xin(99*99-1)/deltx**2-xin(99*99-99)/delty**2
xout(99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99)-xin(98)/deltx**2-xin(198)/delty**2
xout(99*99-98)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(99*99-98)-xin(99*99-97)/deltx**2-xin(99*99-197)/delty**2
do i=2,98
xout(i)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(i)-xin(i-1)/deltx**2-xin(i+1)/deltx**2-xin(i+99)/delty**2
enddo
do i=99*98+2,99*98+98
xout(i)=2.0*(1.0/deltx**2+1.0/delty**2)*xin(i)-xin(i-1)/deltx**2-xin(i+1)/deltx**2-xin(i-99)/delty**2
enddo
do i=2,98
xout((i-1)*99+1)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+1)-xin((i-1)*99+2)/deltx**2-xin((i-1)*99+100)/delty**2-xin((i-1)*99-98)/delty**2
xout((i-1)*99+99)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+99)-xin((i-1)*99+98)/deltx**2-xin((i-1)*99+198)/delty**2-xin((i-1)*99)/delty**2
enddo
do i=2,98
do j=2,98
xout((i-1)*99+j)=2.0*(1.0/deltx**2+1.0/delty**2)*xin((i-1)*99+j)-xin((i-1)*99+j+1)/deltx**2-xin((i-1)*99+j-1)/deltx**2-xin((i-1)*99+j+99)/delty**2-xin((i-1)*99+j-99)/delty**2
enddo
enddo
end subroutine atsub
[em4][em4][em4][em18][em18]
