CHALLENGE using root_fortran to find root of a complex function

The use of dlog(), dsin(), dsqrt(), etc. in new code should be discouraged, just like using D exponents in floating point constants. With the modern fortran KIND system, there are better ways to write all those expressions.

1 Like

Re: “doesn’t make difference but I guess that in other cases it does” - which other context? Unless you are working with a processor does not conform even to Fortran 90 standard revision (from >30 years ago), the use of the generic intrinsics such as LOG is strongly recommended.

Starting with Fortran 2018 from over six years ago, the standard also says:

It is highly questionable as to whether that is “fine”.

Standard-wise, that is acceptable with defined behavior. In fact, you can even use an integer literal, 3 / acos( -1.0_wp )

But mixed-mode arithmetic is something that fails many a programmer, best to stay away.

I see, there is no other context. The recommendation is to use generic name of any intrinsic function so that the compiler chooses the appropriate specific name to match the type of the argument, right?

F77 introduced the idea of generic functions. Then, it was just for intrinsic functions, not functions written by fortran users. Then f90 extended the idea also to user defined functions. In many cases, the result of the function is the same type and kind as its argument. Examples include things like sqrt(x), cos(x), and log(x). In other cases, the result of the function is specified by the value of one of its arguments, sometimes with a default kind if the optional argument is not present. Examples of this include things like real(n,kind=wp).

A surprising example of this second type of function is the cmplx(x,y,kind=wp) intrinsic function. For historical reasons, it is a generic function that always returns a default complex kind value unless the optional kind argument is given, and then it returns a value of that specified kind. Even experienced fortran programmers sometimes forget about that optional argument, expecting the result kind value to match that of its real input arguments. At the time cmplx() was defined, there was only one complex type, which corresponded to default REAL, so it made no sense in f77 to make its output type depend on its input type. KIND values for complex types were added in f90, after the previous convention had been established for over a decade.

Before f77 introduced the idea of generic functions, a specific function of each output type had to be defined in the standard. Up through f77, there was only one integer type, one logical type, and two floating point types, REAL and DOUBLE PRECISION, so that was not much of a burden. That is where the LOG/DLOG, SQRT/DSQRT, and so on functions were defined in fortran. Other languages followed fortran’s lead, with different functions names for each supported combination of input and output types. F90 introduced the KIND system to fortran, which is general and open ended. Now fortran can support any number of kinds, for all the different types. Imagine all the problems that would follow if each specific function had to have its own defined name, with different compilers supporting different numbers of KIND values for each of its types. There are some exceptions, of course, but other languages typically still use the old pre-f77 approach, like they are in the dark ages or something.

1 Like

@RonShepard @FortranFan Thank you for taking the time to write such detailed answers

Just curious: what’s the rationale to never using the double precision exponent? Is it just because as the name says, it doubles the precision of the default real kind, so not necessarily a 64-bit floating-point number, but could be 128-bit if 64-bit is the default real kind?

Or, there are other subtleties that would be worth knowing? :slight_smile:

Yes, that’s the main reason. When a 128-bit floating point intermediate within an expression is generated, and that is combined with other operands within the expression, the mixed-kind promotion rules in fortran require that those other operands also be converted and that the operations then be performed in the higher precision. That’s good if that is what the programmer intended, but not when the programmer really just wanted 64-bit arithmetic.

On the other hand, when the constant is written with a kind designation, 42.0_wp, then there is no question what the programmer intended and no question how the compiler will evaluate the expression.

The other aspect is that if the programmer does want to change precisions at some point, then all that is necessary is to change the value of the wp kind parameter. If the code is written appropriately, that change can be on a single line of code, and it would then be propagated throughout the whole program (which could be millions of lines). If instead the code has E and D exponents embedded throughout, then each of those exponents would need to be changed individually, and literal constants that have no exponent would likely lead to loss of precision. That was the situation fortran programmers faced prior to f90.

1 Like

Thanks a lot for all your valuable comments which have indeed made me learnt a lot. I have indeed dedicated three days refining my code however, it appears that tolerance do play a role. I tried the R script which was written as below

