diff -C1 -r 206b/recipes_f-90/COPYRIGHT.NOTICE 208/recipes_f-90/COPYRIGHT.NOTICE *** 206b/recipes_f-90/COPYRIGHT.NOTICE Tue May 07 11:56:22 1996 --- 208/recipes_f-90/COPYRIGHT.NOTICE Sun Jul 06 08:14:49 1997 *************** *** 1,2 **** The contents of this Numerical Recipes Fortran90 distribution are ! Copyright (C) 1986-1996 Numerical Recipes Software. --- 1,2 ---- The contents of this Numerical Recipes Fortran90 distribution are ! Copyright (C) 1986-1997 Numerical Recipes Software. diff -C1 -r 206b/recipes_f-90/demo/answers/xcaldat.reslt 208/recipes_f-90/demo/answers/xcaldat.reslt *** 206b/recipes_f-90/demo/answers/xcaldat.reslt Tue May 07 11:41:48 1996 --- 208/recipes_f-90/demo/answers/xcaldat.reslt Sun Jul 06 08:11:01 1997 *************** *** 6,8 **** January 1 1 1721424 January 1 1 ! October 14 1582 2299170 October 24 1582 October 15 1582 2299161 October 15 1582 --- 6,8 ---- January 1 1 1721424 January 1 1 ! October 4 1582 2299160 October 4 1582 October 15 1582 2299161 October 15 1582 diff -C1 -r 206b/recipes_f-90/demo/answers/xjulday.reslt 208/recipes_f-90/demo/answers/xjulday.reslt *** 206b/recipes_f-90/demo/answers/xjulday.reslt Tue May 07 11:41:51 1996 --- 208/recipes_f-90/demo/answers/xjulday.reslt Sun Jul 06 08:11:19 1997 *************** *** 5,7 **** January 1 1 1721424 One day later ! October 14 1582 2299170 Day before Gregorian calendar October 15 1582 2299161 Gregorian calendar adopted --- 5,7 ---- January 1 1 1721424 One day later ! October 4 1582 2299160 Day before Gregorian calendar October 15 1582 2299161 Gregorian calendar adopted diff -C1 -r 206b/recipes_f-90/demo/data/DATES.DAT 208/recipes_f-90/demo/data/DATES.DAT *** 206b/recipes_f-90/demo/data/DATES.DAT Tue May 07 11:42:59 1996 --- 208/recipes_f-90/demo/data/DATES.DAT Tue Jul 01 05:56:22 1997 *************** *** 4,6 **** 01 01 1 One day later ! 10 14 1582 Day before Gregorian calendar 10 15 1582 Gregorian calendar adopted --- 4,6 ---- 01 01 1 One day later ! 10 04 1582 Day before Gregorian calendar 10 15 1582 Gregorian calendar adopted diff -C1 -r 206b/recipes_f-90/doc/COPYRIGHT.NOTICE 208/recipes_f-90/doc/COPYRIGHT.NOTICE *** 206b/recipes_f-90/doc/COPYRIGHT.NOTICE Tue May 07 11:56:22 1996 --- 208/recipes_f-90/doc/COPYRIGHT.NOTICE Sun Jul 06 08:14:49 1997 *************** *** 1,2 **** The contents of this Numerical Recipes Fortran90 distribution are ! Copyright (C) 1986-1996 Numerical Recipes Software. --- 1,2 ---- The contents of this Numerical Recipes Fortran90 distribution are ! Copyright (C) 1986-1997 Numerical Recipes Software. diff -C1 -r 206b/recipes_f-90/doc/VERSION 208/recipes_f-90/doc/VERSION *** 206b/recipes_f-90/doc/VERSION Tue May 07 11:55:28 1996 --- 208/recipes_f-90/doc/VERSION Sat Jul 19 09:26:11 1997 *************** *** 1,2 **** Numerical Recipes in FORTRAN90 ! Unix Version 2.06 --- 1,2 ---- Numerical Recipes in FORTRAN90 ! Unix Version 2.08 diff -C1 -r 206b/recipes_f-90/recipes/amoeba.f90 208/recipes_f-90/recipes/amoeba.f90 *** 206b/recipes_f-90/recipes/amoeba.f90 Mon May 06 13:41:53 1996 --- 208/recipes_f-90/recipes/amoeba.f90 Fri May 30 16:02:39 1997 *************** *** 16,17 **** --- 16,18 ---- INTEGER(I4B), PARAMETER :: ITMAX=5000 + REAL(SP), PARAMETER :: TINY=1.0e-10 INTEGER(I4B) :: ihi,ndim *************** *** 35,37 **** y(ihi)=ytmp ! rtol=2.0_sp*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) if (rtol < ftol) then --- 36,38 ---- y(ihi)=ytmp ! rtol=2.0_sp*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))+TINY) if (rtol < ftol) then diff -C1 -r 206b/recipes_f-90/recipes/caldat.f90 208/recipes_f-90/recipes/caldat.f90 *** 206b/recipes_f-90/recipes/caldat.f90 Mon May 06 13:41:54 1996 --- 208/recipes_f-90/recipes/caldat.f90 Fri May 30 16:02:39 1997 *************** *** 10,11 **** --- 10,13 ---- ja=julian+1+jalpha-int(0.25_sp*jalpha) + else if (julian < 0) then + ja=julian+36525*(1-julian/36525) else *************** *** 23,24 **** --- 25,27 ---- if (iyyy <= 0) iyyy=iyyy-1 + if (julian < 0) iyyy=iyyy-100*(1-julian/36525) END SUBROUTINE caldat diff -C1 -r 206b/recipes_f-90/recipes/dfpmin.f90 208/recipes_f-90/recipes/dfpmin.f90 *** 206b/recipes_f-90/recipes/dfpmin.f90 Mon May 06 13:41:54 1996 --- 208/recipes_f-90/recipes/dfpmin.f90 Fri May 30 16:02:39 1997 *************** *** 52,54 **** sumxi=dot_product(xi,xi) ! if (fac**2 > EPS*sumdg*sumxi) then fac=1.0_sp/fac --- 52,54 ---- sumxi=dot_product(xi,xi) ! if (fac > sqrt(EPS*sumdg*sumxi)) then fac=1.0_sp/fac diff -C1 -r 206b/recipes_f-90/recipes/fit.f90 208/recipes_f-90/recipes/fit.f90 *** 206b/recipes_f-90/recipes/fit.f90 Mon May 06 13:41:55 1996 --- 208/recipes_f-90/recipes/fit.f90 Fri May 30 16:02:39 1997 *************** *** 7,9 **** REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig ! INTEGER(I4B) :: ndum REAL(SP) :: sigdat,ss,sx,sxoss,sy,st2 --- 7,9 ---- REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig ! INTEGER(I4B) :: ndata REAL(SP) :: sigdat,ss,sx,sxoss,sy,st2 *************** *** 12,14 **** if (present(sig)) then ! ndum=assert_eq(size(x),size(y),size(sig),'fit') wt=>t --- 12,14 ---- if (present(sig)) then ! ndata=assert_eq(size(x),size(y),size(sig),'fit') wt=>t *************** *** 19,21 **** else ! ndum=assert_eq(size(x),size(y),'fit') ss=real(size(x),sp) --- 19,21 ---- else ! ndata=assert_eq(size(x),size(y),'fit') ss=real(size(x),sp) *************** *** 38,39 **** --- 38,40 ---- t(:)=y(:)-a-b*x(:) + q=1.0 if (present(sig)) then *************** *** 41,43 **** chi2=dot_product(t,t) ! q=gammq(0.5_sp*(size(x)-2),0.5_sp*chi2) else --- 42,44 ---- chi2=dot_product(t,t) ! if (ndata > 2) q=gammq(0.5_sp*(size(x)-2),0.5_sp*chi2) else *************** *** 44,46 **** chi2=dot_product(t,t) - q=1.0 sigdat=sqrt(chi2/(size(x)-2)) --- 45,46 ---- diff -C1 -r 206b/recipes_f-90/recipes/frprmn.f90 208/recipes_f-90/recipes/frprmn.f90 *** 206b/recipes_f-90/recipes/frprmn.f90 Mon May 06 13:41:55 1996 --- 208/recipes_f-90/recipes/frprmn.f90 Fri May 30 16:02:39 1997 *************** *** 37,39 **** if (2.0_sp*abs(fret-fp) <= ftol*(abs(fret)+abs(fp)+EPS)) RETURN ! fp=func(p) xi=dfunc(p) --- 37,39 ---- if (2.0_sp*abs(fret-fp) <= ftol*(abs(fret)+abs(fp)+EPS)) RETURN ! fp=fret xi=dfunc(p) diff -C1 -r 206b/recipes_f-90/recipes/lnsrch.f90 208/recipes_f-90/recipes/lnsrch.f90 *** 206b/recipes_f-90/recipes/lnsrch.f90 Mon May 06 13:41:56 1996 --- 208/recipes_f-90/recipes/lnsrch.f90 Fri May 30 16:02:39 1997 *************** *** 19,22 **** INTEGER(I4B) :: ndum ! REAL(SP) :: a,alam,alam2,alamin,b,disc,f2,fold2,pabs,rhs1,rhs2,slope,& ! tmplam ndum=assert_eq(size(g),size(p),size(x),size(xold),'lnsrch') --- 19,21 ---- INTEGER(I4B) :: ndum ! REAL(SP) :: a,alam,alam2,alamin,b,disc,f2,pabs,rhs1,rhs2,slope,tmplam ndum=assert_eq(size(g),size(p),size(x),size(xold),'lnsrch') *************** *** 26,27 **** --- 25,27 ---- slope=dot_product(g,p) + if (slope >= 0.0) call nrerror('roundoff problem in lnsrch') alamin=TOLX/maxval(abs(p(:))/max(abs(xold(:)),1.0_sp)) *************** *** 42,44 **** rhs1=f-fold-alam*slope ! rhs2=f2-fold2-alam2*slope a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) --- 42,44 ---- rhs1=f-fold-alam*slope ! rhs2=f2-fold-alam2*slope a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) *************** *** 50,53 **** disc=b*b-3.0_sp*a*slope ! if (disc < 0.0) call nrerror('roundoff problem in lnsrch') ! tmplam=(-b+sqrt(disc))/(3.0_sp*a) end if --- 50,58 ---- disc=b*b-3.0_sp*a*slope ! if (disc < 0.0) then ! tmplam=0.5_sp*alam ! else if (b <= 0.0) then ! tmplam=(-b+sqrt(disc))/(3.0_sp*a) ! else ! tmplam=-slope/(b+sqrt(disc)) ! end if end if *************** *** 58,60 **** f2=f - fold2=fold alam=max(tmplam,0.1_sp*alam) --- 63,64 ---- diff -C1 -r 206b/recipes_f-90/recipes/mgfas.f90 208/recipes_f-90/recipes/mgfas.f90 *** 206b/recipes_f-90/recipes/mgfas.f90 Mon May 06 13:41:57 1996 --- 208/recipes_f-90/recipes/mgfas.f90 Fri May 30 16:02:39 1997 *************** *** 75,80 **** v=ut - tau=lop(ut)-rstrct(lop(u)) if (present(rhs)) then ! tau=rstrct(rhs)+tau else trerr=ALPHA*sqrt(sum(tau**2))/size(tau,1) --- 75,80 ---- v=ut if (present(rhs)) then ! tau=lop(ut)-rstrct(lop(u)-rhs) else + tau=lop(ut)-rstrct(lop(u)) trerr=ALPHA*sqrt(sum(tau**2))/size(tau,1) diff -C1 -r 206b/recipes_f-90/recipes/mpdiv.f90 208/recipes_f-90/recipes/mpdiv.f90 *** 206b/recipes_f-90/recipes/mpdiv.f90 Mon May 06 13:41:57 1996 --- 208/recipes_f-90/recipes/mpdiv.f90 Sat Jul 05 00:42:38 1997 *************** *** 3,5 **** USE nr, ONLY : mpinv,mpmul ! USE mpops, ONLY : mpmov,mpsub IMPLICIT NONE --- 3,5 ---- USE nr, ONLY : mpinv,mpmul ! USE mpops, ONLY : mpsad,mpmov,mpsub IMPLICIT NONE *************** *** 9,20 **** INTEGER(I4B), INTENT(IN) :: n,m ! INTEGER(I4B), PARAMETER :: MACC=3 INTEGER(I4B) :: is ! CHARACTER(1), DIMENSION(:), ALLOCATABLE :: s ! CHARACTER(1), DIMENSION(:), ALLOCATABLE, TARGET :: rr ! CHARACTER(1), DIMENSION(:), POINTER :: rr2 ! allocate(rr(2*(n+MACC)),s(n+MACC)) rr2=>rr(2:) ! call mpinv(s,v,n-m+MACC,m) ! call mpmul(rr,s,u,n-m+MACC,n) ! call mpmov(q,rr2,n-m+1) call mpmul(rr,q,v,n-m+1,m) --- 9,21 ---- INTEGER(I4B), INTENT(IN) :: n,m ! INTEGER(I4B), PARAMETER :: MACC=6 INTEGER(I4B) :: is ! CHARACTER(1), DIMENSION(:), ALLOCATABLE, TARGET :: rr,s ! CHARACTER(1), DIMENSION(:), POINTER :: rr2,s3 ! allocate(rr(2*(n+MACC)),s(2*(n+MACC))) rr2=>rr(2:) ! s3=>s(3:) ! call mpinv(s,v,n+MACC,m) ! call mpmul(rr,s,u,n+MACC,n) ! call mpsad(s,rr,n+n+MACC/2,1) ! call mpmov(q,s3,n-m+1) call mpmul(rr,q,v,n-m+1,m) diff -C1 -r 206b/recipes_f-90/recipes/mpinv.f90 208/recipes_f-90/recipes/mpinv.f90 *** 206b/recipes_f-90/recipes/mpinv.f90 Mon May 06 13:41:57 1996 --- 208/recipes_f-90/recipes/mpinv.f90 Fri Aug 15 04:19:01 1997 *************** *** 13,15 **** CHARACTER(1), DIMENSION(:), ALLOCATABLE :: rr,s ! allocate(rr(n+m+1),s(n)) mm=min(MF,m) --- 13,15 ---- CHARACTER(1), DIMENSION(:), ALLOCATABLE :: rr,s ! allocate(rr(max(n,m)+n+1),s(n)) mm=min(MF,m) diff -C1 -r 206b/recipes_f-90/recipes/mpsqrt.f90 208/recipes_f-90/recipes/mpsqrt.f90 *** 206b/recipes_f-90/recipes/mpsqrt.f90 Mon May 06 13:41:57 1996 --- 208/recipes_f-90/recipes/mpsqrt.f90 Fri May 30 16:02:40 1997 *************** *** 24,26 **** call mplsh(r,n) ! call mpmul(s,r,v,n,m) call mplsh(s,n) --- 24,26 ---- call mplsh(r,n) ! call mpmul(s,r,v,n,min(m,n)) call mplsh(s,n) *************** *** 34,36 **** end if ! call mpmul(r,u,v,n,m) call mpmov(w,r(2:),n) --- 34,36 ---- end if ! call mpmul(r,u,v,n,min(m,n)) call mpmov(w,r(2:),n) diff -C1 -r 206b/recipes_f-90/recipes/mrqmin.f90 208/recipes_f-90/recipes/mrqmin.f90 *** 206b/recipes_f-90/recipes/mrqmin.f90 Mon May 06 13:41:57 1996 --- 208/recipes_f-90/recipes/mrqmin.f90 Fri May 30 16:02:40 1997 *************** *** 44,45 **** --- 44,46 ---- call covsrt(covar,maska) + call covsrt(alpha,maska) deallocate(atry,beta,da) diff -C1 -r 206b/recipes_f-90/recipes/powell.f90 208/recipes_f-90/recipes/powell.f90 *** 206b/recipes_f-90/recipes/powell.f90 Mon May 06 13:41:58 1996 --- 208/recipes_f-90/recipes/powell.f90 Fri May 30 16:02:40 1997 *************** *** 18,19 **** --- 18,20 ---- INTEGER(I4B), PARAMETER :: ITMAX=200 + REAL(SP), PARAMETER :: TINY=1.0e-25_sp INTEGER(I4B) :: i,ibig,n *************** *** 34,37 **** call linmin(p,xit,fret) ! if (abs(fptt-fret) > del) then ! del=abs(fptt-fret) ibig=i --- 35,38 ---- call linmin(p,xit,fret) ! if (fptt-fret > del) then ! del=fptt-fret ibig=i *************** *** 39,41 **** end do ! if (2.0_sp*abs(fp-fret) <= ftol*(abs(fp)+abs(fret))) RETURN if (iter == ITMAX) call & --- 40,42 ---- end do ! if (2.0_sp*(fp-fret) <= ftol*(abs(fp)+abs(fret))+TINY) RETURN if (iter == ITMAX) call & diff -C1 -r 206b/recipes_f-90/recipes/qsimp.f90 208/recipes_f-90/recipes/qsimp.f90 *** 206b/recipes_f-90/recipes/qsimp.f90 Mon May 06 13:41:58 1996 --- 208/recipes_f-90/recipes/qsimp.f90 Fri May 30 16:02:40 1997 *************** *** 22,25 **** qsimp=(4.0_sp*st-ost)/3.0_sp ! if (abs(qsimp-os) < EPS*abs(os)) RETURN ! if (qsimp == 0.0 .and. os == 0.0 .and. j > 6) RETURN os=qsimp --- 22,27 ---- qsimp=(4.0_sp*st-ost)/3.0_sp ! if (j > 5) then ! if (abs(qsimp-os) < EPS*abs(os) .or. & ! (qsimp == 0.0 .and. os == 0.0)) RETURN ! end if os=qsimp diff -C1 -r 206b/recipes_f-90/recipes/qtrap.f90 208/recipes_f-90/recipes/qtrap.f90 *** 206b/recipes_f-90/recipes/qtrap.f90 Mon May 06 13:41:58 1996 --- 208/recipes_f-90/recipes/qtrap.f90 Fri May 30 16:02:40 1997 *************** *** 20,23 **** call trapzd(func,a,b,qtrap,j) ! if (abs(qtrap-olds) < EPS*abs(olds)) RETURN ! if (qtrap == 0.0 .and. olds == 0.0 .and. j > 6) RETURN olds=qtrap --- 20,25 ---- call trapzd(func,a,b,qtrap,j) ! if (j > 5) then ! if (abs(qtrap-olds) < EPS*abs(olds) .or. & ! (qtrap == 0.0 .and. olds == 0.0)) RETURN ! end if olds=qtrap diff -C1 -r 206b/recipes_f-90/recipes/rlft2.f90 208/recipes_f-90/recipes/rlft2.f90 *** 206b/recipes_f-90/recipes/rlft2.f90 Mon May 06 13:41:59 1996 --- 208/recipes_f-90/recipes/rlft2.f90 Fri May 30 16:08:37 1997 *************** *** 4,7 **** REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: data ! COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: spec ! COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: speq INTEGER(I4B), INTENT(IN) :: isign --- 4,7 ---- REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: data ! COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: spec ! COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: speq INTEGER(I4B), INTENT(IN) :: isign diff -C1 -r 206b/recipes_f-90/recipes/rlft3.f90 208/recipes_f-90/recipes/rlft3.f90 *** 206b/recipes_f-90/recipes/rlft3.f90 Mon May 06 13:41:59 1996 --- 208/recipes_f-90/recipes/rlft3.f90 Fri May 30 16:02:40 1997 *************** *** 4,7 **** REAL(SP), DIMENSION(:,:,:), INTENT(INOUT) :: data ! COMPLEX(SPC), DIMENSION(:,:,:), INTENT(OUT) :: spec ! COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: speq INTEGER(I4B), INTENT(IN) :: isign --- 4,7 ---- REAL(SP), DIMENSION(:,:,:), INTENT(INOUT) :: data ! COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: spec ! COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: speq INTEGER(I4B), INTENT(IN) :: isign diff -C1 -r 206b/recipes_f-90/recipes/rtsafe.f90 208/recipes_f-90/recipes/rtsafe.f90 *** 206b/recipes_f-90/recipes/rtsafe.f90 Mon May 06 13:41:59 1996 --- 208/recipes_f-90/recipes/rtsafe.f90 Fri May 30 16:02:40 1997 *************** *** 39,41 **** do j=1,MAXIT ! if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) >= 0.0 .or. & abs(2.0_sp*f) > abs(dxold*df) ) then --- 39,41 ---- do j=1,MAXIT ! if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0.0 .or. & abs(2.0_sp*f) > abs(dxold*df) ) then diff -C1 -r 206b/recipes_f-90/recipes/shoot.f90 208/recipes_f-90/recipes/shoot.f90 *** 206b/recipes_f-90/recipes/shoot.f90 Mon May 06 13:42:00 1996 --- 208/recipes_f-90/recipes/shoot.f90 Fri May 30 16:02:40 1997 *************** *** 4,6 **** USE nr, ONLY : odeint,rkqs ! USE sphoot_caller, ONLY: nvar,x1,x2 IMPLICIT NONE --- 4,6 ---- USE nr, ONLY : odeint,rkqs ! USE sphoot_caller, ONLY : nvar,x1,x2; USE ode_path, ONLY : xp,yp IMPLICIT NONE *************** *** 39,40 **** --- 39,41 ---- call load(x1,v,y) + if (associated(xp)) deallocate(xp,yp) call odeint(y,x1,x2,EPS,h1,hmin,derivs,rkqs) diff -C1 -r 206b/recipes_f-90/recipes/shootf.f90 208/recipes_f-90/recipes/shootf.f90 *** 206b/recipes_f-90/recipes/shootf.f90 Mon May 06 13:42:00 1996 --- 208/recipes_f-90/recipes/shootf.f90 Fri May 30 16:02:40 1997 *************** *** 4,6 **** USE nr, ONLY : odeint,rkqs ! USE sphfpt_caller, ONLY : x1,x2,xf,nn2 IMPLICIT NONE --- 4,6 ---- USE nr, ONLY : odeint,rkqs ! USE sphfpt_caller, ONLY : x1,x2,xf,nn2; USE ode_path, ONLY : xp,yp IMPLICIT NONE *************** *** 47,48 **** --- 47,49 ---- call load1(x1,v,y) + if (associated(xp)) deallocate(xp,yp) call odeint(y,x1,xf,EPS,h1,hmin,derivs,rkqs) diff -C1 -r 206b/recipes_f-90/recipes/simplx.f90 208/recipes_f-90/recipes/simplx.f90 *** 206b/recipes_f-90/recipes/simplx.f90 Mon May 06 13:42:00 1996 --- 208/recipes_f-90/recipes/simplx.f90 Fri May 30 16:02:40 1997 *************** *** 9,15 **** REAL(SP), PARAMETER :: EPS=1.0e-6_sp ! INTEGER(I4B) :: ip,k,kh,kp,m12,nl1,nl2,m,n INTEGER(I4B), DIMENSION(size(a,2)) :: l1 ! INTEGER(I4B), DIMENSION(size(a,1)) :: l2,l3 REAL(SP) :: bmax ! LOGICAL(LGT) :: init,proceed m=assert_eq(size(a,1)-2,size(iposv),'simplx: m') --- 9,15 ---- REAL(SP), PARAMETER :: EPS=1.0e-6_sp ! INTEGER(I4B) :: ip,k,kh,kp,nl1,m,n INTEGER(I4B), DIMENSION(size(a,2)) :: l1 ! INTEGER(I4B), DIMENSION(m2) :: l3 REAL(SP) :: bmax ! LOGICAL(LGT) :: init m=assert_eq(size(a,1)-2,size(iposv),'simplx: m') *************** *** 21,28 **** izrov(:)=l1(1:n) ! nl2=m ! l2(1:m)=arth(1,1,m) ! iposv(:)=n+l2(1:m) ! l3(1:m2)=1 init=.true. - main_procedure: do phase1: do --- 21,24 ---- izrov(:)=l1(1:n) ! iposv(:)=n+arth(1,1,m) init=.true. phase1: do *************** *** 29,37 **** if (init) then init=.false. ! if (m2+m3 == 0) then ! proceed=.true. ! exit phase1 ! else ! proceed=.false. ! end if a(m+2,1:n+1)=-sum(a(m1+2:m+1,1:n+1),dim=1) --- 25,29 ---- if (init) then + if (m2+m3 == 0) exit phase1 init=.false. ! l3(1:m2)=1 a(m+2,1:n+1)=-sum(a(m1+2:m+1,1:n+1),dim=1) *************** *** 38,41 **** end if ! kp=l1(imaxloc(a(m+2,l1(1:nl1)+1))) ! bmax=a(m+2,kp+1) phase1a: do --- 30,37 ---- end if ! if (nl1 > 0) then ! kp=l1(imaxloc(a(m+2,l1(1:nl1)+1))) ! bmax=a(m+2,kp+1) ! else ! bmax=0.0 ! end if phase1a: do *************** *** 45,52 **** else if (bmax <= EPS .and. a(m+2,1) <= EPS) then ! m12=m1+m2+1 ! do ip=m12,m if (iposv(ip) == ip+n) then ! kp=l1(imaxloc(abs(a(ip+1,l1(1:nl1)+1)))) ! bmax=a(ip+1,kp+1) ! if (bmax > 0.0) exit phase1a end if --- 41,51 ---- else if (bmax <= EPS .and. a(m+2,1) <= EPS) then ! do ip=m1+m2+1,m if (iposv(ip) == ip+n) then ! if (nl1 > 0) then ! kp=l1(imaxloc(abs(a(ip+1,l1(1:nl1)+1)))) ! bmax=a(ip+1,kp+1) ! else ! bmax=0.0 ! end if ! if (bmax > EPS) exit phase1a end if *************** *** 53,58 **** end do ! proceed=.true. ! m12=m12-1 ! where (spread(l3(1:m12-m1),2,n+1) == 1) & ! a(m1+2:m12+1,1:n+1)=-a(m1+2:m12+1,1:n+1) exit phase1 --- 52,55 ---- end do ! where (spread(l3(1:m2),2,n+1) == 1) & ! a(m1+2:m1+m2+1,1:n+1)=-a(m1+2:m1+m2+1,1:n+1) exit phase1 *************** *** 66,85 **** end do phase1a ! phase1b: do ! call simp2(m+1,n) ! if (iposv(ip) >= n+m1+m2+1) then ! k=ifirstloc(l1(1:nl1) == kp) ! nl1=nl1-1 ! l1(k:nl1)=l1(k+1:nl1+1) ! else ! if (iposv(ip) < n+m1+1) exit phase1b ! kh=iposv(ip)-m1-n ! if (l3(kh) == 0) exit phase1b ! l3(kh)=0 end if ! a(m+2,kp+1)=a(m+2,kp+1)+1.0_sp ! a(1:m+2,kp+1)=-a(1:m+2,kp+1) ! exit phase1b ! end do phase1b call swap(izrov(kp),iposv(ip)) - if (proceed) exit phase1 end do phase1 --- 63,80 ---- end do phase1a ! call simp2(m+1,n) ! if (iposv(ip) >= n+m1+m2+1) then ! k=ifirstloc(l1(1:nl1) == kp) ! nl1=nl1-1 ! l1(k:nl1)=l1(k+1:nl1+1) ! else ! kh=iposv(ip)-m1-n ! if (kh >= 1) then ! if (l3(kh) /= 0) then ! l3(kh)=0 ! a(m+2,kp+1)=a(m+2,kp+1)+1.0_sp ! a(1:m+2,kp+1)=-a(1:m+2,kp+1) ! end if end if ! end if call swap(izrov(kp),iposv(ip)) end do phase1 *************** *** 86,90 **** phase2: do ! kp=l1(imaxloc(a(1,l1(1:nl1)+1))) ! bmax=a(1,kp+1) ! if (bmax <= 0.0) then icase=0 --- 81,89 ---- phase2: do ! if (nl1 > 0) then ! kp=l1(imaxloc(a(1,l1(1:nl1)+1))) ! bmax=a(1,kp+1) ! else ! bmax=0.0 ! end if ! if (bmax <= EPS) then icase=0 *************** *** 99,104 **** call swap(izrov(kp),iposv(ip)) - if (.not. proceed) cycle main_procedure end do phase2 - exit main_procedure - end do main_procedure CONTAINS --- 98,100 ---- *************** *** 107,109 **** IMPLICIT NONE ! INTEGER(I4B) :: i,ii,k REAL(SP) :: q,q0,q1,qp --- 103,105 ---- IMPLICIT NONE ! INTEGER(I4B) :: i,k REAL(SP) :: q,q0,q1,qp *************** *** 110,121 **** ip=0 ! i=ifirstloc(a(l2(1:nl2)+1,kp+1) < -EPS) ! if (i > nl2) RETURN ! q1=-a(l2(i)+1,1)/a(l2(i)+1,kp+1) ! ip=l2(i) ! do i=i+1,nl2 ! ii=l2(i) ! if (a(ii+1,kp+1) < -EPS) then ! q=-a(ii+1,1)/a(ii+1,kp+1) if (q < q1) then ! ip=ii q1=q --- 106,116 ---- ip=0 ! i=ifirstloc(a(2:m+1,kp+1) < -EPS) ! if (i > m) RETURN ! q1=-a(i+1,1)/a(i+1,kp+1) ! ip=i ! do i=ip+1,m ! if (a(i+1,kp+1) < -EPS) then ! q=-a(i+1,1)/a(i+1,kp+1) if (q < q1) then ! ip=i q1=q *************** *** 124,126 **** qp=-a(ip+1,k+1)/a(ip+1,kp+1) ! q0=-a(ii+1,k+1)/a(ii+1,kp+1) if (q0 /= qp) exit --- 119,121 ---- qp=-a(ip+1,k+1)/a(ip+1,kp+1) ! q0=-a(i+1,k+1)/a(i+1,kp+1) if (q0 /= qp) exit *************** *** 127,129 **** end do ! if (q0 < qp) ip=ii end if --- 122,124 ---- end do ! if (q0 < qp) ip=i end if