Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
29422 Discussions

convert where to if else style code in fortran.

aketh_t_
Beginner
2,118 Views

hi below is the code I rewrote(with do loops) and further down(with where statements) is the original code. The code seems to be producing differnt values for WORK arrays over time. How to resolve this issue?? 

do j=1,ny_block
           do i=1,nx_block

           if(LMASK(i,j) .eq. .true. )then

            WORK1(i,j,kk) =  KAPPA_THIC(i,j,kbt,k,bid)  &
                           * SLX(i,j,kk,kbt,k,bid) * dz(k)

            WORK2(i,j,kk) = c2 * dzwr(k) * ( WORK1(i,j,kk)            &
              - KAPPA_THIC(i,j,ktp,k+1,bid) * SLX(i,j,kk,ktp,k+1,bid) &
                                            * dz(k+1) )

            WORK2_NEXT(i,j) = c2 * ( &
              KAPPA_THIC(i,j,ktp,k+1,bid) * SLX(i,j,kk,ktp,k+1,bid) - &
              KAPPA_THIC(i,j,kbt,k+1,bid) * SLX(i,j,kk,kbt,k+1,bid) )

            WORK3(i,j,kk) =  KAPPA_THIC(i,j,kbt,k,bid)  &
                           * SLY(i,j,kk,kbt,k,bid) * dz(k)

            WORK4(i,j,kk) = c2 * dzwr(k) * ( WORK3(i,j,kk)            &
              - KAPPA_THIC(i,j,ktp,k+1,bid) * SLY(i,j,kk,ktp,k+1,bid) &
                                            * dz(k+1) )

            WORK4_NEXT(i,j) = c2 * ( &
              KAPPA_THIC(i,j,ktp,k+1,bid) * SLY(i,j,kk,ktp,k+1,bid) - &
              KAPPA_THIC(i,j,kbt,k+1,bid) * SLY(i,j,kk,kbt,k+1,bid) )

            endif

           enddo
          enddo

         

           where ( LMASK )

            WORK1(:,:,kk) =  KAPPA_THIC(:,:,kbt,k,bid)  &
                           * SLX(:,:,kk,kbt,k,bid) * dz(k)
            WORK2(:,:,kk) = c2 * dzwr(k) * ( WORK1(:,:,kk)            &
              - KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) &
                                            * dz(k+1) )

            WORK2_NEXT = c2 * ( &
              KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) - &
              KAPPA_THIC(:,:,kbt,k+1,bid) * SLX(:,:,kk,kbt,k+1,bid) )

            WORK3(:,:,kk) =  KAPPA_THIC(:,:,kbt,k,bid)  &
                           * SLY(:,:,kk,kbt,k,bid) * dz(k)
            WORK4(:,:,kk) = c2 * dzwr(k) * ( WORK3(:,:,kk)            &
              - KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) &
                                            * dz(k+1) )

            WORK4_NEXT = c2 * ( &
              KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) - &
              KAPPA_THIC(:,:,kbt,k+1,bid) * SLY(:,:,kk,kbt,k+1,bid) )

          endwhere

Are they equivalent??

 

Coz I am having issues of correctness here.

0 Kudos
14 Replies
TimP
Honored Contributor III
2,118 Views

I would say they aren't equivalent, as the explicit DO loops forbid the compiler from switching loop nest unless it can "prove" correctness, while the multiple-rank syntax gives the compiler less guidance, and gets into insufficiently tested territory.

0 Kudos
aketh_t_
Beginner
2,118 Views

I am so sorry the second one is the original code and first is changed code.

at the end of loops assuming inputs are the same do they produce same value for work arrays???

0 Kudos
aketh_t_
Beginner
2,118 Views

At the end I find work arrays are producing different values. I dont know how to ensure correctness while rewriting my code.

0 Kudos
FortranFan
Honored Contributor III
2,118 Views

Try this:

   where ( LMASK )
      
      WORK1(:,:,kk) =  KAPPA_THIC(:,:,kbt,k,bid)  &
         * SLX(:,:,kk,kbt,k,bid) * dz(k)
         
      WORK3(:,:,kk) =  KAPPA_THIC(:,:,kbt,k,bid)  &
         * SLY(:,:,kk,kbt,k,bid) * dz(k)
         
   end where
   
   where ( LMASK )
      
      WORK2(:,:,kk) = c2 * dzwr(k) * ( WORK1(:,:,kk)            &
         - KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) &
         * dz(k+1) )
      
      WORK2_NEXT = c2 * ( &
         KAPPA_THIC(:,:,ktp,k+1,bid) * SLX(:,:,kk,ktp,k+1,bid) - &
         KAPPA_THIC(:,:,kbt,k+1,bid) * SLX(:,:,kk,kbt,k+1,bid) )
      
      WORK4(:,:,kk) = c2 * dzwr(k) * ( WORK3(:,:,kk)            &
         - KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) &
         * dz(k+1) )
      
      WORK4_NEXT = c2 * ( &
         KAPPA_THIC(:,:,ktp,k+1,bid) * SLY(:,:,kk,ktp,k+1,bid) - &
         KAPPA_THIC(:,:,kbt,k+1,bid) * SLY(:,:,kk,kbt,k+1,bid) )
      
   end where