if(DEBUG_pH) { #this is on! it calculates pH for each time step
#---------------------------------------
#pH module

          if(ti==2)
            
          {
            
            y10f <- 0 
            y20f <- 0 
            y30f <- 0 
            y40f <- 0 
            
            x10f <- 0   
            x20f <- 0 
            x30f <- 0 
            
          }
          
          InCf <- TANproduction[ti] * (ASL/WSL) /14/2
          OutNf <- (SoilEmissionFlux[ti-1] * (PatchArea / 100))/14 + (((TANBudget[ti-1]) * (1-LOSS))/14)
          
          InNf <- TANproduction[ti] * (ASL/WSL) /14
          
          
          
          
          
          
          Tfix <- T[ti]+273.15  #K
          
          Volumefix <- WaterBudget2cm[ti] * (ASL/WSL)   #dm3
          Volumegfix <- (Porosity * PatchArea * ASL) - Volumefix   #dm3
          
          
          
          
          OutCf <- 0
          
          
          
          zH0f <- 10^(-pH[ti-1])

Conc_calc ←
function(zHV)

{
T=Tfix
V=Volumefix
Vg=Volumegfix
x10=x10f
x20=x20f
x30=x30f
y10=y10f
y20=y20f
y30=y30f
y40=y40f
InC=InCf
InN=InNf
OutN=OutNf
zHV0=zH0f
OutC=OutCf

A <- x10 + x20 + x30 + InN - OutN
B <- y10 + y20 + y30 + y40 + InC 


Kw <- 10^(-14)

K1 <- 5.67 * (10 ^ (-10)) * exp(-6286 * ((1 / T) - (1 / 298.15))) 

K2 <- (10 ^ (-1 * (((2902.39 / T) + (0.02379 * T) - 6.4980))))

K3 <- (10 ^ (-1 * (((3404.71 / T) + (0.032786 * T) - 14.8435))))

H1 <- (56 * exp(4092 * ((1 / T) - (1 / 298.15)))) / ((1 / 8.314) * 0.001 * (1 / T) * (1.013 * 10 ^ 5))

H2 <- ((1/29.41) * exp(2400 * ((1 / T) - (1 / 298.15)))) / ((1 / 8.314) * 0.001 * (1 / T) * (1.013 * 10 ^ 5))




zH <- zHV * V
vec1 <- c(0, 0, 0, 0, 0, 1, 0, -zH,	0,	0,	(1+H1*(V/Vg)),	0 			) 			#x1
vec2 <- c(0, 0, 0, 0, 1, 0, 0, K1*V,	0,	0,	H1*(V/Vg),		0			)			#x2
vec6 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)	#x3
vec3 <- c(0, 0, 1, 0, 0, 0, 0, 0, 	-zH, 	0, 	0, 			H2*(V/Vg) 		)		#y1
vec4 <- c(0, 1, 0, 0, 0, 0, 0, 0, 	K2*V,	-zH, 	0, 			H2*(V/Vg)		)		#y2
vec5 <- c(1, 0, 0, 0, 0, 0, 0, 0, 	0, 	K3*V, 0, 			(1+H2*(V/Vg))	)		#y3
vec7 <- c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0)	#y4
vec8 <- c(1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0)	#RH2
vec9 <- c(-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)	#RK3
vec10 <- c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0)	#RK2
vec11 <- c(0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0)	#RK1
vec12 <- c(0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0)	#RH1


Amatrix <- cbind(vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, vec11, vec12)
colnames(Amatrix)<- c("x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", "RK2", "RK1", "RH1")


Bvector <- c(y30, y20+InC, y10, y40-OutC, x20+InN, x10, x30-OutN, 0,	0,	0,	H1*A*(V/Vg),	H2*B*(V/Vg))	#const


solve(Amatrix,Bvector)

roots <- solve(Amatrix,Bvector)

zH0 <- zHV0*V

pH0 <- -log10(zH0/V)


pH <- -log10(zH/V)

(zH0 -InC + (-roots[9]+roots[10]+roots[11]) + (buff * (pH-pH0)))-((zH)) 
}

