program main implicit none integer::i,N,info,Nt real*16::h,tol real*16::ta,tb,tx real*16,allocatable::y(:),t(:) real*16,external::grk !$call mkl_set_num_threads(1) !$call omp_set_num_threads(2) N=2 ! Number of equations allocate(y(1:N)) Nt=1001 ! Time grids allocate(t(1:Nt)) ta=0q0 ! Start time tb=50q0 ! End time do i=1,Nt t(i)=(i-1)*(tb-ta)/dble(Nt-1)+ta enddo !initial conditions y(1 )=1q0 ! y (ta) y(2 )=0q0 ! y'(ta) write(10,'(3e36.20e3)')t(1),y(1),y(2) tol=1q-15 do i=2,Nt info=-1 tx=t(i-1) h=t(i)-tx do while(info.le.0) call qirk64(grk,tx,h,N,y,t(i),info,tol) enddo write(10,'(3e36.20e3)')t(i),y(1),y(2) enddo stop end program main function grk(N,x,y,s) implicit none integer,intent(in)::N,s real*16,intent(in)::x real*16,intent(in)::y(1:N) real*16::grk grk=0q0 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 qirk64(grk,x,h,N,y,xbound,info,tol) implicit none integer,intent(in)::N real*16,intent(in)::xbound,tol real*16,intent(inout)::x,h,y(1:N) integer,intent(inout)::info real*16,external::grk real*16,allocatable::y4(:),y6(:) real*16,allocatable::dy4(:),dy6(:) integer::i,FLAG,key real*16::R0,R,delta,tx,tx0,Sy,err real*16,parameter::hmin=1q-25,hmax=0.1q0 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(1q0,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=0q0 do i=1,N Sy=Sy+abs(y(i)) enddo if(Sy.ge.1q0)then err=tol*Sy else err=tol endif do while(FLAG.eq.1) R=0q0 R0=0q0 tx=x y6(1:N)=y(1:N) tx0=x y4(1:N)=y(1:N) !$omp parallel !$omp sections !$omp section call qirk6_sub(grk,tx,h,N,y6,dy6) !$omp section call qirkL3C4_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.1q-40)then delta=(err/(2q0*R))**0.20q0 else delta=4q0 endif !step 7 if(delta.le.0.1q0)then !function changes dramatically. h=0.1q0*h elseif(delta.ge.4q0)then !function changes loosely. h=4q0*h else !function changes moderately. h=delta*h endif !step 8 if(abs(h).ge.hmax)then h=sign(1q0,h)*hmax elseif(abs(h).lt.hmin)then h=sign(1q0,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.0q0.and.xbound-x.ge.0q0)then info=1 FLAG=0 elseif(h.ge.0q0.and.xbound-x.le.0q0)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 qirk64 subroutine qirk6_sub(grk,x,h,N,y,dy) implicit none integer,intent(in)::N real*16,intent(in)::h real*16,intent(inout)::x,y(1:N) real*16,intent(out)::dy(1:N) real*16,external::grk ! ! Gauss-Legendre 3 stage 6 order ! integer::i,j,l,p,q,Ns,m,nn,info real*16::tx,sx,sk,gxy0,dfy integer,parameter::s=3 integer,parameter::itrmax=20 real*16,parameter::delta=2q-16 real*16::a(1:s,1:s),b6(1:s),c(1:s) real*16,allocatable::Jmat(:,:),bvec(:) real*16,allocatable::K(:),ty(:) integer::ISW real*16,parameter::eps=-1q0 integer,allocatable::ipivp(:),ipivq(:) real*16,allocatable::vp(:),vs(:),vw(:) real*16,parameter::s15=sqrt(15q0) real*16,parameter::c536=5q0/36q0 Ns=N*s c(1:3)=(/0.5q0-0.1q0*s15, 0.5q0, 0.5q0+0.1q0*s15/) a(1,1:3)=(/c536, 2q0/9q0-s15/15q0, c536-s15/30q0/) a(2,1:3)=(/c536+s15/24q0, 2q0/9q0, c536-s15/24q0/) a(3,1:3)=(/c536+s15/30q0, 2q0/9q0+s15/15q0, c536/) b6(1:3)=(/5q0/18q0, 4q0/9q0, 5q0/18q0/) allocate(Jmat(1:Ns,1:Ns),bvec(1:Ns)) allocate(K(1:Ns),ty(1:N)) allocate(ipivp(1:Ns),ipivq(1:Ns),vp(1:Ns),vs(1:Ns),vw(1:Ns)) ! 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)=1q0+Jmat(i,i) enddo ! call dgetrf(Ns,Ns,Jmat,Ns,ipiv,info) ISW=0 call GECP_Q(Jmat,Ns,Ns,bvec,eps,ISW,ipivp,ipivq,vp,vs,vw,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)) = K(i+s*(p-1))-h*grk(N,tx,ty,p) enddo enddo !call dgetrs('N',Ns,1,Jmat,Ns,ipiv,bvec,Ns,info) ISW=2 call GECP_Q(Jmat,Ns,Ns,bvec,eps,ISW,ipivp,ipivq,vp,vs,vw,info) do i=1,Ns K(i) = K(i)-bvec(i) enddo sx=0q0 sk=0q0 do i=1,Ns sx = sx+abs(bvec(i)) sk = sk+abs(K(i)) enddo if(sx/sk.lt.N*1q-30)exit enddo x=x+h dy(1:N)=0q0 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 qirk6_sub subroutine qirkL3C4_sub(grk,x,h,N,y,dy) implicit none integer,intent(in)::N real*16,intent(in)::h real*16,intent(inout)::x,y(1:N) real*16,intent(out)::dy(1:N) real*16,external::grk ! ! Gauss-Lobatto 3 stage 4 order ! integer::i,j,l,p,q,Ns,m,nn,info real*16::tx,sx,sk,gxy0,dfy integer,parameter::s=3 integer,parameter::itrmax=20 real*16,parameter::delta=2q-16 real*16::a(1:s,1:s),b4(1:s),c(1:s) real*16,allocatable::Jmat(:,:),bvec(:) real*16,allocatable::K(:),ty(:) real*16,parameter::c16=1q0/6q0 integer::ISW real*16,parameter::eps=-1q0 integer,allocatable::ipivp(:),ipivq(:) real*16,allocatable::vp(:),vs(:),vw(:) Ns=N*s c(1:3)=(/0q0,0.5q0,1q0/) a(1,1:3)=(/c16, -1q0/3q0, c16/) a(2,1:3)=(/c16, 5q0/12q0, -1q0/12q0/) a(3,1:3)=(/c16, 2q0/3q0, c16/) b4(1:3)=(/c16, 2q0/3q0, c16/) allocate(Jmat(1:Ns,1:Ns),bvec(1:Ns)) allocate(K(1:Ns),ty(1:N)) allocate(ipivp(1:Ns),ipivq(1:Ns),vp(1:Ns),vs(1:Ns),vw(1:Ns)) ! 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)=1q0+Jmat(i,i) enddo ISW=0 call GECP_Q(Jmat,Ns,Ns,bvec,eps,ISW,ipivp,ipivq,vp,vs,vw,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)) = K(i+s*(p-1))-h*grk(N,tx,ty,p) enddo enddo ISW=2 call GECP_Q(Jmat,Ns,Ns,bvec,eps,ISW,ipivp,ipivq,vp,vs,vw,info) do i=1,Ns K(i) = K(i)-bvec(i) enddo sx=0q0 sk=0q0 do i=1,Ns sx = sx+abs(bvec(i)) sk = sk+abs(K(i)) enddo if(sx/sk.lt.N*1q-30)exit enddo x=x+h dy(1:N)=0q0 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 qirkL3C4_sub ! ==================================================================== ! ! -------------------------------------------------------------------- ! ! ! SUBROUTINE GECP_Q(A,NX,N,B,EPS,ISW,P,Q,VP,VS,VW,ICON) ! ! ! ! -------------------------------------------------------------------- ! ! <> ! ! MAY 31ST, 2002 ! ! <> ! ! WATANABE, YOSHITAKA ! ! COMPUTING AND COMMUNICATIONS CENTER, KYUSHU UNIERSITY ! ! <> -1 ! ! SUBROUTINE `GECP_Q' SOLVES THE REAL LINEAR SYSTEM X=A B ! ! BY GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING AND SCALING. ! ! <> ! ! QUADRUPLE ! ! <> ! ! A: (INPUT) REAL COEFFICIENT MARIX. ! ! (OUTPUT) LU DECOMPOSED MATRIX (SEE REMARK). ! ! NX: (INPUT) AJUSTED SIZE OF MATRIX A. ! ! N: (INPUT) DIMENSION OF THE LINEAR SYSTEM (NX>=N>1). ! ! B: (INPUT) RIGHT-HAND-SIDE VECTOR ! ! (OUTPUT) NUMERICAL APPLOXIMATE SOLUTION. ! ! EPS: (INPUT) PARAMETER FOR CHECKING THE SINGULALITY (>=0). ! ! IF EPS<=0 THEN EPS IS SET BY MACHINE EPSILON. ! ! ISW: (INPUT) ISW=0 ...EXECUTE ONLY LU DECOMPOSITION. ! ! ISW=1 ...SOLVE LINEAR SYSTEM A*X=B USUALLY. ! ! ISW=2 ...SKIP LU DECOMPOSITION AND SOLVE A*X=B. ! ! ISW=3 ...SOLVE LINEAR SYSTEM (A^T)*X=B. ! ! ISW=4 ...SKIP LU DECOMPOSITION AND SOLVE A(^T)*X=B. ! ! ISW=OTHERS .... EXECUTE THE SAME AS ISW=1. ! ! P: (OUTPUT) PIVOTING INFORMATION VECTOR FOR COLUMN. ! ! Q: (OUTPUT) PIVOTING INFORNATION VECTOR FOR ROW. ! ! VP: (OUTPUT) SCALING INFORMATION VECTOR FOR NORMALIZATION. ! ! VS: (OUTPUT) SCALING INFORMATION VECTOR OF EQUATIONS. ! ! VW: (OUT) WORKING VECTOR. ! ! ICON: (OUTPUT) CONDITION CODE. ! ! 0 ... PROGRAM ENDED NORMALY. ! ! 2000 ... ALL ELEMENTS OF A COLUMN OR ROW ARE ZERO ! ! OR PIVOT LESS THAN EPS. ! ! 3000 ... PARAMETER ERROR. ! ! <> ! ! A IS DECOMPOSED LU WHERE ! ! L ... THE LOWER TRIANGULAR MATRIX WHOSE ALL DIAGONAL ELEMENT ! ! ARE 1. ! ! U ... THE UPPER TRIANGULAR MATRIX WHOSE DIAGONAL ELEMENTS ARE ! ! STORED AS 1/U(I,I). ! ! -------------------------------------------------------------------- ! ! ! ! ==================================================================== ! ! DECRALATION ! ! -------------------------------------------------------------------- ! IMPLICIT NONE ! ! INTEGER,INTENT(IN) :: N,NX,ISW REAL(KIND(1.0Q0)),DIMENSION(NX,NX),INTENT(INOUT) :: A REAL(KIND(1.0Q0)),DIMENSION(NX),INTENT(INOUT) :: B,VS,VP REAL(KIND(1.0Q0)),DIMENSION(NX),INTENT(OUT) :: VW REAL(KIND(1.0Q0)),INTENT(IN) :: EPS INTEGER,DIMENSION(NX),INTENT(INOUT) :: P,Q INTEGER,INTENT(OUT) :: ICON ! ! REAL(KIND(1.0Q0)) :: D,EPSZ INTEGER :: I,J,K,L,C ! ! INTRINSIC EPSILON,MAX,ABS ! ! ! ==================================================================== ! ! CHECKING PARAMETERS ! ! -------------------------------------------------------------------- ! ICON=0 IF(EPS<=0.0Q0) THEN EPSZ=EPSILON(1.0Q0) ELSE EPSZ=EPS END IF IF(N<1.OR.N>NX) THEN ICON=3000 RETURN END IF IF(ICON/=0.AND.ICON/=2.AND.ICON/=3.AND.ICON/=4) ICON=1 ! ! ! ==================================================================== ! ! NORMALIZATION ! ! -------------------------------------------------------------------- ! IF(ISW/=2.AND.ISW/=4) THEN ! ! ---------- FOR EACH EQUATION --------------------------------------- ! DO I=1,N D=0.0Q0 DO J=1,N D=MAX(ABS(A(I,J)),D) END DO IF(D==0.0Q0) THEN ICON=2000 RETURN END IF DO J=1,N A(I,J)=A(I,J)/D END DO VP(I)=D END DO ! ---------- FOR EACH UNKNOWN VALUE----------------------------------- ! DO J=1,N D=0.0Q0 DO I=1,N D=MAX(ABS(A(I,J)),D) END DO IF(D==0.0Q0) THEN ICON=2000 RETURN END IF DO I=1,N A(I,J)=A(I,J)/D END DO VS(J)=D END DO ! ==================================================================== ! ! COMPLETE PIVOTING ! ! -------------------------------------------------------------------- ! DO K=1,N P(K)=K ! INITIAL VALUE Q(K)=K END DO ! ---------- SERCH PIVOT --------------------------------------------- ! DO K=1,N D=0.0Q0 DO J=K,N DO I=K,N IF(ABS(A(I,J))>D) THEN D=ABS(A(I,J)) L=I ! COLUMN NUMBER C=J ! ROW NUMBER END IF END DO END DO IF(D==0.0Q0) THEN ICON=2000 RETURN END IF ! ---------- CHANGE INDEX -------------------------------------------- ! IF(K/=L) THEN DO J=1,N D=A(K,J) A(K,J)=A(L,J) A(L,J)=D END DO J=P(K) P(K)=P(L) P(L)=J END IF ! ! IF(K/=C) THEN DO I=1,N D=A(I,K) A(I,K)=A(I,C) A(I,C)=D END DO J=Q(K) Q(K)=Q(C) Q(C)=J END IF ! ! ! ==================================================================== ! ! LU DECOMPOSITION ! ! -------------------------------------------------------------------- ! IF(ABS(A(K,K))