0

Let me provide some background information for the problem I’m facing in Fortran.

• When I run the scripts on my laptop without parallelizing, no problem occurs.

When I run the scripts on my laptop using OpenMP (USE OMP_LIB), a crash occurs. • After doing some inspection, it looks like with OpenMP the code crashes when calling function zbrent. My guess is that the function funcz is not handled properly with the INTERFACE to run the codes with OpenMP. • For clarity purposes, I reduce the function zbrent to the following script below, in which the crash still happens. • Also, below it you can find the script where I call zbrent.

!*******************************************************************                    
function zbrent(funcz,x1,x2,tol )

        implicit none
        real(8), intent(in) :: x1,x2,tol
        real(8) :: zbrent
        
        INTERFACE
            function funcz(x)

                implicit none
              
                real(8), intent(in) :: x
                real(8) :: funcz
       
            END function funcz
        END INTERFACE

        integer, parameter :: ITMAX=200
        real(8), parameter :: EPS=epsilon(x1)
        integer :: iter_z
        real(8) :: a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
                        
        a=x1
        b=x2
        
        fa=funcz(a)
        fb=funcz(b)
        
        if ((fa>0d0 .and. fb>0d0) .or. (fa<0d0 .and. fb<0d0)) then
            print*,'root must be bracketed for zbrent ' 
            stop
        endif
        c=b
        fc=fb
        
        zbrent=0.5d0*(x1+x2)

end function 


!*******************************************************************  



    DOUBLE PRECISION function a_fun(b,a)

        USE ZBREN_INT

        IMPLICIT NONE

        REAL(8), INTENT(IN)::b,a
        INTEGER, PARAMETER:: MAXFN = 1000
        DOUBLE PRECISION, PARAMETER:: ERRREL = 1D-10, &
            ERRABS = 1D-10
        DOUBLE PRECISION :: left, right, a_max1,a_min1, dif_right, dif_left

        a_max1 = rhoa*a + (1-rhoa)*mua + width*epsmax
        a_min1 = rhoa*a + (1-rhoa)*mua + width*epsmin

        dif_left  = dif_fun(a_min1,b)
        dif_right = dif_fun(a_max1,b)

        if (dif_right*dif_left<0) then  

            IF (dif_left.LE.-1d3+100d0) THEN

                if (dif_left>0) then
                    a_fun = a_min1
                ELSE
                    left =  a_min1
                    right = a_max1
                    IF (SIGN(1d0,left*right).NE.-1d0) THEN
                        PRINT*,'dif_left,dif_right',dif_left,dif_right
                        PAUSE
                        STOP
                    ENDIF
                    a_fun = zbrent(dif_fun2,left,right,ERRABS)
                endif
            else

                left =  a_min1
                right = a_max1
                IF (SIGN(1d0,left*right).NE.-1d0) THEN
                    PRINT*,'dif_left,dif_right',dif_left,dif_right
                    PAUSE
                    STOP
                ENDIF
                 a_fun = zbrent(dif_fun2,left,right,ERRABS)
            endif

        else if (dif_left>0) then
            a_fun = a_min1
        else 
            a_fun = a_max1
        end if

    RETURN

    CONTAINS

        DOUBLE PRECISION FUNCTION dif_fun2(x)

            DOUBLE PRECISION, INTENT(IN)::x

            dif_fun2 = dif_fun(x,b)

        end function

    end function a_fun

DOUBLE PRECISION function dif_fun(a,b)

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::a,b
DOUBLE PRECISION::dif_fun_temp

dif_fun = v0_fun(b,a)-vaut_fun(a)

end function dif_fun
!!*********************************************************************************

As mentioned, this script has been run without OpenMP and I believe I have narrowed down the issue to the code chunk above. Any insight is appreciated, let me know if further information is necessary.

bellman_1913
  • 11
  • 3
  • The line `use omp_lib` does not enable OpenMP per se, it only imports the interfaces of the subroutines and functions defined in that module. Very often, this `use omp_lib` is not used at all. Also, the magic number 8 in `real(8)` is quite ugly, see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Vladimir F Героям слава Apr 28 '23 at 21:12
  • There is no threading in the functions that you show and there is no way for us to test the code. See [ask] and [mcve]. Show us at least the code that calls the functions including all relevant declarations. Best, a full program that we can test. – Vladimir F Героям слава Apr 28 '23 at 21:15
  • @VladimirFГероямслава apologies, am new to Fortran and Stack Overflow! i think the code is too long to post to SO despite my efforts to shorten it. I think that the problem can come from some conflict between IMSL functions and others when parallelizing. Is there any way to post the vfproj file or the three source files for it to stack overflow or does the code have to be in a text box? – bellman_1913 May 01 '23 at 19:52

0 Answers0