对于李建平老师的EOF中JACOBI子程序的疑问

centaurus

想问一下知道的人，我在用李老师的EOF程序做SST时，发现程序到了JACOBI那进入死循环了 ，底下是部分程序。我发现我的l=0，这导致程序一直回到60那重复。因为对JACOBI的数学表达不是很了解。想知道这个l具体是什么？还有怎么解决这种问题呢？subroutine jacobi(m,a,s,d,eps)                implicit none
                integer*4                ::        m
                real*4                        ::        a(m,m),s(m,m),d(m),eps
                integer*4                ::        i,j,i1,l,iq,iq1,ip
                real*4                        ::        g,s1,s2,s3,v1,v2,v3,u,st,ct
                        s=0.0
                        do i=1,m
                                s(i,i)=1.0
                        enddo
                        g=0.
                        do i=2,m
                                i1=i-1
                                do j=1,i1
                                g=g+2.0*a(i,j)*a(i,j)
                                enddo
                        enddo

                        s1=sqrt(g)
                        s2=eps/float(m)*s1
                        s3=s1
                        l=0
50                        s3=s3/float(m)
60                        do iq=2,m
                                iq1=iq-1                                       
                                do ip=1,iq1
                                        if(abs(a(ip,iq)).lt.s3) exit
                                        l=1
                                        v1=a(ip,ip)
                                        v2=a(ip,iq)
                                        v3=a(iq,iq)
                                        u=0.5*(v1-v3)
                                        if(u.eq.0.0) g=1.
                                        if(abs(u).ge.1e-10) g=-sign(1.,u)*v2/sqrt(v2*v2+u*u)
                                        st=g/sqrt(2.*(1.+sqrt(1.-g*g)))
                                        ct=sqrt(1.-st*st)
                                        do i=1,m
                                                g=a(i,ip)*ct-a(i,iq)*st
                                                a(i,iq)=a(i,ip)*st+a(i,iq)*ct
                                                a(i,ip)=g
                                                g=s(i,ip)*ct-s(i,iq)*st
                                                s(i,iq)=s(i,ip)*st+s(i,iq)*ct
                                                s(i,ip)=g
                                        enddo
                                        do i=1,m
                                                a(ip,i)=a(i,ip)
                                                a(iq,i)=a(i,iq)
                                        enddo
                                        g=2.*v2*st*ct
                                        a(ip,ip)=v1*ct*ct+v3*st*st-g
                                        a(iq,iq)=v1*st*st+v3*ct*ct+g
                                        a(ip,iq)=(v1-v3)*st*ct+v2*(ct*ct-st*st)
                                        a(iq,ip)=a(ip,iq)
                                enddo
                          enddo
                          write(*,*)l
                          stop
                        if((l-1)==0)then
                                l=0
                                go to 60
                        else
                                go to 150
                        endif

150               

centaurus

错了，我写错了，是l=1，然后IF语句那l-1=0，就回头重复去了

棋子灯花

。。。二爷碉堡了！！

centaurus

棋子灯花 发表于 2014-3-28 14:48
。。。二爷碉堡了！！

艰难。。。。

言深深

JCB是求解行列式的···
fortran是不会死循环的；
l=0，你看看是不是哪里没有执行···

wnyh

是吗？不会吧

L-S-Y

jab是行列式，计算距阵的

scwjx520

做EOF是不是很难啊

mmwj1990

我也出现了同样的问题，请问您的最后怎么解决了？非常感谢！

centaurus

mmwj1990 发表于 2014-9-13 11:02
我也出现了同样的问题，请问您的最后怎么解决了？非常感谢！

不好意思，过去了将近半年也给忘了。