aketh t. wrote:

At the end I find work arrays are producing different values. I dont know how to ensure correctness while rewriting my code.

You can create a minimal working example using your multi-dimensional arrays and try it out both ways, with DO loop and WHERE construct.

0 Kudos
aketh_t_
Beginner
2,118 Views

err I am sorry.

I will change the main thread contents.

The first code is the changed second is the original. I want do loops in my code for openmp support.

0 Kudos
mecej4
Honored Contributor III
2,118 Views

Apart from the other issues, you may have created a subtle bug with the seemingly benign test in the line

     if(LMASK(i,j) .eq. .true. )then

Presumably, LMASK is of type LOGICAL.

The issue is covered by Doctor Fortran in "To .EQV. or to .NEQV., that is the question", or "It's only LOGICAL", see https://software.intel.com/en-us/forums/topic/275071#comment-1548435 .

0 Kudos
aketh_t_
Beginner
2,118 Views

wont 

if(LMASK(i,j))then

.

.

.

endif

be sufficient

0 Kudos
mecej4
Honored Contributor III
2,118 Views

Replying to #8:

     Yes; sufficient, but also LOGICAL!

0 Kudos
Ron_Green
Moderator
2,118 Views

the 2 versions may have the expressions performed in different orders, and hence you can get differing results.  DO loops are optimization/vectorization candidates.

you may get closer to identical results if you use -O0

If you insist on using optimization, you can look at the optimization report via -qopt-report5 and you'll see differences in how these blocks are treated.

I can quite easily see where the results would vary slightly at optimization.

Here is a good read on this topic:  https://software.intel.com/sites/default/files/managed/f1/35/fp-control-2013-03.pdf

ron

0 Kudos
jimdempseyatthecove
Honored Contributor III
2,118 Views

aketh,

The point you are missing from mecej4's post is: We know not what your LMASK is, and equally important, how you set the values into LMASK.

In Fortran (default settings) .true. is indicated by the least significant bit in the variable being set to 1. .false. is indicated by the least significant bit of the variable being 0.

Therefore, for logicals, use .EQV. and .NEQV., do not use == as when the value held in the logical variable is numerically 1, the variable may be logically .true. but not equal to logical .true..

IF(someLogical == .true. ) xxx

May (may) require the integer equivalent value within your someLogical to be -1 in order to evaluate to .true.

Whereas

IF(someLogical .EQV. .true.) xxx

Makes the decision based only on the lsb of what is in your someLogical.

As long as you do not set a logical variable with an integer value, then == and .EQV. are (likely) interchangeable..

The link that mecj4 gave you is worth reading closely.

If by chance (or intention) LMASK contains /=.true. values that are .true. then use:

if(LMASK(i,j))then

As this will evaluate only the lsb of what is inside LMASK(i,j)

Jim Dempsey

0 Kudos
aketh_t_
Beginner
2,118 Views

Hi, 

To me variation in data is okay, but what I would like is if the data are EQUIVALENT.

1.222222 and 1.222223 are okay.

0 Kudos
aketh_t_
Beginner
2,118 Views

LMASK is of type logical

logical :: LMASK(x,y)

0 Kudos
jimdempseyatthecove
Honored Contributor III
2,118 Views

aketh t. wrote:

Hi, 

To me variation in data is okay, but what I would like is if the data are EQUIVALENT.

1.222222 and 1.222223 are okay.

Assuming the left number is in A and the right number is in B

The expression A*EPSILON(A) will yield a number that is approximately equivalent to the value of the least significant bit in the mantissa of A. You may need (want) to use ABS(A)*EPSILON(A)

DELTA =  ABS(A)*EPSILON(A)
IF(((A-DELTA).LE.B).AND.((A+DELTA).GE.B)) ...

Caution, you may need a DELTA that is larger that obtained above. It is your responsibility to determine a suitable value.

You can also use the intrinsic function NEAREST

IF((NEAREST(A,-1.0) .LE. B).AND.(NEAREST(A,1.0).GE.B)) ...

Again, convergence routines may require a larger degree of uncertainty (or certainty).

Jim Dempsey

0 Kudos
aketh_t_
Beginner
2,118 Views

hi all,

I managed to solve my issue with all your help. thanks you all.

Very often it happens we solve our problems and simply dont have the courtesy to thank you all once done.

therefore here is my special thanks to you all, without whom I would never be able go along with my efforts at parallelization and vectorization of CESM climate code.

I realize how you cater to the a huge community, which must be a thankless job.

 

0 Kudos
Reply