Blockquote

I see that when I print zHV in the R code and I didn’t specify the tolerance I see that the value of zHV changes. The author of the R code added this in the R code (related to tolerance i.e., tol in R)

Hconc[ti] ← uniroot(Conc_calc, interval=c(0,1), tol=10^-25)$root

When I set the tolerance in R I see that the result looks as below

[1] 0
[1] 1
[1] 0.5
[1] 0.25
[1] 0.125
[1] 0.0625
[1] 0.03125
[1] 0.015625
[1] 0.0078125
[1] 0.00390625
[1] 0.001953125
[1] 0.0009765625
[1] 0.0004882812
[1] 0.0002441406
[1] 0.0001220703
[1] 6.103516e-05
[1] 3.051758e-05
[1] 1.525879e-05
[1] 7.629395e-06
[1] 9.983235e-06
[1] 9.500752e-06
[1] 9.447661e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
[1] 9.447761e-06
> [1] pH = 5.024671 (i.e., -log10 H+ concentration) !which is the default result from the R code

However, in the same R code if the tolerance is not set in the R code the result is as follows

Hconc[ti] ← -uniroot(Conc_calc, interval=c(0,1))$root
pH[ti] ← -log10(Hconc[ti])

if(ti == 2 & ti3 == 1){

print(uniroot(Conc_calc, interval=c(0,1))$root)
print(pH)

[1] 0
[1] 1
[1] 0.5
[1] 0.25
[1] 0.125
[1] 0.0625
[1] 0.03125
[1] 0.015625
[1] 0.0078125
[1] 0.00390625
[1] 0.001953125
[1] 0.0009765625
[1] 0.0004882812
[1] 0.0002441406
[1] 0.0001220703
[1] 0.0001220703
[1] 0.0001220703
[1] 3.91339 (i.e., -log10 H+ concentration)

I have corrected the code implementing major comments for now and printed out variables out here is the updated code in Fortran which there was no setting of tol like in R: and what i printed.

subroutine pH_calculation
use root_module, only: root_scalar
use iso_fortran_env, only: wp => real64
character(len=5), parameter :: method = ‘brent’
integer :: ti3, ti
real(wp) :: zHV, Tf, V, Vg, x10, x20, x30, y10, y20, y30
real (wp) :: y40, InC, InCf, InN, OutN, zHV0, OutC
real(wp) :: outnf, InNf, outcf, Zh0f
real(wp) :: Tfix, Volumefix, volumegfix, x10f
real(wp) :: x20f, x30f, y10f, y20f, y30f, y40f
integer:: iflag
!real(wp), dimension(12) :: roots
real(wp) :: result
real(wp) :: roots(12)

        Do ti3=1, cv 
            Do ti = 2, tn
                if (DEBUG_pH) then ! pH module
    
                    if (ti == 2) then
    
                    y10f = 0.0_wp
    
                    y20f = 0.0_wp
    
                    y30f = 0.0_wp
    
                    y40f = 0.0_wp
    
                    x10f = 0.0_wp
    
                    x20f = 0.0_wp
    
                    x30f = 0.0_wp
                end if
    
                InCf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0_wp / 2.0_wp
                
                OutNf = (SoilEmissionFlux(ti-1, ti3) * (PatchArea_L/100.0_wp)) / 14.0_wp + (((TANBudget(ti-1, ti3)) * &
                        (1.0_wp-LOSS)) / 14.0_wp)

                InNf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0_wp

                Tfix = Ts(ti) + 273.15_wp ! K

                Volumefix = WaterBudget2cm(ti, ti3) * (ASL/WSL) ! dm3
                Volumegfix = (Porosity * PatchArea_L * ASL) - Volumefix ! dm3

                OutCf = 0.0_wp

                zH0f = 10.0_wp ** (-pH(ti-1, ti3)) !Issue with this line
                

                Tf=Tfix !To be passed to the function
                
                V=Volumefix
            
                Vg=Volumegfix
            
                x10=x10f 
            
                x20=x20f
            
                x30=x30f
            
                y10=y10f
            
                y20=y20f
            
                y30=y30f
            
                y40=y40f
            
                InC=InCf
            
                InN=InNf
            
                OutN=OutNf
            
                zHV0=zH0f
            
                OutC=OutCf
                !print *, 'Vg is',  InCf,  ti, ti3
                !print *, 'outN', zH0f, ti, ti3
                call root_scalar('bisection',Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag)
                
                Hconc(ti, ti3)= zHV ! the root of the function is the H+ concentration

                 
                !print *, 'Vg', ZHV, ti
                pH(ti, ti3)= -log10(Hconc(ti, ti3)) 
                 
                roots =Conc_calc2(zHV) ! the root of the function is the H+ concentration which the author of the R code used
                 !MOA called it roots
                y10f = roots(4)
                y20f = roots(5)
                y30f = roots(6)
                y40f = roots(7)
                x10f = roots(1)
                x20f = roots(2) 
                x30f =roots(3)
                
                !test = reshape([roots(8:12), InC], shape(test))
              
            else
                pH(ti, ti3) = constpH
                Hconc(ti, ti3) = 10.0_wp ** (-pH(ti,ti3))
            end if
       
        end do 
    end do 
    contains 
    function Conc_calc(zHV) result(result)
        use iso_fortran_env, only: wp => real64            
        IMPLICIT NONE
        real(wp), intent(in)  :: zHV
        real(wp) :: A, B, Kw, K1, K2, K3, H1, H2, zH
        real(wp) :: zH0, pH0, pHI, buff_L
        real(wp), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7
        real(wp), dimension(12) :: vec8, vec9, vec10, vec11, vec12
        real(wp) :: Amatrix(12,12), Bvector(12) 
        real(wp) :: roots(12)
        real(wp) :: result
        real(wp), parameter :: TREF=298.15_wp
        
        integer :: n=12 ! specify the size of the matrix
    
    
        A = x10 + x20 + x30 + InN - OutN
    
        B = y10 + y20 + y30 + y40 + InC
    
               
    
        Kw = 10.0_wp**(-14.0_wp)   ! Kw = 10.0**(-14.0) at 25 C
    
    
        !Print *, 'zHV', Tf, ti
                                    
        K1 = 5.67_wp * (10.0_wp**(-10.0_wp)) * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / TREF))) !defines the first H2CO3 dissociation into HCO3-
    
        K2 = (10.0_wp**(-1.0_wp * (((2902.39_wp/ Tf) + (0.02379_wp * Tf) - 6.4980_wp)))) !defines the second HCO3- dissociation into CO3--
    
        K3 = (10.0_wp**(-1.0_wp * (((3404.71_wp / Tf) + (0.032786_wp * Tf) - 14.8435_wp)))) !defines the third H2CO3-- dissociation into H2O and co2                   
    
        H1 = (56.0_wp * exp(4092.0_wp * ((1.0_wp/ Tf) - (1.0_wp / TREF)))) / ((1.0_wp / 8.314_wp) &
    
        * 0.001_wp * (1.0_wp / Tf) *  (1.013_wp * 10.0_wp ** 5.0_wp)) 
    
        H2 = ((1.0_wp / 29.41_wp) * exp(2400.0_wp * ((1.0_wp / Tf) - (1.0_wp / TREF))))  &
    
        / ((1.0_wp / 8.314_wp) * 0.001_wp * (1.0_wp / Tf) * (1.013_wp * 10.0_wp ** 5.0_wp))

                                             
    
          zH = zHV * V
       
         

        vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 
        0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1
        
        vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2
        vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
        vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1
    
        vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2
    
        vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &  
        K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3
    
        vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4
    
        vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2
    
        vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3
    
        vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2
    
        vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1
    
        vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1
    
                
    
        Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, &
    
        vec11, vec12], shape(Amatrix))
    
                        ! Set the column names
    
        !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", &
    
                        ! "RK2", "RK1", "RH1"]
    
        Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, &
        H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const
    
      
        
        call solve_system(n, Amatrix, Bvector, roots)   !I called lapack to solve the system         
    
        zH0 = zHV0 * V
        pH0 = -log10(zH0 / V)
        pHI = -log10(zH / V) !I renamed the line pH to PHI      
        result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH))
    end function Conc_calc 

