module GBL implicit none integer::count end module GBL program main use GBL implicit none integer::N,Ns,istep,ih0,Jup,info,ih,iJac real*16::h,h0,err0,eta0,tol real*16,allocatable::y(:),Rtol(:),Atol(:) real*16,allocatable::Jmat(:,:),errJ(:,:),z0(:),Jac(:,:) integer,allocatable::ipivp(:),ipivq(:),epivp(:),epivq(:) real*16,allocatable::vp(:),vs(:),evp(:),evs(:) integer::i,Nx real*16::xa,xb,tx real*16,allocatable::x(:) external::grk N=2 ! Number of 1st order ODEs allocate(y(1:N),Rtol(1:N),Atol(1:N)) y=0q0; Rtol=0q0; Atol=0q0 Nx=2 allocate(x(1:Nx)) x=0q0 xa=0q0 ! Initial value of x xb=2q0 ! End value of x do i=1,Nx ! Separate equal interval x. x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa enddo y(1)=2q0 ! Initial values at x of y_1 y(2)=-0.66q0 ! Initial values at x of y_2 tol=1q-4 ! Tolerance h=1q-6 ! Initial step size Rtol(1:N)=tol ! Relative tolerance Atol(1:N)=tol ! Absolute tolerance !------------ Initial set up ------------ Ns=N*3 ! Ns=N*s, s means s-stage IRK. allocate(Jmat(1:Ns,1:Ns),errJ(1:N,1:N),Jac(1:N,1:N),z0(1:Ns)) allocate(ipivp(1:Ns),ipivq(1:Ns),epivp(1:N),epivq(1:N)) allocate(vp(1:Ns),vs(1:Ns),evp(1:N),evs(1:N)) ipivp=0; ipivq=0; epivp=0; epivq=0 vp=0q0; vs=0q0; evp=0q0; evs=0q0 Jmat=0q0; errJ=0q0; Jup=0; iJac=0; Jac=0q0 ih0=0; z0=0q0; eta0=0q0; h0=0q0; err0=0q0 !---------------------------------------- count=0 istep=0 do i=2,Nx info=0 tx=x(i-1) do while(info.le.0) call irkgl(istep,grk,N,tx,h,x(i),y,info,Atol,Rtol & ,ih,z0,ih0,h0,eta0,err0,Jup & ,ipivp,ipivq,vp,vs,Jmat,epivp,epivq,evp,evs,errJ,iJac,Jac) write(20,'(4e25.15e3)')tx,y(1),y(2),h enddo h=h0 write(30,'(4e25.15e3)')x(i),y(1),y(2),h enddo write(6,*)count stop end program main subroutine grk(N,x,y,f) use GBL implicit none integer,intent(in)::N real*16,intent(in)::x,y(1:N) real*16,intent(out)::f(1:N) ! Write right-hand-side of 1st order ODEs f(1)=y(2) f(2)=1q6*((1q0-y(1)**2)*y(2)-y(1)) count=count+1 return end subroutine grk !============================= subroutine irkgl(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol & ,ih,z0,ih0,h0,eta0,err0,Jup & ,ipivp,ipivq,vp,vs,Jmat,epivp,epivq,evp,evs,errJ,iJac,Jac) implicit none integer,parameter::s=3 integer,intent(in)::Neq integer,intent(inout)::istep,info,Jup,iJac real*16,intent(in)::xend,Atol(1:Neq),Rtol(1:Neq) real*16,intent(inout)::x,h,y(1:Neq) integer,intent(inout)::ih0,ih real*16,intent(inout)::z0(1:Neq*s),h0,eta0,err0 integer,intent(inout)::ipivp(1:Neq*s),ipivq(1:Neq*s),epivp(1:Neq),epivq(1:Neq) real*16,intent(inout)::vp(1:Neq*s),vs(1:Neq*s),evp(1:Neq),evs(1:Neq) real*16,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq) real*16,intent(inout)::Jac(1:Neq,1:Neq) external::grk ! ! Implicit Runge-Kutta method based on ! the Gauss-Legendre 3-stage 6-order ! ! Properties of this routine: ! 1. A-stable ! 2. Symplectic ! 3. Symmetric ! 4. Step size control ! Note, Gauss-Legendre IRK method is Symplectic ! even if we change the step size. ! ! Meaning of parameters ! istep : Number of IRK step ! grk : Right hand Side of ODEs ! Neq : Number of 1st-order ODEs ! x : Integral parameter (automatically updated) ! h : Step size (automatically updated) ! xend : End point of the x range ! y : Values of ODEs ! info : Information of the IRK process ! Atol : Absolute tolerance ! Rtol : Relative tolerance ! ! Other parameters are work parameters, ! referenced for istep >= 1 ! ***0 : Previous parameters ! Jmat : LU decomposited J' (= I-hAJ) matrix ! ipiv : Pivot information for Jmat ! errJ : LU decomposited (I-h\gamma J) matrix ! epiv : Pivot information for errJ ! Jup : Decide the update of Jmat and errJ, ! Jup = 0 --> No update ! Jup = 1 --> Update ! Jac : Jacobian matrix of the grk function ! iJac : Decide the update of Jac, ! iJac = 0 --> No update when Jup=1 ! iJac = 1 --> Update when Jup=1 ! ! How to use? ! 1. Call this routine with istep = 0 and info = 0. ! 2. Call and Loop this routine while info <= 0. ! ! ==== Example ===== ! istep=0 ! info=0 ! do while(info.le.0) ! call irk(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol & ! ,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac) ! enddo ! ================== ! ! After starting computation with istep=0, ! you must not touch WORK parameters. ! ! ! istep = 0 : when you start computation, ! set WORK parameters like ; ! ih = 0 ! z0(1:Ns) = 0d0 ! ih0 = 0 ! h0 = 0d0 ! eta0 = 0d0 ! err0 = 0d0 ! Jup = 0 ! ipiv(1:Ns) = 0 ! Jmat(1:Ns,1:Ns) = 0d0 ! epiv(1:Neq) = 0 ! errJ(1:Neq,1:Neq) = 0d0 ! iJac = 0 ! Jac(1:Neq,1:Neq) = 0d0 ! ! Author : sikino ! URL : http://slpr.sakura.ne.jp/qp/ ! Date : 2019/01/14 (yyyy/mm/dd) ! 2019/01/21 keep Jacobian matrix ! 2019/01/22 did Something ! 2019/01/23 initial value estimation ! real*16::tx,th real*16,allocatable::ty(:),tz0(:) integer,parameter::kmax=7 ! Newton iteration limit real*16,parameter::hmin=1q-20,hmax=1q0 integer::kexit,Ns,key,FLAG,Newt real*16::theta,err,fac,del,del1,del2,teta0,th0 if(istep.le.-1)then write(6,*)"**** Error, unexpected istep" stop endif Ns=Neq*s if(abs(h).ge.hmax)then h=sign(1q0,h)*hmax endif if(h.ge.abs(xend-x))h=xend-x FLAG=1 if(abs(x-xend).le.hmin)then info=1 FLAG=0 endif if(istep.eq.0)then ih=0; z0=0q0; ih0=0; h0=0q0; eta0=0q0; err0=0q0 Jup=0; ipivp=0; ipivq=0; Jmat=0q0; epivp=0; epivq=0; errJ=0q0 vp=0q0; vs=0q0; evp=0q0; evs=0q0; iJac=0; Jac=0q0 call discrete_h(h,ih,th,hmin,hmax) h=th ih0 = ih+1 Jup = 1 iJac = 1 else endif key=0 allocate(ty(1:Neq),tz0(1:Ns)) ty=0q0; tz0=0q0 do while(FLAG.eq.1) if(ih.ne.ih0)then Jup=1 endif tx=x ty(1:Neq)=y(1:Neq) tz0(1:Ns)=z0(1:Ns) teta0 = eta0 th0=h0 call qirk6(istep,grk,Neq,tx,h,ty,Jup & ,Atol,Rtol,err,kmax,kexit,Newt,theta & ,tz0,th0,teta0,ipivp,ipivq,vp,vs,Jmat,epivp,epivq,evp,evs,errJ,iJac,Jac) ! Even if the step is fail, ipiv,Jmat,epiv,errJ are updated if Jup=1. if(Jup.eq.1)then Jup=0 endif if(kexit.eq.1)then ! Change step size !fac = 0.9q0*(2q0*kmax+1q0)/(2q0*kmax+dble(Newt-1)) fac = 0.95q0*(2q0*kmax+1q0)/(2q0*kmax+dble(Newt)) if(err.ge.1q-60)then if(istep.eq.0)then del = fac*((1q0/err)**(0.33q0)) !(27) else del1 = fac*((1q0/err)**(0.33q0)) !(27) del2 = del1*(h/h0)*((err0/err)**(0.33q0)) !(27) del = del1 if(del2.lt.del)del=del2 if(del.gt.1q0)then del=del1 else del=del2 endif endif else del=100q0 endif elseif(kexit.eq.2)then del=0q0 else write(6,*)" **** detect unexpected kexit" stop endif ! Accept or Reject if(err.le.1q0.or.key.eq.1)then FLAG=0 ! This step with h is accepted x=x+h y(1:Neq)=ty(1:Neq) z0(1:Ns)=tz0(1:Ns) h0=h ih0=ih eta0=teta0 err0=err iJac = 1 Jup=1 ! Don't update Jacobian for next step if(Newt.le.2.or.theta.lt.1q-3)Jup=0 !if(Newt.le.1.or.theta.lt.1q-3)Jup=0 endif if(del.le.0.1q0)then !function changes dramatically. h=0.1q0*h elseif(del.ge.4q0)then !function changes loosely. h=4q0*h else !function changes moderately. h=del*h endif 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 size alignment if(abs(xend-x).le.abs(h))then h=xend-x Jup=1 if(abs(h).le.hmin)then info=1 FLAG=0 endif else call discrete_h(h,ih,th,hmin,hmax) h=th endif if(h.le.0q0.and.xend-x.ge.0q0)then info=1 FLAG=0 elseif(h.ge.0q0.and.xend-x.le.0q0)then info=1 FLAG=0 endif if(key.eq.1)then write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x info=-9 endif enddo istep=istep+1 return end subroutine irkgl subroutine qirk6(istep,grk,Neq,x,h,y,Jup & ,Atol,Rtol,err,kmax,kexit,Newt,theta & ,z0,h0,eta0,ipivp,ipivq,vp,vs,Jmat,epivp,epivq,evp,evs,errJ,iJac,Jac) implicit none integer,parameter::s=3 integer,intent(in)::istep,Neq,Jup,kmax integer,intent(out)::kexit,Newt real*16,intent(in)::h,Atol(1:Neq),Rtol(1:Neq) real*16,intent(inout)::x,y(1:Neq),z0(1:Neq*s),h0,eta0 real*16,intent(out)::err,theta integer,intent(inout)::ipivp(1:Neq*s),ipivq(1:Neq*s),epivp(1:Neq),epivq(1:Neq) integer,intent(inout)::iJac real*16,intent(inout)::vp(1:Neq*s),vs(1:Neq*s),evp(1:Neq),evs(1:Neq) real*16,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq) real*16,intent(inout)::Jac(1:Neq,1:Neq) external::grk ! ! istep >= 0 ! ! Input ! Jup = 0 : Don't update Jacobian ! = 1 : Update Jacobian ! Output ! kexit = 1 : Newton iteration converge. ! = 2 : Newton iteration didn't converge. ! ! Below parameters have meaning when kexit = 1. ! Newt : Number of Newton iteration till converge. ! theta : Convergion rate, \theta < 1. ! err : Estimated err, if err < 1, satisfied required tol. ! integer::i,j,k,n,m,p,q,Ns,info,ISW real*16,parameter::kappa=5q-2 real*16,parameter::Uround=5q-32 real*16,parameter::sq15=sqrt(15q0) real*16,parameter::zero=0q0 ! Real eigenvalue of A matrix in butcher table for Gauss-Legendre real*16,parameter::gamma=0.215314423116112178244733530380696q0 real*16::a(1:s,1:s),c(1:s),d(1:s),d2(1:s),dc(1:s) real*16::c12,c23,c31,xc1,xc2,xc3,xx,omega real*16::x0,tx real*16,allocatable::z(:),y0(:),dy(:),f(:),tf(:),ty(:),tf0(:) real*16,allocatable::w(:),ew0(:),e(:),vw(:),ew(:),evw(:) real*16::Ntol,sc,sdz,sdz0,tmp,eta Ns=Neq*s allocate(z(1:Ns)) z=0q0 ! 3-stage Gauss-Legendre c(1:3)=(/0.5q0-0.1q0*sq15, 0.5q0, 0.5q0+0.1q0*sq15/) a(1,1:3)=(/5q0/36q0, 2q0/9q0-sq15/15q0, 5q0/36q0-sq15/30q0/) a(2,1:3)=(/5q0/36q0+sq15/24q0, 2q0/9q0, 5q0/36q0-sq15/24q0/) a(3,1:3)=(/5q0/36q0+sq15/30q0, 2q0/9q0+sq15/15q0, 5q0/36q0/) d(1:3)=(/5q0/3q0,-4q0/3q0,5q0/3q0/) d2(1:3)=(/-15q0-10q0*sqrt(5q0/3q0),12q0,-15q0+10q0*sqrt(5q0/3q0)/) dc(1:3)=(/(5q0+sq15)*10q0/3q0, -40q0/3q0 ,(5q0-sq15)*10q0/3q0/) if(istep.eq.0)then z(1:Ns)=0q0 else allocate(dy(1:Neq)) dy=0q0 do n=1,Neq do j=1,s dy(n)=dy(n)+d(j)*z0((n-1)*s+j) enddo enddo ! Lagrange interpolation omega=h/h0 c12=1q0/(c(1)-c(2)) c23=1q0/(c(2)-c(3)) c31=1q0/(c(3)-c(1)) do n=1,Neq do p=1,s xx=1q0+omega*c(p) xc1=xx-c(1) xc2=xx-c(2) xc3=xx-c(3) z((n-1)*s+p)=& -z0((n-1)*s+1)*xc2*xc3*(c12*c31)*xx/c(1) & -z0((n-1)*s+2)*xc3*xc1*(c12*c23)*xx/c(2) & -z0((n-1)*s+3)*xc1*xc2*(c31*c23)*xx/c(3) & -dy(n) enddo enddo deallocate(dy) endif allocate(w(1:Ns),vw(1:Ns),ew(1:Neq),evw(1:Neq)) w=0q0; vw=0q0; ew=0q0; evw=0q0 if(Jup.eq.1)then ! LU decomposition of J' matrix if(iJac.eq.1)then call Jacobian(Neq,x,y,grk,Jac) iJac=0 endif do m=1,Neq do n=1,Neq do q=1,s do p=1,s Jmat(s*(n-1)+p,s*(m-1)+q) = - h*Jac(n,m)*a(p,q) enddo enddo enddo enddo do n=1,Neq do m=1,Neq errJ(n,m) = - h*gamma*Jac(n,m) enddo enddo do i=1,Ns Jmat(i,i) = 1q0 + Jmat(i,i) enddo do i=1,Neq errJ(i,i) = 1q0 + errJ(i,i) enddo ! LU factorization for main part of IRK ISW=0 call GECP_Q(Jmat,Ns,Ns,w,zero,ISW,ipivp,ipivq,vp,vs,vw,info) ! LU factorization for estimate error ISW=0 call GECP_Q(errJ,Neq,Neq,ew,zero,ISW,epivp,epivq,evp,evs,evw,info) endif allocate(f(1:Ns),tf(1:Neq),ty(1:Neq)) f=0q0; tf=0q0; ty=0q0 !=============================== Ntol=sqrt(Rtol(1)) !if(0.03q0.le.Ntol)Ntol=0.03q0 !if(1q-6.le.Ntol)Ntol=1q-6 !if(1q-9.le.Ntol)Ntol=1q-9 ! <-- double precision if(1q-12.le.Ntol)Ntol=1q-12 !if(1q-15.le.Ntol)Ntol=1q-15 !if(1q-20.le.Ntol)Ntol=1q-20 !=============================== sdz0=0q0 ! Initialize ! Simple Newton iteration kexit=0 do k=1,kmax do j=1,s tx = x + c(j)*h do n=1,Neq ty(n) = y(n) + z((n-1)*s+j) enddo call grk(Neq,tx,ty,tf) do n=1,Neq f((n-1)*s+j) = tf(n) enddo enddo w(1:Ns)=0q0 do n=1,Neq do p=1,s do j=1,s w((n-1)*s+p) = w((n-1)*s+p) + a(p,j)*f((n-1)*s+j) enddo w((n-1)*s+p) = z((n-1)*s+p) - h*w((n-1)*s+p) enddo enddo ! Solve J' \delta z = - w do i=1,Ns w(i) = - w(i) enddo ISW=2 call GECP_Q(Jmat,Ns,Ns,w,zero,ISW,ipivp,ipivq,vp,vs,vw,info) ! --> Now, w is \Delta z ! z^{(k+1)} = z^{(k)} + Delta z do i=1,Ns z(i) = z(i) + w(i) enddo sdz=0q0 do i=1,Ns sdz = sdz + w(i)**2 enddo sdz=sqrt(sdz) ! Stop iteration criteria if(istep.eq.0.and.k.eq.1)then ! Do nothing kexit=0 elseif(istep.eq.0.and.k.ge.2)then if(sdz0-sdz.lt.0q0)then ! Convergion rate > 1, must change small step size. kexit=2 elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then ! Rough convergion estimation fail, must change small step size. kexit=2 elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then ! Good behavior. Iteration finish. kexit=1 endif elseif(istep.ge.1.and.k.eq.1)then tmp = Uround if(eta0.gt.tmp)tmp=eta0 tmp=tmp**0.8q0 if(tmp*sdz.lt.kappa*Ntol)then kexit=1 endif elseif(istep.ge.1.and.k.ge.2)then if(sdz0-sdz.lt.0q0)then ! Convergion rate > 1, must change small step size. kexit=2 elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then ! Rough convergion estimation fail, must change small step size. kexit=2 elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then ! Good behavior. Iteration finish. kexit=1 endif else write(6,*)" *****Unexpected parameters" stop endif if(kexit.ne.0)exit sdz0 = sdz enddo ! ! 0 < \eta < \infty --> Good. ! if \eta ~ 0, good behavior ! -\infty < \eta < -1, --> Bad. ! Error increase as iteration increase. ! if(k.eq.kmax+1.or.kexit.eq.2.or.kexit.eq.0)then ! Did not converge k_max iteration. kexit=2 eta=100q0 Newt=kmax err=100q0 theta=1000q0 ! Convergion ratio, \theta ~ 0 is Good. ! h0, z0, x, y are don't updated return else ! (sdz0-sdz) > 0 if(k.eq.1)then ! No sdz0 case. eta = Uround if(eta0.gt.eta)eta=eta0 eta = eta**0.8q0 theta = 100q0 ! Here, \theta cannot evaluate because k=1. else if(sdz0.eq.sdz)then ! z does not change --> converge enough. eta = Uround theta = 0q0 ! \theta=0. else ! General case. eta = sdz/(sdz0-sdz) theta = sdz/sdz0 endif endif Newt = k endif x0=x allocate(y0(1:Neq)) y0(1:Neq)=y(1:Neq) ! Update x and y(1:Neq) x=x+h do n=1,Neq do j=1,s y(n) = y(n) + d(j)*z((n-1)*s+j) enddo enddo ! Error estimation allocate(ew0(1:Neq),tf0(1:Neq),e(1:Neq)) ew=0q0; ew0=0q0; tf0=0q0; e=0q0 do n=1,Neq do j=1,s ew(n) = ew(n) + dc(j)*z((n-1)*s+j) enddo enddo ew0(1:Neq)=ew(1:Neq) ISW=2 call GECP_Q(errJ,Neq,Neq,ew,zero,ISW,epivp,epivq,evp,evs,evw,info) do n=1,Neq e(n)=ew(n) enddo err=0q0 do n=1,Neq sc=abs(y0(n)) if(abs(y(n)).gt.y0(n))sc=abs(y(n)) sc=Atol(n)+sc*Rtol(n) err=err+(e(n)/sc)**2 enddo err=sqrt(err/dble(Neq)) if(err.ge.1q0)then call grk(Neq,x0,y0,tf0) do n=1,Neq ty(n) = y0(n) + e(n) enddo call grk(Neq,x0,ty,tf) do n=1,Neq ew(n) = ew0(n) + gamma*h*(tf(n)-tf0(n)) enddo ISW=2 call GECP_Q(errJ,Neq,Neq,ew,zero,ISW,epivp,epivq,evp,evs,evw,info) do n=1,Neq e(n)=ew(n) enddo err=0q0 do n=1,Neq sc=abs(y0(n)) if(abs(y(n)).gt.y0(n))sc=abs(y(n)) sc=Atol(n)+sc*Rtol(n) err=err+(e(n)/sc)**2 enddo err=sqrt(err/dble(Neq)) endif z0(1:Ns)=z(1:Ns) h0=h eta0=eta return end subroutine qirk6 subroutine discrete_h(h,ih,th,hmin,hmax) implicit none real*16,intent(in)::h,hmin,hmax integer,intent(out)::ih real*16,intent(out)::th ! ! Discrete h to th with the index ih ! real*16::dx,hmin1,hmax1 integer,parameter::imax=300 hmin1=0.10q0*hmin hmax1=1q0*hmax dx=(log10(hmax1)-log10(hmin1))/dble(imax) do ih=0,imax th=10q0**(-ih*dx+log10(hmax1)) if(th.le.abs(h))exit enddo if(h.lt.0q0)th=-th return end subroutine discrete_h subroutine Jacobian(Neq,x,y,fxy,Jac) implicit none integer,intent(in)::Neq real*16,intent(in)::x,y(1:Neq) real*16,intent(out)::Jac(1:Neq,1:Neq) external::fxy ! ! Compute Jacobian matrix at (x,y_1,y_2,...y_{Neq}) ! Obtained by forward difference ! integer::n,m real*16::dy real*16,parameter::delta=4q-16 real*16,allocatable::f0(:),f1(:),ty(:) allocate(f0(1:Neq),f1(1:Neq),ty(1:Neq)) f0=0q0; f1=0q0; ty=0q0 call fxy(Neq,x,y,f0) do m=1,Neq ty(1:Neq) = y(1:Neq) dy=sqrt(abs(y(m))) if(dy.lt.1q0)dy=1q0 dy=delta*dy ty(m) = ty(m)+dy call fxy(Neq,x,ty,f1) do n=1,Neq Jac(n,m) = (f1(n)-f0(n)) / dy enddo enddo return end subroutine Jacobian ! ==================================================================== ! ! -------------------------------------------------------------------- ! ! ! 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. ! ! change by sikino at 2019/01/21 to eliminate warning ! ! <> ! ! 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 L=0; C=0 ! added by sikino ! ! ! ==================================================================== ! ! 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))