Spurious results when compiling with /O2 and OpenMP, version 15.0.1

[Kevin_McGrattan]

We recently upgraded our 64 bit Fortran compilers on Windows, OSX and linux to 15.0.1.148 Build 20141023. One of our nightly verification cases failed when using the Windows release compilation (/O2 /Qopenmp). The failure was substantial -- a totally spurious result. All other builds worked (linux, OSX, Windows debug). This is a very large computational fluids code broken up into about 30 source files. The code is fully compliant with Fortran 2003 and we run the thread checker on hundreds of test cases each day. The thread checker did not detect a problem, and we have run these test cases with previous versions from roughly 9 through 13. I tracked down the problem to one source file, which contained a single OpenMP construct. When I removed the OpenMP comments, but still compiled with /Qopenmp, the case still failed. When I compiled this one source file without the /Qopenmp, the case worked. I even got the case to work with OpenMP when I added a write statement to one of the loops -- a loop that was not parallelized with OpenMP. This leads me to suspect that there is a problem with the /O2 optimization (/O1 works fine) combined with /Qopenmp. It would be difficult to submit the entire code and test case, but can you think of something that might have changed in version 15 vis a vis OpenMP and optimization?

[Martyn_C_Intel]

According to the documentation, /Qtrapuv  sets /Qinit:snan  and so implicitly also /fpe:0.