!There is another function here
end subroutine

If I print the root from the fortran code here is what I obtain. could you please give a feedback.

PS C:\Users\GAG> .\my_program.exe
root is 0.0000000000000000
root is 1.0000000000000000
root is 0.50000000000000000
root is 0.25000000000000000
root is 0.12500000000000000
root is 6.2500000000000000E-002
root is 3.1250000000000000E-002
root is 1.5625000000000000E-002
root is 7.8125000000000000E-003
root is 3.9062500000000000E-003
root is 1.9531250000000000E-003
root is 9.7656250000000000E-004
root is 4.8828125000000000E-004
root is 2.4414062500000000E-004
root is 1.2207031250000000E-004
root is 6.1035156250000000E-005
root is 3.0517578125000000E-005
root is 1.5258789062500000E-005
root is 7.6293945312500000E-006
root is 3.8146972656250000E-006
root is 1.9073486328125000E-006
root is 9.5367431640625000E-007
root is 4.7683715820312500E-007
root is 2.3841857910156250E-007
root is 1.1920928955078125E-007
root is 5.9604644775390625E-008
root is 2.9802322387695312E-008
root is 1.4901161193847656E-008
root is 7.4505805969238281E-009
root is 3.7252902984619141E-009
root is 1.8626451492309570E-009
root is 9.3132257461547852E-010
root is 4.6566128730773926E-010
root is 2.3283064365386963E-010
root is 3.4924596548080444E-010
root is 2.9103830456733704E-010
root is 3.2014213502407074E-010
root is 3.3469405025243759E-010
root is 3.2741809263825417E-010
root is 3.2378011383116245E-010
root is 3.2196112442761660E-010
root is 3.2287061912938952E-010
pH is 9.4909714735041533

