Cube Root Error

[Frank_M]

    integer(4)                                  :: i                ! a counter
    real(8),    dimension(pdlen)                :: tmp1             ! temp array for calculating
    real(8),    dimension(pdlen)                :: tmp2             ! temp array for calculating 
    
    <bla, bla bla >
 
    
    tmp1 = ABS(tmp2)
    
    do i = 1, pdlen
        skew(i) = tmp1(i)**(1/3)
    end do

   In the above code skew is being returned as an array of 1.0

   The values in tmp1 are small, but SQRT(tmp1) returns the correct answers.

  What am I doing wrong here?

[mecej4]

1/3 = 0, according to the rules of integer arithmetic.

[Frank_M]

 do i = 1, pdlen
        skew(i) = tmp1(i)**DBLE(1/3)
    end do

Changed it to the above, it still returns 1.0

[Frank_M]

     Lacking a CBRT() function, is there any way short of  Newton-Raphson of taking a cube root in Intel fortran?

[Frank_M]

Never mind.

I found the error.

missing a dam dot.

[mecej4]

If you replace (1/3) with (1./3) or (1/3.) or (1./3.), the expression will be evaluated in single-precision (24 bits). You can make the expression to come out in double precision by writing 'd0' after the decimal point, or by writing dble(1)/3, etc.

[Jacob_Williams · Accepted solution]

The correct and portable way is to use a function like this (will give the real solution with the correct sign).  

!cube root : real solution
pure function cube_root(x) result(y)
    
    use, intrinsic :: iso_fortran_env, only: wp => real64  !double precision
    
    implicit none
    
    real(wp),intent(in) :: x
    real(wp) :: y
    
    real(wp),parameter :: one_third = 1.0_wp / 3.0_wp
    
    y = sign( abs(x) ** one_third , x )
    
end function cube_root

[Frank_M]

Yeah thanks for the help.

That final bit of code looks like this:

    tmp1 = ABS(tmp2)
    
    do i = 1, pdlen
        skew(i) = tmp1(i)**(1/3.0_8)
        skew(i) = SIGN(skew(i),tmp2(i))
    end do

The real long term solution is to start wearing my glasses.

[jimdempseyatthecove]

>>The real long term solution is to start wearing my glasses.

Ah, unlike experience, that problem tends to get worse with age.

Jim Dempsey

[Jacob_Williams]

FYI: note that if you replace "pure" in my cube_root function with "pure elemental", all your code can be replaced with:

skew = cube_root(tmp2)

No need for the loop and temporary variable.

[TimP]

Frank_M wrote:

     Lacking a CBRT() function, is there any way short of  Newton-Raphson of taking a cube root in Intel fortran?

The examples given in this thread, when vectorized, will call the library cube root function which may be expected to use Newton iterations.  When not optimized, **(1/3.) would be expected to involve log() and exp().  The underlying library functions might be expected to be the same with Intel Fortran and C++.  If you wish, you could call the cbrt function via C interop, but this would miss opportunities for optimization.