https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109865
Bug ID: 109865 Summary: different results when routine moved inside the contains statement Product: gcc Version: og12 (devel/omp/gcc-12) Status: UNCONFIRMED Severity: normal Priority: P3 Component: fortran Assignee: unassigned at gcc dot gnu.org Reporter: Gary.White at ColoState dot edu Target Milestone: --- Created attachment 55087 --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=55087&action=edit set of subroutines where moving mc11ad inside the contains statement produces incorrect results In the following code, when the subroutine mc11ad is moved inside the contains statement, incorrect results are produced. Works fine outside the contains statement as provided here. This subroutine is only called from the main subroutine va09ad, nowhere else. Other clues are that if I put this routine inside a module, incorrect results are produced. This set of routines is used in a much larger code, so I'm not able to isolate the problem down to a simple example. I have verified that this is a gfortran bug because the code produces correct results with the Intel compiler with mc11ad inside or outside the contains statement. gfortran compiler I'm using: Using built-in specs. COLLECT_GCC=gfortran COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/12/lto-wrapper OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa OFFLOAD_TARGET_DEFAULT=1 Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion='Ubuntu 12.1.0-2ubuntu1~22.04' --with-bugurl=file:///usr/share/doc/gcc-12/README.Bugs --enable-languages=c,ada,c++,go,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-12 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --enable-cet --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-12-sZcx2y/gcc-12-12.1.0/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-12-sZcx2y/gcc-12-12.1.0/debian/tmp-gcn/usr --enable-offload-defaulted --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 12.1.0 (Ubuntu 12.1.0-2ubuntu1~22.04) Options being used to compile the code: COPTIONS = -cpp -std=f2018 -c -D ieee -D dbleprecision -m64 -fsignaling-nans -ffpe-summary='invalid','zero','overflow','underflow' -O3 -funroll-loops -ffast-math subroutine va09ad(functn,n,x,f,g,h,w,dfn,eps2,mode,maxfn,iprint,iexit,itn,parm,beta,design) use status_module, only : int32, wp, sysout, syscrn, nlogit, ncovs_nlbetas, ncovs, & zero, one, two, modnam, num_to_str, outtext implicit none integer(kind=int32), intent(in) :: n, mode, maxfn, iprint integer(kind=int32), intent(out) :: iexit, itn real(kind=wp), intent(in) :: dfn, eps2(n) real(kind=wp), intent(out) :: x(n),g(n),h(n*(n+1)/2),w(n*3), f real(kind=wp), intent(out) :: parm(nlogit), beta(ncovs_nlbetas), design(nlogit,ncovs) external functn integer(kind=int32) ig, igg, is, ir, mk, ij, i, ifn, icon real(kind=wp) alphalocal, z, gs, gys, df, gs0, tot, fy, zz, dgs, sigma, epsln character(len=:), allocatable :: frmt interface subroutine functn(nparm, xbeta, xloglk, g, parm, beta, design) use status_module use data_module implicit none integer(kind=int32), intent(in) :: nparm real(kind=wp), intent(in out) :: xbeta(nparm), g(nparm), parm(nlogit), beta(ncovs_nlbetas), design(nlogit,ncovs) real(kind=wp), intent(out) :: xloglk end subroutine functn end interface ! This interface causes problems with prcisub.f90, ! just as does including mc11ad inside the contains statement. !interface ! subroutine mc11ad(a,n,z,sigma,w,ir,mk,eps2) ! use status_module, only : wp, int32, zero, one ! implicit none ! integer(kind=int32), intent(in) :: n, mk ! integer(kind=int32), intent(out) :: ir ! real(kind=wp), intent(in out) :: a(n*(n+1)/2),z(n),w(n*3) ! real(kind=wp), intent(in) :: eps2 ! real(kind=wp), intent(in out) :: sigma ! end subroutine mc11ad !end interface if(iprint/=0) then write(unit=sysout,fmt='('' ENTERING MINIMIZATION''/)') end if ig=n igg=n+n is=igg iexit=0 ir=n if(mode==1) then ! Initialize h matrix h=zero ij=(n*(n+1)/2)+1 do i=1,n ij=ij-i h(ij)=one end do else if(mode==2) then call mc11bd(h,n,ir) if(ir<n)return end if f=zero z=f itn=0 call functn(n,x,f,g,parm,beta,design) ifn=1 df=dfn if(dfn==zero)df=f-z if(dfn<zero) df=abs(df*f) if(df<=zero)df=one 20 continue if(iprint==0) goto 21 if(mod(itn,iprint)/=0)goto 21 write(unit=sysout,fmt='(24I5)')itn,ifn,maxfn,iexit write(unit=syscrn,fmt='(24I5)')itn,ifn,maxfn,iexit write(unit=sysout,fmt='(1X,''Current Function Value = '',G15.7)') f write(unit=syscrn,fmt='(1X,''Current Function Value = '',G15.7)') f if(iprint<0) goto 21 write(unit=sysout,fmt='(1X,''Current Parameter Values:'')') write(unit=sysout,fmt='(5G15.7)') (x(i),i=1,n) write(unit=syscrn,fmt='(1X,''Current Parameter Values:'')') write(unit=syscrn,fmt='(5G15.7)') (x(i),i=1,n) write(unit=sysout,fmt='(1X,''Current Gradient:'')') write(unit=sysout,fmt='(5G15.7)') (g(i),i=1,n) write(unit=syscrn,fmt='(1X,''Current Gradient:'')') write(unit=syscrn,fmt='(5G15.7)') (g(i),i=1,n) flush(unit=sysout) 21 frmt=modnam(1:min(len(modnam),64))//' Iteration '//num_to_str(itn) call outtext(frmt) itn=itn+1 w(ig+1:ig+n)=g(1:n) call mc11ed(h,n,g,w,ir) w(is+1:is+n)=-g(1:n) gs=sum(-g(1:n)*w(ig+1:ig+n)) iexit=2 if(gs >= zero) goto 92 gs0=gs alphalocal=-two*df/gs if(alphalocal>one)alphalocal=one df=f tot=zero 30 continue iexit=3 if(ifn == maxfn) goto 92 icon=0 iexit=1 do i=1,n z=alphalocal*w(is+i) if(abs(z) >= eps2(i)) icon=1 x(i)=x(i)+z end do call functn(n,x,fy,g,parm,beta,design) ifn=ifn+1 gys=sum(g(1:n)*w(is+1:is+n)) if(fy >= f) goto 40 if(abs(gys/gs0) <= 0.9_wp) goto 50 if(gys > zero) goto 40 tot=tot+alphalocal z=10._wp if(gs < gys) z=gys/(gs-gys) if(z > 10._wp) z=10._wp alphalocal=alphalocal*z f=fy gs=gys goto 30 40 continue x(1:n)=x(1:n)-alphalocal*w(is+1:is+n) if(icon == 0) goto 92 z=3._wp*(f-fy)/alphalocal+gys+gs zz=sqrt(max(zero,z**2-gs*gys)) z=one-(gys+zz-z)/(two*zz+gys-gs) alphalocal=alphalocal*z goto 30 50 continue alphalocal=tot+alphalocal f=fy if(icon == 0) goto 90 df=df-f dgs=gys-gs0 w(igg+1:igg+n)=g(1:n) g(1:n)=-w(ig+1:ig+n) if(dgs+alphalocal*gs0<=zero) then ! Complementary dfp formula sigma=one/gs0 ir=-ir mk=1 epsln=zero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) g(1:n)=w(igg+1:igg+n)-w(ig+1:ig+n) sigma=one/(alphalocal*dgs) ir=-ir mk=0 epsln=zero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) else ! dfp formula zz=alphalocal/(dgs-alphalocal*gs0) sigma=-zz mk=1 epsln=epsilon(sigma) ! Was 1.d-16 call mc11ad(h,n,g,sigma,w,ir,mk,epsln) z=dgs*zz-one g(1:n)=w(igg+1:igg+n)+z*w(ig+1:ig+n) sigma=one/(zz*dgs**2) mk=0 epsln=zero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) end if g(1:n)=w(igg+1:igg+n) goto 20 92 continue g(1:n)=w(ig+1:ig+n) 90 continue if(iprint/=0) then write(unit=sysout,fmt='(24I5)')itn,ifn,iexit write(unit=sysout,fmt='(1X,''Current Function Value = '',G15.7)') f write(unit=sysout,fmt='(1X,''Current Parameter Values:'')') write(unit=sysout,fmt='(5G15.7)') x(1:n) write(unit=sysout,fmt='(1X,''Current Gradient:'')') write(unit=sysout,fmt='(5G15.7)') g(1:n) flush(unit=sysout) end if contains !----------------------------------------------------------------------- ! Multiply a vector by the inverse of the factors given in a subroutine mc11ed(a,n,z,w,ir) use status_module, only : wp, int32 implicit none integer(kind=int32), intent(in) :: n integer(kind=int32), intent(in) :: ir real(kind=wp) a(n*(n+1)/2),z(n),w(n*3) real(kind=wp) v integer(kind=int32) j, ij, ii, nip, i if(ir < n) return w(1)=z(1) if(n<= 1) then z(1)=z(1)/a(1) return end if do i=2,n ij=i v=z(i) do j=1,i-1 v=v-a(ij)*z(j) ij=ij+n-j end do w(i)=v z(i)=v end do z(n)=z(n)/a(ij) do nip=2,n i=n+1-nip ii=ij-nip v=z(i)/a(ii) ij=ii do j=i+1,n ii=ii+1 v=v-a(ii)*z(j) end do z(i)=v end do end subroutine mc11ed !----------------------------------------------------------------------- ! Factorize a matrix given in a subroutine mc11bd(a,n,ir) use status_module, only : wp, int32, zero implicit none integer(kind=int32), intent(in) :: n integer(kind=int32), intent(out) :: ir real(kind=wp), intent(out) :: a(n*(n+1)/2) real(kind=wp) aa, v integer(kind=int32) ii, i, jk, ij, ni, ip, ik ir=n if(n <= 1) then if(a(1) > zero)return a(1)=zero ir=0 else ii=1 do i=2,n aa=a(ii) ni=n+1+ii-i if(aa > zero)then ip=ii+1 ii=ni+1 jk=ii do ij=ip,ni v=a(ij)/aa do ik=ij,ni a(jk)=a(jk)-a(ik)*v jk=jk+1 end do a(ij)=v end do else a(ii)=zero ir=ir-1 ii=ni+1 end if end do if(a(ii) <= zero) then a(ii)=zero ir=ir-1 end if end if end subroutine mc11bd end subroutine va09ad !----------------------------------------------------------------------- subroutine mc11ad(a,n,z,sigma,w,ir,mk,eps2) use status_module, only : wp, int32, zero, one implicit none integer(kind=int32), intent(in) :: n, mk integer(kind=int32), intent(out) :: ir real(kind=wp), intent(in out) :: a(n*(n+1)/2),z(n),w(n*3) real(kind=wp), intent(in) :: eps2 real(kind=wp), intent(in out) :: sigma real(kind=wp) ti, v, tim, al, r, b, gm, y integer(kind=int32) np, ij, i, j, ip, mm ! Update factors given in a by sigma*z*zt ti=zero; v=zero; tim=zero; al=zero; r=zero; b=zero; gm=zero; y=zero; np=0; ij=0; i=0; j=0; ip=0; mm=0; if(n <= 1) then a(1)=a(1)+sigma*z(1)**2 ir=1 if(a(1) > zero) return a(1)=zero ir=0 return end if np=n+1 if(sigma > zero)goto 40 if(sigma == zero .or. ir == 0)return ti=one/sigma ij=1 if(mk /= 0) then do i=1,n if(a(ij) /= zero) ti=ti+w(i)**2/a(ij) ij=ij+np-i end do else w(1:n)=z(1:n) do i=1,n ip=i+1 v=w(i) if(a(ij) <= zero) then w(i)=zero ij=ij+np-i cycle end if ti=ti+v**2/a(ij) if(i /= n) then do j=ip,n ij=ij+1 w(j)=w(j)+v*a(ij) end do end if ij=ij+1 end do end if if(ir <= 0) goto 21 if(ti > zero) goto 22 if(mk-1 <= 0) then goto 40 else goto 23 end if 21 ti=zero ir=-ir-1 goto 23 22 ti=eps2/sigma if(eps2 == zero) ir=ir-1 23 continue mm=1 tim=ti do i=1,n j=np-i ij=ij-i if(a(ij) /= zero) tim=ti-w(j)**2/a(ij) w(j)=ti ti=tim end do goto 41 40 continue mm=0 tim=one/sigma 41 continue ij=1 do i=1,n ip=i+1 v=z(i) if(a(ij) <= zero) then if(ir > 0 .or. sigma < zero .or. v == zero) then ti=tim ij=ij+np-i cycle else ir=1-ir a(ij)=v**2/tim if(i == n) return do j=ip,n ij=ij+1 a(ij)=z(j)/v end do return end if ti=tim ij=ij+np-i cycle end if al=v/a(ij) if(mm <= 0) then ti=tim+v*al else ti=w(i) end if r=ti/tim a(ij)=a(ij)*r if(r == zero) exit if(i == n) exit b=al/ti if(r > 4._wp) then gm=tim/ti do j=ip,n ij=ij+1 y=a(ij) a(ij)=b*z(j)+y*gm z(j)=z(j)-v*y end do else do j=ip,n ij=ij+1 z(j)=z(j)-v*a(ij) a(ij)=a(ij)+b*z(j) end do end if tim=ti ij=ij+1 end do if(ir < 0) ir=-ir end subroutine mc11ad