program main implicit none integer::N,info double precision::h,tol double precision::ta,tb,t double precision,allocatable::y(:) double precision,external::grk !$call mkl_set_num_threads(1) !$call omp_set_num_threads(2) N=2 ! Number of equations allocate(y(1:N)) ta=0d0 ! Start time tb=50d0 ! End time !initial conditions y(1 )=1d0 ! y (ta) y(2 )=0d0 ! y'(ta) t=ta tol=1d-10 info=-1 h=tb-ta write(10,'(4e30.17e3)')t,y(1),y(2),h do while(info.le.0) call irk64(grk,t,h,N,y,tb,info,tol) write(10,'(4e30.17e3)')t,y(1),y(2),h enddo stop end program main function grk(N,x,y,s) implicit none integer,intent(in)::N,s double precision,intent(in)::x double precision,intent(in)::y(1:N) double precision::grk grk=0d0 if(s.eq.1)then grk = y(2) elseif(s.eq.2)then grk = - y(1) else write(6,*)"***Error grk"; stop endif return end function grk subroutine irk64(grk,x,h,N,y,xbound,info,tol) implicit none integer,intent(in)::N double precision,intent(in)::xbound,tol double precision,intent(inout)::x,h,y(1:N) integer,intent(inout)::info double precision,external::grk double precision,allocatable::y4(:),y6(:) double precision,allocatable::dy4(:),dy6(:) integer::i,FLAG,key double precision::R0,R,delta,tx,tx0,Sy,err double precision,parameter::hmin=1d-14,hmax=0.1d0 allocate(y6(1:N),y4(1:N)) y6(1:N)=y(1:N) y4(1:N)=y(1:N) allocate(dy4(1:N),dy6(1:N)) !----------- key=0 if(abs(h).ge.hmax)then h=sign(1d0,h)*hmax endif if(h.ge.abs(xbound-x))h=xbound-x FLAG=1 if(abs(x-xbound).le.hmin)then info=1 FLAG=0 endif Sy=0d0 do i=1,N Sy=Sy+abs(y(i)) enddo if(Sy.ge.1d0)then err=tol*Sy else err=tol endif do while(FLAG.eq.1) R=0d0 R0=0d0 tx=x y6(1:N)=y(1:N) tx0=x y4(1:N)=y(1:N) !$omp parallel !$omp sections !$omp section call dirk6_sub(grk,tx,h,N,y6,dy6) !$omp section call dirkL3C4_sub(grk,tx0,h,N,y4,dy4) !$omp end sections !$omp end parallel do i=1,N R = R + dy6(i) R0 = R0 + dy4(i) enddo R = abs((R-R0)) !step 5 if(R.le.err.or.key.eq.1)then x=x+h y(1:N)=y6(1:N) FLAG=0 endif if(R.ge.1d-20)then delta=(err/(2d0*R))**0.20d0 else delta=4d0 endif !step 7 if(delta.le.0.1d0)then !function changes dramatically. h=0.1d0*h elseif(delta.ge.4d0)then !function changes loosely. h=4d0*h else !function changes moderately. h=delta*h endif !step 8 if(abs(h).ge.hmax)then h=sign(1d0,h)*hmax elseif(abs(h).lt.hmin)then h=sign(1d0,h)*hmin key=1 endif !step 9 if(abs(xbound-x).le.abs(h))then h=xbound-x if(abs(h).le.hmin)then info=1 FLAG=0 endif end if if(h.le.0d0.and.xbound-x.ge.0d0)then info=1 FLAG=0 elseif(h.ge.0d0.and.xbound-x.le.0d0)then info=1 FLAG=0 endif enddo if(key.eq.1)then write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x info=-9 endif return end subroutine irk64 subroutine dirk6_sub(grk,x,h,N,y,dy) implicit none integer,intent(in)::N double precision,intent(in)::h double precision,intent(inout)::x,y(1:N) double precision,intent(out)::dy(1:N) double precision,external::grk ! ! Gauss-Legendre 3 stage 6 order ! integer::i,j,l,p,q,Ns,m,nn,info double precision::tx,sx,sk,gxy0,dfy integer,parameter::s=3 integer,parameter::itrmax=20 double precision,parameter::delta=2d-8 double precision::a(1:s,1:s),b6(1:s),c(1:s) double precision,parameter::s15=sqrt(15d0) double precision,parameter::c536=5d0/36d0 integer,allocatable::ipiv(:) double precision,allocatable::Jmat(:,:),bvec(:,:) double precision,allocatable::K(:),ty(:) Ns=N*s c(1:3)=(/0.5d0-0.1d0*s15, 0.5d0, 0.5d0+0.1d0*s15/) a(1,1:3)=(/c536, 2d0/9d0-s15/15d0, c536-s15/30d0/) a(2,1:3)=(/c536+s15/24d0, 2d0/9d0, c536-s15/24d0/) a(3,1:3)=(/c536+s15/30d0, 2d0/9d0+s15/15d0, c536/) b6(1:3)=(/5d0/18d0, 4d0/9d0, 5d0/18d0/) allocate(ipiv(1:Ns),Jmat(1:Ns,1:Ns),bvec(1:Ns,1:1)) allocate(K(1:Ns),ty(1:N)) ! Initial guess of K do p=1,N K(1+s*(p-1))=grk(N,x,y,p) enddo do i=2,s do p=1,N K(i+s*(p-1))=K(1+s*(p-1)) enddo enddo ! ! construct J matrix do nn=1,N gxy0 = grk(N,x,y,nn) do p=1,N ty(1:N) = y(1:N) ty(p) = ty(p)+delta dfy = (grk(N,x,ty,nn)-gxy0) / delta do q=1,s do m=1,s Jmat(m+s*(nn-1),q+s*(p-1)) = -h*a(m,q)*dfy enddo enddo enddo enddo do i=1,Ns Jmat(i,i)=1d0+Jmat(i,i) enddo call dgetrf(Ns,Ns,Jmat,Ns,ipiv,info) do l=1,itrmax do i=1,s tx = x+c(i)*h ty(1:N) = y(1:N) do q=1,N do j=1,s ty(q) = ty(q)+a(i,j)*K(j+s*(q-1)) enddo enddo do p=1,N bvec(i+s*(p-1),1) = K(i+s*(p-1))-h*grk(N,tx,ty,p) enddo enddo call dgetrs('N',Ns,1,Jmat,Ns,ipiv,bvec,Ns,info) do i=1,Ns K(i) = K(i)-bvec(i,1) enddo sx=0d0 sk=0d0 do i=1,Ns sx = sx+abs(bvec(i,1)) sk = sk+abs(K(i)) enddo if(sx/sk.lt.N*1d-15)exit enddo x=x+h dy(1:N)=0d0 do p=1,N do i=1,s dy(p)=dy(p)+b6(i)*K(i+s*(p-1)) enddo y(p)=y(p)+dy(p) enddo return end subroutine dirk6_sub subroutine dirkL3C4_sub(grk,x,h,N,y,dy) implicit none integer,intent(in)::N double precision,intent(in)::h double precision,intent(inout)::x,y(1:N) double precision,intent(out)::dy(1:N) double precision,external::grk ! ! Gauss-Lobatto 3 stage 4 order ! integer::i,j,l,p,q,Ns,m,nn,info double precision::tx,sx,sk,gxy0,dfy integer,parameter::s=3 integer,parameter::itrmax=20 double precision,parameter::delta=2d-8 double precision::a(1:s,1:s),b4(1:s),c(1:s) double precision,parameter::c16=1d0/6d0 integer,allocatable::ipiv(:) double precision,allocatable::Jmat(:,:),bvec(:,:) double precision,allocatable::K(:),ty(:) Ns=N*s c(1:3)=(/0d0,0.5d0,1d0/) a(1,1:3)=(/c16, -1d0/3d0, c16/) a(2,1:3)=(/c16, 5d0/12d0, -1d0/12d0/) a(3,1:3)=(/c16, 2d0/3d0, c16/) b4(1:3)=(/c16, 2d0/3d0, c16/) allocate(ipiv(1:Ns),Jmat(1:Ns,1:Ns),bvec(1:Ns,1:1)) allocate(K(1:Ns),ty(1:N)) ! Initial guess of K do p=1,N K(1+s*(p-1))=grk(N,x,y,p) enddo do i=2,s do p=1,N K(i+s*(p-1))=K(1+s*(p-1)) enddo enddo ! ! construct J matrix do nn=1,N gxy0 = grk(N,x,y,nn) do p=1,N ty(1:N) = y(1:N) ty(p) = ty(p)+delta dfy = (grk(N,x,ty,nn)-gxy0) / delta do q=1,s do m=1,s Jmat(m+s*(nn-1),q+s*(p-1)) = -h*a(m,q)*dfy enddo enddo enddo enddo do i=1,Ns Jmat(i,i)=1d0+Jmat(i,i) enddo call dgetrf(Ns,Ns,Jmat,Ns,ipiv,info) do l=1,itrmax do i=1,s tx = x+c(i)*h ty(1:N) = y(1:N) do q=1,N do j=1,s ty(q) = ty(q)+a(i,j)*K(j+s*(q-1)) enddo enddo do p=1,N bvec(i+s*(p-1),1) = K(i+s*(p-1))-h*grk(N,tx,ty,p) enddo enddo call dgetrs('N',Ns,1,Jmat,Ns,ipiv,bvec,Ns,info) do i=1,Ns K(i) = K(i)-bvec(i,1) enddo sx=0d0 sk=0d0 do i=1,Ns sx = sx+abs(bvec(i,1)) sk = sk+abs(K(i)) enddo if(sx/sk.lt.N*1d-15)exit enddo x=x+h dy(1:N)=0d0 do p=1,N do i=1,s dy(p)=dy(p)+b4(i)*K(i+s*(p-1)) enddo y(p)=y(p)+dy(p) enddo return end subroutine dirkL3C4_sub