I think form the R code it is obvious the tolerance matter especially for a sensitive variable I am simulating? Please could you all give a feedback from my code snippet for the subroutine especially with the root_fortran module? and also any possible solution.
Thanks a lot for your help.

Blockquote

FYI root_scalar has three optional tolerance inputs. See: root_scalar – roots-fortran

1 Like

Thanks for the reply there are three tolerances indeed.

ftol !! absolute tolerance for f=0, !!
relative tol for x
atol !! absolute tol for x
there is also tol1 which is defined as tol1= eps+1.0_wp

The R own has only one called the convergence tolerance as in this documentation uniroot function - RDocumentation which is by default set to tol = .Machine$double.eps^0.25 However a user can override this default by specifying the tol like in the case of the code above the author used tol=10^-25,

which of the tolerance (ftol, atol, rtol and tol1) will perform similar operation to the R function (the uniroot function in R utilizes the brent method)?

Thanks for all the comments. Although I am new to Fortran, I have been stuck for the last 10 days trying to find root of a function.

I have been trying to try to link root fortran to my work to find the root a function. I use other libraries(such as lapack and also solved root problems with root_fortran on other codes). but I dont know why the root fortran does not work with my code, I also tried debugging. Therefore, I have resolved to writing a similar function to the uniroot function in R. Here is the link to the zeroin https://www.netlib.org/fmm/zeroin.f. but again I guess I am missing a step. could someone please help me look at this in relation to the zeroin function and provide a feedback? I will be very grateful for help

module root_finder
use iso_fortran_env, only: wp => real64
implicit none
contains

