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.
2 Likes
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:
1 Like
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? 
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 moduleif(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=OutCfA <- 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
containsfunction 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
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. (
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: 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_wpinstead of10.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) == 0or the change inxfor one step of the algorithm is less thantol(plus an allowance for representation error inx).[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