From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 48) id A95213857720; Mon, 15 May 2023 18:19:15 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org A95213857720 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1684174755; bh=CIhVRDRvduHIquj/L95YGlW73jCpQ+/gpDXQRCycEyQ=; h=From:To:Subject:Date:From; b=p1bL4RjAUQfTb5dtbHBGaw3SjSTEqIKlf8VSUGZLmCgtgypWV7a7nHqAv/dC2cXlj q7jUTOmbf1uQmpARME/5aPPIbucqtmdDyzHBBzvYoVhiqGfXW70OVQhhXX39ncaipy zqfyp3++52K1SCuZHGuAl4gNayWFx+i3gi/zAf/k= From: "Gary.White at ColoState dot edu" To: gcc-bugs@gcc.gnu.org Subject: [Bug fortran/109865] New: different results when routine moved inside the contains statement Date: Mon, 15 May 2023 18:19:14 +0000 X-Bugzilla-Reason: CC X-Bugzilla-Type: new X-Bugzilla-Watch-Reason: None X-Bugzilla-Product: gcc X-Bugzilla-Component: fortran X-Bugzilla-Version: og12 (devel/omp/gcc-12) X-Bugzilla-Keywords: X-Bugzilla-Severity: normal X-Bugzilla-Who: Gary.White at ColoState dot edu X-Bugzilla-Status: UNCONFIRMED X-Bugzilla-Resolution: X-Bugzilla-Priority: P3 X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org X-Bugzilla-Target-Milestone: --- X-Bugzilla-Flags: X-Bugzilla-Changed-Fields: bug_id short_desc product version bug_status bug_severity priority component assigned_to reporter target_milestone attachments.created Message-ID: Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/ Auto-Submitted: auto-generated MIME-Version: 1.0 List-Id: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D109865 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=3D55087&action=3Dedit set of subroutines where moving mc11ad inside the contains statement produc= es incorrect results In the following code, when the subroutine mc11ad is moved inside the conta= ins 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 u= sed in a much larger code, so I'm not able to isolate the problem down to a sim= ple example.=20=20 I have verified that this is a gfortran bug because the code produces corre= ct results with the Intel compiler with mc11ad inside or outside the contains statement. gfortran compiler I'm using: Using built-in specs. COLLECT_GCC=3Dgfortran COLLECT_LTO_WRAPPER=3D/usr/lib/gcc/x86_64-linux-gnu/12/lto-wrapper OFFLOAD_TARGET_NAMES=3Dnvptx-none:amdgcn-amdhsa OFFLOAD_TARGET_DEFAULT=3D1 Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion=3D'Ubuntu 12.1.0-2ubuntu1~22.04' --with-bugurl=3Dfile:///usr/share/doc/gcc-12/README.= Bugs --enable-languages=3Dc,ada,c++,go,d,fortran,objc,obj-c++,m2 --prefix=3D/usr --with-gcc-major-version-only --program-suffix=3D-12 --program-prefix=3Dx86_64-linux-gnu- --enable-shared --enable-linker-build-= id --libexecdir=3D/usr/lib --without-included-gettext --enable-threads=3Dposix --libdir=3D/usr/lib --enable-nls --enable-clocale=3Dgnu --enable-libstdcxx-= debug --enable-libstdcxx-time=3Dyes --with-default-libstdcxx-abi=3Dnew --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=3Drelea= se --with-target-system-zlib=3Dauto --enable-objc-gc=3Dauto --enable-multiarch --disable-werror --enable-cet --with-arch-32=3Di686 --with-abi=3Dm64 --with-multilib-list=3Dm32,m64,mx32 --enable-multilib --with-tune=3Dgeneric --enable-offload-targets=3Dnvptx-none=3D/build/gcc-12-sZcx2y/gcc-12-12.1.0/= debian/tmp-nvptx/usr,amdgcn-amdhsa=3D/build/gcc-12-sZcx2y/gcc-12-12.1.0/deb= ian/tmp-gcn/usr --enable-offload-defaulted --without-cuda-driver --enable-checking=3Drelease --build=3Dx86_64-linux-gnu --host=3Dx86_64-linux-gnu --target=3Dx86_64-linu= x-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 =3D -cpp -std=3Df2018 -c -D ieee -D dbleprecision -m64 -fsignaling-nans -ffpe-summary=3D'invalid','zero','overflow','underflow' -O3 -funroll-loops -ffast-math=20 subroutine va09ad(functn,n,x,f,g,h,w,dfn,eps2,mode,maxfn,iprint,iexit,itn,parm,beta,de= sign) use status_module, only : int32, wp, sysout, syscrn, nlogit, ncovs_nlbetas, ncovs, & zero, one, two, modnam, num_to_str, outtext implicit none integer(kind=3Dint32), intent(in) :: n, mode, maxfn, iprint integer(kind=3Dint32), intent(out) :: iexit, itn real(kind=3Dwp), intent(in) :: dfn, eps2(n) real(kind=3Dwp), intent(out) :: x(n),g(n),h(n*(n+1)/2),w(n*3), f real(kind=3Dwp), intent(out) :: parm(nlogit), beta(ncovs_nlbetas), design(nlogit,ncovs) external functn integer(kind=3Dint32) ig, igg, is, ir, mk, ij, i, ifn, icon real(kind=3Dwp) alphalocal, z, gs, gys, df, gs0, tot, fy, zz, dg= s, sigma, epsln character(len=3D:), allocatable :: frmt interface=20 subroutine functn(nparm, xbeta, xloglk, g, parm, beta, design) use status_module=20=20 use data_module implicit none integer(kind=3Dint32), intent(in) :: nparm real(kind=3Dwp), intent(in out) :: xbeta(nparm), g(nparm), parm(nlogit), beta(ncovs_nlbetas), design(nlogit,ncovs)=20=20=20 real(kind=3Dwp), intent(out) :: xloglk end subroutine functn end interface ! This interface causes problems with prcisub.f90,=20 ! 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=3Dint32), intent(in) :: n, mk ! integer(kind=3Dint32), intent(out) :: ir ! real(kind=3Dwp), intent(in out) :: a(n*(n+1)/2),z(n),w(n*3) ! real(kind=3Dwp), intent(in) :: eps2 ! real(kind=3Dwp), intent(in out) :: sigma ! end subroutine mc11ad !end interface if(iprint/=3D0) then write(unit=3Dsysout,fmt=3D'('' ENTERING MINIMIZATION''/)') end if ig=3Dn igg=3Dn+n is=3Digg iexit=3D0 ir=3Dn if(mode=3D=3D1) then ! Initialize h matrix h=3Dzero ij=3D(n*(n+1)/2)+1 do i=3D1,n ij=3Dij-i h(ij)=3Done end do else if(mode=3D=3D2) then call mc11bd(h,n,ir) if(ir=3D zero) goto 92 gs0=3Dgs alphalocal=3D-two*df/gs if(alphalocal>one)alphalocal=3Done df=3Df tot=3Dzero 30 continue iexit=3D3 if(ifn =3D=3D maxfn) goto 92 icon=3D0 iexit=3D1 do i=3D1,n z=3Dalphalocal*w(is+i) if(abs(z) >=3D eps2(i)) icon=3D1 x(i)=3Dx(i)+z end do call functn(n,x,fy,g,parm,beta,design) ifn=3Difn+1 gys=3Dsum(g(1:n)*w(is+1:is+n)) if(fy >=3D f) goto 40 if(abs(gys/gs0) <=3D 0.9_wp) goto 50 if(gys > zero) goto 40 tot=3Dtot+alphalocal z=3D10._wp if(gs < gys) z=3Dgys/(gs-gys) if(z > 10._wp) z=3D10._wp alphalocal=3Dalphalocal*z f=3Dfy gs=3Dgys goto 30 40 continue x(1:n)=3Dx(1:n)-alphalocal*w(is+1:is+n) if(icon =3D=3D 0) goto 92 z=3D3._wp*(f-fy)/alphalocal+gys+gs zz=3Dsqrt(max(zero,z**2-gs*gys)) z=3Done-(gys+zz-z)/(two*zz+gys-gs) alphalocal=3Dalphalocal*z goto 30 50 continue alphalocal=3Dtot+alphalocal f=3Dfy if(icon =3D=3D 0) goto 90 df=3Ddf-f dgs=3Dgys-gs0 w(igg+1:igg+n)=3Dg(1:n) g(1:n)=3D-w(ig+1:ig+n) if(dgs+alphalocal*gs0<=3Dzero) then ! Complementary dfp formula sigma=3Done/gs0 ir=3D-ir mk=3D1 epsln=3Dzero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) g(1:n)=3Dw(igg+1:igg+n)-w(ig+1:ig+n) sigma=3Done/(alphalocal*dgs) ir=3D-ir mk=3D0 epsln=3Dzero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) else ! dfp formula zz=3Dalphalocal/(dgs-alphalocal*gs0) sigma=3D-zz mk=3D1 epsln=3Depsilon(sigma) ! Was 1.d-16 call mc11ad(h,n,g,sigma,w,ir,mk,epsln) z=3Ddgs*zz-one g(1:n)=3Dw(igg+1:igg+n)+z*w(ig+1:ig+n) sigma=3Done/(zz*dgs**2) mk=3D0 epsln=3Dzero call mc11ad(h,n,g,sigma,w,ir,mk,epsln) end if g(1:n)=3Dw(igg+1:igg+n) goto 20 92 continue g(1:n)=3Dw(ig+1:ig+n) 90 continue if(iprint/=3D0) then write(unit=3Dsysout,fmt=3D'(24I5)')itn,ifn,iexit write(unit=3Dsysout,fmt=3D'(1X,''Current Function Value =3D '',= G15.7)') f write(unit=3Dsysout,fmt=3D'(1X,''Current Parameter Values:'')') write(unit=3Dsysout,fmt=3D'(5G15.7)') x(1:n) write(unit=3Dsysout,fmt=3D'(1X,''Current Gradient:'')') write(unit=3Dsysout,fmt=3D'(5G15.7)') g(1:n) flush(unit=3Dsysout) end if contains =20=20=20=20=20=20=20=20 !----------------------------------------------------------------------- ! 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=3Dint32), intent(in) :: n integer(kind=3Dint32), intent(in) :: ir real(kind=3Dwp) a(n*(n+1)/2),z(n),w(n*3) real(kind=3Dwp) v integer(kind=3Dint32) j, ij, ii, nip, i if(ir < n) return w(1)=3Dz(1) if(n<=3D 1) then z(1)=3Dz(1)/a(1) return end if do i=3D2,n ij=3Di v=3Dz(i) do j=3D1,i-1 v=3Dv-a(ij)*z(j) ij=3Dij+n-j end do w(i)=3Dv z(i)=3Dv end do z(n)=3Dz(n)/a(ij) do nip=3D2,n i=3Dn+1-nip ii=3Dij-nip v=3Dz(i)/a(ii) ij=3Dii do j=3Di+1,n ii=3Dii+1 v=3Dv-a(ii)*z(j) end do z(i)=3Dv end do end subroutine mc11ed =20=20=20=20=20=20=20=20 !----------------------------------------------------------------------- ! Factorize a matrix given in a subroutine mc11bd(a,n,ir) use status_module, only : wp, int32, zero implicit none integer(kind=3Dint32), intent(in) :: n integer(kind=3Dint32), intent(out) :: ir real(kind=3Dwp), intent(out) :: a(n*(n+1)/2) real(kind=3Dwp) aa, v integer(kind=3Dint32) ii, i, jk, ij, ni, ip, ik ir=3Dn if(n <=3D 1) then if(a(1) > zero)return a(1)=3Dzero ir=3D0 else ii=3D1 do i=3D2,n aa=3Da(ii) ni=3Dn+1+ii-i if(aa > zero)then ip=3Dii+1 ii=3Dni+1 jk=3Dii do ij=3Dip,ni v=3Da(ij)/aa do ik=3Dij,ni a(jk)=3Da(jk)-a(ik)*v jk=3Djk+1 end do a(ij)=3Dv end do else a(ii)=3Dzero ir=3Dir-1 ii=3Dni+1 end if end do if(a(ii) <=3D zero) then a(ii)=3Dzero ir=3Dir-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=3Dint32), intent(in) :: n, mk integer(kind=3Dint32), intent(out) :: ir real(kind=3Dwp), intent(in out) :: a(n*(n+1)/2),z(n),w(n*3) real(kind=3Dwp), intent(in) :: eps2 real(kind=3Dwp), intent(in out) :: sigma real(kind=3Dwp) ti, v, tim, al, r, b, gm, y integer(kind=3Dint32) np, ij, i, j, ip, mm ! Update factors given in a by sigma*z*zt ti=3Dzero; v=3Dzero; tim=3Dzero; al=3Dzero; r=3Dzero; b=3Dzero; gm= =3Dzero; y=3Dzero; np=3D0; ij=3D0; i=3D0; j=3D0; ip=3D0; mm=3D0; if(n <=3D 1) then a(1)=3Da(1)+sigma*z(1)**2 ir=3D1 if(a(1) > zero) return a(1)=3Dzero ir=3D0 return end if np=3Dn+1 if(sigma > zero)goto 40 if(sigma =3D=3D zero .or. ir =3D=3D 0)return ti=3Done/sigma ij=3D1 if(mk /=3D 0) then do i=3D1,n if(a(ij) /=3D zero) ti=3Dti+w(i)**2/a(ij) ij=3Dij+np-i end do else w(1:n)=3Dz(1:n) do i=3D1,n ip=3Di+1 v=3Dw(i) if(a(ij) <=3D zero) then w(i)=3Dzero ij=3Dij+np-i cycle end if ti=3Dti+v**2/a(ij) if(i /=3D n) then do j=3Dip,n ij=3Dij+1 w(j)=3Dw(j)+v*a(ij) end do end if ij=3Dij+1 end do end if if(ir <=3D 0) goto 21 if(ti > zero) goto 22 if(mk-1 <=3D 0) then goto 40 else goto 23 end if 21 ti=3Dzero ir=3D-ir-1 goto 23 22 ti=3Deps2/sigma if(eps2 =3D=3D zero) ir=3Dir-1 23 continue mm=3D1 tim=3Dti do i=3D1,n j=3Dnp-i ij=3Dij-i if(a(ij) /=3D zero) tim=3Dti-w(j)**2/a(ij) w(j)=3Dti ti=3Dtim end do goto 41 40 continue mm=3D0 tim=3Done/sigma 41 continue ij=3D1 do i=3D1,n ip=3Di+1 v=3Dz(i) if(a(ij) <=3D zero) then if(ir > 0 .or. sigma < zero .or. v =3D=3D zero) then ti=3Dtim ij=3Dij+np-i cycle else ir=3D1-ir a(ij)=3Dv**2/tim if(i =3D=3D n) return do j=3Dip,n ij=3Dij+1 a(ij)=3Dz(j)/v end do return end if ti=3Dtim ij=3Dij+np-i cycle end if al=3Dv/a(ij) if(mm <=3D 0) then ti=3Dtim+v*al else ti=3Dw(i) end if r=3Dti/tim a(ij)=3Da(ij)*r if(r =3D=3D zero) exit if(i =3D=3D n) exit b=3Dal/ti if(r > 4._wp) then gm=3Dtim/ti do j=3Dip,n ij=3Dij+1 y=3Da(ij) a(ij)=3Db*z(j)+y*gm z(j)=3Dz(j)-v*y end do else do j=3Dip,n ij=3Dij+1 z(j)=3Dz(j)-v*a(ij) a(ij)=3Da(ij)+b*z(j) end do end if tim=3Dti ij=3Dij+1 end do if(ir < 0) ir=3D-ir end subroutine mc11ad=