public inbox for gcc-bugs@sourceware.org
help / color / mirror / Atom feed
* [Bug fortran/109865] New: different results when routine moved inside the contains statement
@ 2023-05-15 18:19 Gary.White at ColoState dot edu
  2023-05-15 18:32 ` [Bug fortran/109865] " anlauf at gcc dot gnu.org
                   ` (16 more replies)
  0 siblings, 17 replies; 18+ messages in thread
From: Gary.White at ColoState dot edu @ 2023-05-15 18:19 UTC (permalink / raw)
  To: gcc-bugs

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

^ permalink raw reply	[flat|nested] 18+ messages in thread
end of thread, other threads:[~2023-05-17 16:55 UTC | newest]

Thread overview: 18+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2023-05-15 18:19 [Bug fortran/109865] New: different results when routine moved inside the contains statement Gary.White at ColoState dot edu
2023-05-15 18:32 ` [Bug fortran/109865] " anlauf at gcc dot gnu.org
2023-05-15 18:50 ` kargl at gcc dot gnu.org
2023-05-15 19:11 ` Gary.White at ColoState dot edu
2023-05-15 19:24 ` sgk at troutmask dot apl.washington.edu
2023-05-15 20:32 ` Gary.White at ColoState dot edu
2023-05-15 20:41 ` kargl at gcc dot gnu.org
2023-05-15 21:54 ` Gary.White at ColoState dot edu
2023-05-15 22:01 ` Gary.White at ColoState dot edu
2023-05-16 15:55 ` Gary.White at ColoState dot edu
2023-05-16 21:51 ` sgk at troutmask dot apl.washington.edu
2023-05-16 22:32 ` Gary.White at ColoState dot edu
2023-05-17  0:56 ` Gary.White at ColoState dot edu
2023-05-17  0:56 ` Gary.White at ColoState dot edu
2023-05-17  1:09 ` Gary.White at ColoState dot edu
2023-05-17  3:08 ` jvdelisle at gcc dot gnu.org
2023-05-17 16:49 ` Gary.White at ColoState dot edu
2023-05-17 16:55 ` pinskia at gcc dot gnu.org

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).