If you can provide just the source for rerelhole.f, your command line and possibly one or two important declarations if they are hidden in include files or modules, (your previous source didn't need any), we can try to reproduce. That should probably be a new thread, however.

If the problem is cause by /Qtrapuv, you could try using /check:uninit instead.

[Michel]

I have updated the compiler to Update 2 (2015.2.179) and recompiled all source.

I still have this problem with /O2

Is there any view on when this issue will be solved?
Is there an internal tracking ID?

Michel.

[Steven_L_Intel1]

The issue ID is DPD200364973  It is fixed for a major release later this year.

[Michel]

Thank you Steve.

Hopefully it would not take too long, we are not really happy to release our product with degraded performance.
And we don't want to go back to the last 2014 compiler.

Michel.

[Steven_L_Intel1]

I have asked the developers if they would be willing to put the fix into an earlier update. A lot depends on how complex the fix was, as it may be that they want the extended beta test of a major release to make sure there are no problems. I will let you know what they tell me.

[Martyn_C_Intel]

Michel,

            The developers have tried to put a fix into the 15.0 compiler, but found that this had consequences for many other loops. They don't think it is safe to make such a change in a released product. They suggest that you post-increment the index ap in your loops, instead of pre-incrementing it. This avoids the issues that you have been seeing, but still allows the inner loop to vectorize. The simplest way to do this:

      subroutine prms2d3l(fnc,prms,bp,ep,num,tmat,vw,d3l,ap)

      implicit none

      include '../com/panel3d.i'

      integer*2 fnc(*),num,ap,i,j,acp
      real*4 prms(3,*),bp(2,*),ep(2,*),vw,d3l(3,*),cp(2,360)
      real*8 tmat(3,3),bp8(2),ep8(2),mp8(2),r8

      cirseg=5.0
      ap=0
      do i=1,num
       ap=ap+1
       d3l(1,ap)=bp(1,i)*tmat(1,1)+vw*tmat(2,1)+bp(2,i)*tmat(3,1)
       d3l(2,ap)=bp(1,i)*tmat(1,2)+vw*tmat(2,2)+bp(2,i)*tmat(3,2)
       d3l(3,ap)=bp(1,i)*tmat(1,3)+vw*tmat(2,3)+bp(2,i)*tmat(3,3)
       if (fnc(i).eq.1.and.abs(prms(3,i)).gt.1e-6) then
        bp8(1)=bp(1,i)
        bp8(2)=bp(2,i)
        ep8(1)=ep(1,i)
        ep8(2)=ep(2,i)
        mp8(1)=prms(1,i)
        mp8(2)=prms(2,i)
        r8=prms(3,i)
        call prms2cp4(bp8,ep8,mp8,r8,acp,cp)
        ap=ap+1
        do j=1,acp
         d3l(1,ap)=cp(1,j)*tmat(1,1)+vw*tmat(2,1)+cp(2,j)*tmat(3,1)
         d3l(2,ap)=cp(1,j)*tmat(1,2)+vw*tmat(2,2)+cp(2,j)*tmat(3,2)
         d3l(3,ap)=cp(1,j)*tmat(1,3)+vw*tmat(2,3)+cp(2,j)*tmat(3,3)
         ap=ap+1
        enddo
        ap=ap-1       
       endif
      enddo
      ap=ap+1
      d3l(1,ap)=ep(1,num)*tmat(1,1)+vw*tmat(2,1)+ep(2,num)*tmat(3,1)
      d3l(2,ap)=ep(1,num)*tmat(1,2)+vw*tmat(2,2)+ep(2,num)*tmat(3,2)
      d3l(3,ap)=ep(1,num)*tmat(1,3)+vw*tmat(2,3)+ep(2,num)*tmat(3,3)
      return
      end

Or if you choose to post-increment everywhere, including the outer loop, it might look like:

     subroutine prms2d3l(fnc,prms,bp,ep,num,tmat,vw,d3l,ap)

      implicit none

      include '../com/panel3d.i'

      integer*2 fnc(*),num,ap,i,j,acp
      real*4 prms(3,*),bp(2,*),ep(2,*),vw,d3l(3,*),cp(2,360)
      real*8 tmat(3,3),bp8(2),ep8(2),mp8(2),r8

      cirseg=5.0
      ap=1
      do i=1,num
       d3l(1,ap)=bp(1,i)*tmat(1,1)+vw*tmat(2,1)+bp(2,i)*tmat(3,1)
       d3l(2,ap)=bp(1,i)*tmat(1,2)+vw*tmat(2,2)+bp(2,i)*tmat(3,2)
       d3l(3,ap)=bp(1,i)*tmat(1,3)+vw*tmat(2,3)+bp(2,i)*tmat(3,3)
       if (fnc(i).eq.1.and.abs(prms(3,i)).gt.1e-6) then
        bp8(1)=bp(1,i)
        bp8(2)=bp(2,i)
        ep8(1)=ep(1,i)
        ep8(2)=ep(2,i)
        mp8(1)=prms(1,i)
        mp8(2)=prms(2,i)
        r8=prms(3,i)
        call prms2cp4(bp8,ep8,mp8,r8,acp,cp)
        ap=ap+1
        do j=1,acp
         d3l(1,ap)=cp(1,j)*tmat(1,1)+vw*tmat(2,1)+cp(2,j)*tmat(3,1)
         d3l(2,ap)=cp(1,j)*tmat(1,2)+vw*tmat(2,2)+cp(2,j)*tmat(3,2)
         d3l(3,ap)=cp(1,j)*tmat(1,3)+vw*tmat(2,3)+cp(2,j)*tmat(3,3)
         ap=ap+1
        enddo
        ap=ap-1       
       endif
       ap=ap+1       
      enddo
      d3l(1,ap)=ep(1,num)*tmat(1,1)+vw*tmat(2,1)+ep(2,num)*tmat(3,1)
      d3l(2,ap)=ep(1,num)*tmat(1,2)+vw*tmat(2,2)+ep(2,num)*tmat(3,2)
      d3l(3,ap)=ep(1,num)*tmat(1,3)+vw*tmat(2,3)+ep(2,num)*tmat(3,3)
      return
      end

Hopefully, this will allow you to release your product without any degradation in performance.

[Martyn_C_Intel]

An update to the version 16.0 beta compiler will contain a fix for this issue. The beta test program for the 16.0 compiler has just started. If you wish to learn more about the beta program, and perhaps enroll and test this issue when the fix becomes available, you may do so at

https://software.intel.com/articles/intel-parallel-studio-xe-2016-beta

[Michel]

I have just registered for the version 16 beta and will start checking when I have the time.

We are still only compiling with /O1

Michel.

[Martyn_C_Intel]

The fix is not in the initial beta, it is expected to be in update 1, which should be available in a couple of weeks.

I would not start testing until then.

[Michel]

I ran a full build with the beta compiler using /O3 and ran a quick test for one identified problem.

The result was positive. But I will hold the full test until after update 1 is released.

[Steven_L_Intel1]

I would not assume that you are seeing the same problem. But I do recommend you install 15.0.4 at least and see what you get.

16.0 should be out by the end of August.

[John_Campbell]

Although the error has been identified as a compiler problem, I wonder if a code solution would be to clean up some of the vectorization problems.

A change to the problem loop in post #18 could be:

real*8 vec(3)
  ..
  vec(2) = vw
  do j=1,acp 
   ap=ap+1 
!   d3l(1,ap)=cp(1,j)*tmat(1,1)+vw*tmat(2,1)+cp(2,j)*tmat(3,1) 
!   d3l(2,ap)=cp(1,j)*tmat(1,2)+vw*tmat(2,2)+cp(2,j)*tmat(3,2) 
!   d3l(3,ap)=cp(1,j)*tmat(1,3)+vw*tmat(2,3)+cp(2,j)*tmat(3,3) 
   vec(1) = cp(1,j)
   vec(3) = cp(2,j)
   d3l(1,ap) = dot_product (vec, tmat(1:3,1)) 
   d3l(2,ap) = dot_product (vec, tmat(1:3,2))
   d3l(3,ap) = dot_product (vec, tmat(1:3,3))
  enddo 

!or
  vec(2) = vw
  do j=1,acp 
   ap=ap+1 
   vec(1) = cp(1,j)
   vec(3) = cp(2,j)
   do k = 1,3
     d3l(k,ap) = dot_product (vec, tmat(1:3,k)) 
   end do
  enddo 

This would provide dot_products with no mixing of real*4 and *8. If successful, this could be applied similar loops. It could remove some of the complexities for required optimisation.

John

[mecej4]

This simplification can be taken a step further; lines 22-24 of #34 can be replaced by:

     d3l(1:3,ap) = matmul (vec, tmat(1:3,1:3))

to give d = v^TT, where d is the vector in d3l, v is the vector in vec and T is the 3 X 3 matrix.

One of the MKL/Lapack ?gemv routines can be called instead. Testing would be required to ascertain which version gives the best performance.