function zeroin(ax,bx,f,tol) result(root)
! Inputs
real(wp), intent(in) :: ax, bx
interface
real(wp) function f(x)
import :: wp
real(wp), intent(in) :: x
end function f
end interface
real( wp), intent(in) :: tol
real( wp) :: root
integer :: i, maxiter = 1000

    ! Local variables
    real(wp) :: a,b,c,d,e,eps,fa,fb,fc
    real(wp) :: p,q,r,s,tol1,xm
    
    ! Initialize
    eps = epsilon(1.0_wp)
    eps = eps/2.0_wp
    tol1 = 1.0_wp + eps

    ! Set brackets
    a = ax
    b = bx
    fa = f(a)
    fb = f(b)

    ! Bisection loop
    do i=1,maxiter
      c = a 
      fc = fa
      d = b - a
      e = d

      ! Update brackets
      if (abs(fc) >= abs(fb)) then
         a = b
         b = c
         c = a
         fa = fb
         fb = fc
         fc = fa
      end if

      ! Convergence test
      tol1 = 2.0_wp*eps*abs(b) + 0.5_wp*tol
      xm = 0.5_wp*(c - b)
      if (abs(xm) <= tol1 .or. fb == 0.0_wp) exit

      ! Try interpolation
      if (abs(e) < tol1 .and. abs(fa) < abs(fb)) then
        d = xm
        e = d
        else if (a /= c) then
           q = fa/fc
           r = fb/fc
           s = fb/fa
           p = s*(2.0_wp*xm*q*(q - r) - (b - a)*(r - 1.0_wp)) 
           q = (q - 1.0_wp)*(r - 1.0_wp)*(s - 1.0_wp)
           
        else 
           s = fb/fa
           p = 2.0_wp*xm*s
           q = 1.0_wp - s
         ! Use bisection
         
      end if

      ! Adjust signs and complete step
      if (p > 0.0_wp) q = -q
      p = abs(p)

      if ( (2.0_wp*p) >= (3.0_wp*xm*q - abs(tol1*q)) .and. &
           p >= abs(0.5_wp*e*q) ) then
         ! Interpolation acceptable
         e = d  
         d = p/q     
      else
         ! Use bisection
         d = xm
         e = d
      end if

      ! Update brackets
      a = b
      fa = fb
      if (abs(d) > tol1) then
         b = b + d
      else
         b = b + sign(tol1,xm)
      end if
      fb = f(b)

      if (fb*(fc/abs(fc)) > 0.0_wp) cycle

    end do
    
    ! Return result
    root = b
    
  end function
end module

@Mubaraqlanre ,

Sorry to read of the challenges you have been facing.

First, can you please provide a little bit of background on this work? Is this homework assignment or term project in a technical (chemistry) course?

According to the source code of uniroot (https://github.com/SurajGupta/r-source/blob/master/src/library/stats/man/uniroot.Rd), it just calls zeroin, which uses an absolute tolerance on the interval.

The equivalent option in roots-fortran is atol. In addition you should set rtol to 0. tol1 is not part of the public interface, it is just a variable used internally. (:warning: see edit below)

The default tolerance in R’s uniroot is relatively high. The value of eps^(1/4) is about 0.0001 for 64-bit double precision. On the other hand 10^-25 is too low. I’d expect 2*eps to be the “best” you can achieve. I’d expect this accuracy to be higher than the rest of the experimentally-measured parameters in your chemical equilibrium model.

Edit 1: I’d recommend consulting an introductory textbook on numerical methods for an explanation of terms such as machine precision and absolute or relative error. One of my favorites is Computer Methods for Mathematical Computations (1977) by Forsythe, Malcolm and Moler, published by Prentice Hall. But most newer textbooks also contain an introductory chapter on how machines represent numbers and what this means for a model analyst.

Edit 2: :warning: the statements earlier are wrong, zeroin uses a relative tolerance. The option rtol in the roots-fortran method “brent” is equivalent to tol in the classic zeroin.

1 Like

By the way, I don’t thing anyone else pointed this out, but the values in your original post are the same to 8 eight decimal places. Note that the Fortran output uses scientific notation with the appended E to signify the exponent.

If you were using the default tolerance value tol = .Machine$double.eps^0.25 or approximately 0.00012, it means your values were already “the same” effectively. Nevertheless, you should follow all the previously given advice concerning:

  • consistent usage of kind specifiers (the _wp)
  • simplified literals (1.0E-14_wp instead of 10.0**(-14))
  • transformed logarithm expressions

What strikes me as odd is in your original post you were calling root_scalar as follows:

call root_scalar(‘brent’,Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag, &
    atol = 1.0_wp, &
    Ftol = 0.10_wp**(0.25_wp))

Both the absolute tolerance of 1.0 and the function value tolerance of 0.1**0.25 ~ 0.562 appear quite high for your problem.

As you’ve said, the code is fairly complex - it contains chemistry, unit conversions, and numerics, all mixed together. It’s common to start with code like this when establishing a model and investigating it interactively in a language like R or MATLAB. Switching to a compiled language is a good opportunity to trim, refactor and document the code.

1 Like

Perhaps an idea this Community might want to brainstorm is to have periodic virtual meetings of certain duration that are support sessions for Fortranners facing thorny issues, kind of like a companion to the (erstwhile) monthly Fortran calls.

If members with specific skills with computing and user training such as @ivanpribec , @everythingfunctional , @certik , @milancurcic chime in to such support sessions, specific guidance can be provided to Fortranners with their challenges and they may gain further in their productive use of Fortran in their endeavors, such as @Mubaraqlanre here with the PhD work.

It appears OP would be keen to share the “screen” and have someone take a look at the code together given the feeling of being stuck with some issue, whether it be tolerance or some other numerical aspect. One would expect OP would be able to do so in such as support session.

This Community providing such direct, live guidance to each other might take the productive application of Fortran in scientific and technical computing to the next level, it will be back to the basics of peer-to-peer human interaction, I would hope (and prefer) it beats online forums and ChatGPT, etc.!

1 Like

I actually do already host a monthly “Office Hours” for the Fortran Users of NERSC (FUN) group.

It’s supposed to be more of a “NERSC Users Only” event, but in the case a non-NERSC user participation would be beneficial to NERSC users, I could easily justify making an exception. They’re recorded and posted publicly on the NERSC YouTube Channel anyways.

1 Like

Thanks a lot to everyone for the recommendation.

I did implement the major suggestions in the last weekend and was due to the issue of simplified literals for instance Kw= 1E-14_wp instead of Kw = 10.0**(-14.0) as mentioned by @jacobwilliams , @FortranFan @ivanpribec and many others.

I tried to isolate the subroutine, printing all the variables and checking differences. I then used sample data for both the R and fortran it seems the brent method which I was using was not converging (probably because of complexity?).
The rtol is what controls the tolerance in the brent method in root fortran as specified here in https://www.netlib.org/fmm/zeroin.f.

40 tol1 = 2.0d0epsdabs(b) + 0.5d0*tol

and in root fortran here:

I then started to try to use the brenth, brentq, bisection and Zhang methods.

All these methods work and match the R one even without specifying rtol, atol and or ftol.
The speed with which the function runs with these methods also differs and were rapid.

I am thankful for the inputs of everyone.

1 Like

Great catch! I’m very sorry for the egregious error in my previous post.

The method brent in roots-fortran or the classic zeroin, uses a relative tolerance convergence test, and not an absolute one as I mistakenly claimed earlier. I was misled by the first part of the documentation of zeroin which states:

c  tol    desired length of the interval of uncertainty of the
c         final result ( .ge. 0.0d0)

c  ... 

c      it is assumed  that   f(ax)   and   f(bx)   have  opposite  signs
c  without  a  check.  zeroin  returns a zero  x  in the given interval
c  ax,bx  to within a tolerance  4*macheps*abs(x) + tol, where macheps
c  is the relative machine precision.

I incorrectly interpreted “interval of uncertainty” as an absolute interval.

Upon more accurate reading I see uniroot makes this clearer:

uniroot() uses Fortran subroutine “zeroin” (from Netlib) based on algorithms given in the reference below. They assume a continuous function (which then is known to have at least one root in the interval).

Convergence is declared either if f(x) == 0 or the change in x for one step of the algorithm is less than tol (plus an allowance for representation error in x ).

[emphasis mine]

It’s easy to miss that the method "brent" in roots-fortran, ignores the atol parameter, when the interface offers it. (cc @jacobwilliams)

1 Like

Yes, brent in roots-fortran is basically identical to the classic the zeroin method, and only uses rtol.

1 Like