Precision Error.

[Manan_G_]

Hello,

I have been using MATLAB for quite some years now but am totally new to Fortran. I am getting some precision errors which I never encountered in MATLAB. The code has the Plank's Constant, h = 6.626068e-34, which has to be cubed. There are multiple files in the whole project where this operation is used. Ideally, the cubed answer must be a very small number. But I am getting an answer of 0.0000000000E+000, because of which all my further calculations result in NaN.

I have tried all sorts of different specifiers like REAL, DOUBLE PRECISION but I still get the same error. I would be really appreciate if you can help me. Thank you in advance.

[TimP]

We can't guess in what way you have tried double precision, e.g. the following should be OK: double precision result, h result = 6.626068d-34 ** 3 ! or even h = 6.626068e-34 result = h ** 3 Note that you would need to promote a real to double before operating on it, to gain the additional range. matlab would default to double throughout; you would need to make double precision explicit where you need it in Fortran.

[Manan_G_]

Hi TimP, Thank you for your quick response. I really appreciate it. :) However, the problem persists. If I use 6.626068d-34 ** 3 directly, I get a correct answer. But if I use h ** 3, I get the same error. Why do you think is this happening? Also, I must inform you that the variable 'h' as been declared common to the entire project. I do not know what the default data type it assumes. Is that whats causing the problem? Also, I have a similar problem with other expressions, where all the variables have been declared as common. here is the code: subroutine sc_table() include 'scattering_table_parameters.h' common &/pi/pi,two_pi &/temp/tem,Vt &/fund_const/q,h,kb,am0,eps_0 &/mass/amd,amc,smh,hhm,tm(3),hm(3) &/nonp/af,af2,af4 &/scatt_par/emax,de,w(n_scatt_max,n_regions_scatt), 1 tau_max(n_latt,n_latt,n_regions_scatt), 2 max_scatt_mech(n_regions_scatt) &/lattice_temp_par/tem_latt_max,dlatt &/scatt_par2/flag_mech(n_scatt_max,n_regions_scatt) &/table/ 1scatt_table(n_lev,n_latt,n_latt,n_scatt_max,n_regions_scatt) &/dielec_func/eps_sc,eps_oxide,gamma &/density/density,sound_velocity &/select_acouctic/acoustic_scattering &/select_Coulomb/Coulomb_scattering &/select_intervalley_g1/intervalley_TA_g &/select_intervalley_g2/intervalley_LA_g &/select_intervalley_g3/intervalley_LO_g &/select_intervalley_f1/intervalley_TA_f &/select_intervalley_f2/intervalley_LO_f &/select_intervalley_f3/intervalley_TO_f &/select_intervalley_3/intervalley_first_g &/select_intervalley_4/intervalley_first_f &/sigma_acoustic/sigma_acoustic &/doping_regions/n_region &/coulomb/doping_density(n_doping_regions), 1 Energy_debye(n_latt,n_regions_scatt) &/Deformation_potential/DefPot(6) &/energy_intervalley/wq_phonon(6) &/dispersion_function/w0_phonon(6),sound_phonon(6),c_phonon(6) &/acoustic/sigma &/intervalley1/coupling_constant &/intervalley2/delta_fi,final_valleys,i_mech integer n_region character*30 out_file_1, out_file_2 real kb if(intervalley_LO_g.eq.1)then wq = wq_phonon(3) coupling_constant = DefPot(3) w0=w0_phonon(3) vs=sound_phonon(3) c=c_phonon(3) delta_fi = 0. final_valleys = 1. i_mech = 1 out_file_1 = 'intervalley_LO_g_ab' out_file_2 = 'intervalley_LO_g_em' C call intervalley_LO_phonon(i_count,i_g,wq,w0, C 1 vs,c,out_file_1,out_file_2) call intervalley_LO_dispersion(i_count,i_g,wq,w0, 1 vs,c,out_file_1,out_file_2) endif subroutine intervalley_LO_dispersion(i_count,i_region,wq,w0, 1 vs,c, out_file_1,out_file_2) include 'scattering_table_parameters.h' common &/pi/pi,two_pi &/temp/tem,Vt &/fund_const/q,h,kb,am0,eps_0 &/dri/qh &/mass/amd,amc,smh,hhm,tm(3),hm(3) &/nonp/af,af2,af4 &/scatt_par/emax,de,w(n_scatt_max,n_regions_scatt), 1 tau_max(n_latt,n_latt,n_regions_scatt), 2 max_scatt_mech(n_regions_scatt) &/lattice_temp_par/tem_latt_max,dlatt &/scatt_par2/flag_mech(n_scatt_max,n_regions_scatt) &/table/ 1scatt_table(n_lev,n_latt,n_latt,n_scatt_max,n_regions_scatt) &/density/density,sound_velocity &/intervalley1/coupling_constant &/intervalley2/delta_fi,final_valleys,i_mech &/temp_param/band_gap(n_latt),af_temp(n_latt), 1 af2_temp(n_latt),af4_temp(n_latt) integer i_count,i_region real kb, wq,w0,vs,c real k_vector DOUBLE PRECISION tau_const_nom, tau_const_denom DOUBLE PRECISION q_log_ab_nom, q_log_ab_denom DOUBLE PRECISION q_log_em_nom, q_log_em_denom character*30 out_file_1, out_file_2 C Defining some constants for absorption and emission tau_const_nom = final_valleys*amc*((coupling_constant*q)**2.) tau_const_denom=8*pi*(h**2)*density*abs(c) tau_const=tau_const_nom/tau_const_denom PRINT *, 'tau_const_nom = ',tau_const_nom A=1+2.*amc*abs(c)/h E=1-2.*amc*abs(c)/h q_const_ab=smh/A !smh=sqrt(2*amc*q)/h q_const_em=smh/E w0_ekv_ab=(w0/qh)*A !energy in eV w0_ekv_em=(w0/qh)*E C (a) Scattering rate - absorption i_count = i_count + 1 C open(unit=10, file=out_file_1, status='unknown') C write(10,*)'energy ',out_file_1 C **************************************************** C * changing the optical phonon and the lattice temeperature C **************************************************** do k=1,n_latt !changing of the LO phonon temperature t_lo=dlatt*float(k-1) Vt1=Vt*(1+t_lo/tem) rnq = 1./(exp(wq/Vt1)-1.) do j = 1, n_latt t_latt=dlatt*float(j-1) do i = 1, n_lev ee = de*float(i) fe = ee*(1.+af_temp(j)*ee) C k_vector=smh*sqrt(fe) k_vector=smh*sqrt(ee) if(k_vector.le.0)then absorption = 0. else factor_max_ab = sqrt(ee+w0_ekv_ab)+sqrt(ee) factor_min_ab = sqrt(ee+w0_ekv_ab)-sqrt(ee) q_max_ab=q_const_ab*factor_max_ab q_min_ab=q_const_ab*factor_min_ab q_log_ab_nom=q_max_ab**2.-w0/abs(c) q_log_ab_denom=q_min_ab**2.-w0/abs(c) absorption = tau_const*rnq 1 *abs(log(q_log_ab_nom/q_log_ab_denom))/k_vector endif scatt_table(i,j,k,i_count,i_region) = absorption ! PRINT *, 'Intervalley LO Abs = ',scatt_table(i,j,k, ! 1 i_count,i_region) C write(10,*)ee,(t_latt+tem),(t_lo+tem),absorption enddo enddo enddo C close(10) flag_mech(i_count,i_region) = i_mech w(i_count,i_region) = wq - delta_fi C (b) Scattering rate - emission i_count = i_count + 1 C open(unit=11, file=out_file_2, status='unknown') C write(11,*)'energy ',out_file_2 C **************************************************** C * changing the lattice temeperature C **************************************************** do k=1,n_latt !changing of the LO phonon temperature t_lo=dlatt*float(k-1) Vt1=Vt*(1+t_lo/tem) rnq = 1./(exp(wq/Vt1)-1.) do j = 1, n_latt t_latt=dlatt*float(j-1) do i = 1, n_lev ee = de*float(i) fe = ee*(1.+af_temp(j)*ee) C k_vector=smh*sqrt(fe) k_vector=smh*sqrt(ee) if((ee-w0_ekv_em).le.0)then emission = 0. else factor_max_em = sqrt(ee-w0_ekv_em)+sqrt(ee) factor_min_em = sqrt(ee)-sqrt(ee-w0_ekv_em) q_max_em=q_const_em*factor_max_em q_min_em=q_const_em*factor_min_em q_log_em_nom=q_max_em**2.-w0/abs(c) q_log_em_denom=q_min_em**2.-w0/abs(c) emission = tau_const*(rnq+1.) 1 *abs(log(q_log_em_nom/q_log_em_denom))/k_vector endif scatt_table(i,j,k,i_count,i_region) = emission C write(11,*)ee,(t_latt+tem),(t_lo+tem),emission enddo enddo enddo C close(11) flag_mech(i_count,i_region) = i_mech w(i_count,i_region) = - wq - delta_fi return end -------------------------------------------------------------------------------------------------------- All the variables are initialized in mat_par_initialization.f and I have attached it. The problem occurs in 'tau_const_nom' and 'tau_const_denom' which give a result of 0.00000000000000000E+000, because of my calculations go wrong. I have tried all possible ways that I could think of to solve this. Can you predict why does this happen? I would really appreciate it if you can help me. Thank you in advance.

[TimP]

It looks as if you let all your common variables go to default typing, which would be real(kind=4) unless you set a compiler option to change default precision. I suppose the compiler's map list option should help to confirm it.

[John_Campbell]

Manam G, . What hardware environment was used to develop this software ? It would appear that it might not be a default 32-bit real. You need to use 64-bit (8-byte) reals to achieve what you want, which appears to have been the default for the initial development environment. If you select default 8-byte reals, you may also need to consider the precision of all constants declared. Even h**2 fails for default reals. . You are ignoring the precision and exponent range available from the real variable kind you are using. Read the documentation of real KIND and SELECTED_REAL_KIND. . I wrote a simple program to demonstrate the overflow of "Result" for different real precisions available with ifort and most other win32/64 compilers. This should demonstrate the problem you are having. ! test real Kind ! integer (kind(1)) :: i do i = 0,40 write (*,*) i, selected_real_kind(i,1) end do ! call h_4 call h_8 call h_16 end subroutine h_4 ! real (kind=selected_real_kind(6,3)) :: h, result ! p = 1 to 6 integer (kind(1)) :: i ! write (*,*) ' Kind Test =',kind(h), ' huge=',huge(h) do i = 1,3 h = 6.626068e-34 result = h ** I write (*,*) ' i=',i, ' h =',h, ' result=',result end do end subroutine h_8 ! real (kind=selected_real_kind(12,3)) :: h, result ! p = 7 to 15 integer (kind(1)) :: i ! write (*,*) ' Kind Test =',kind(h), ' huge=',huge(h) do i = 1,3 h = 6.626068e-34 result = h ** I write (*,*) ' i=',i, ' h =',h, ' result=',result end do end subroutine h_16 ! real (kind=selected_real_kind(18,3)) :: h, result ! p = 16 + integer (kind(1)) :: i ! write (*,*) ' Kind Test =',kind(h), ' huge=',huge(h) do i = 1,3 h = 6.626068e-34 result = h ** I write (*,*) ' i=',i, ' h =',h, ' result=',result end do end

[Les_Neilson]

Manan Just as a comment on the coding style. You say you are new to Fortran - is this your own code or some you have inhertited? I say this because it is written in F77 style. If this is your own code then I strongly reccommend you use the newer f95+ style. Use modules instead of common. Use implicit none in your code Variables such as pi, two_pi, and h itself do not change value and so should be declared as named constants (parameter) e.g. have a seperate module where you declare all of the scientific constants and then USE that module throughout your code There are other suggestions such as having a KINDS module where you declare three KINDS : dp and sp for double and single precision and wp for working precision where you set wp to one of either sp or dp. Use kind wp in your code then when you want to change precision much of your changes is limited to re-assigning wp to the required kind. Les

[Manan_G_]

Hi Everyone, My problem has been solved. Thank you very much TimP, John and Les. Les: I have inherited this code and will keep your suggestion in mind. Thank you once again everyone. I really appreciate you taking time out to help me. Cheers..!! :) :)

[John_Campbell]

I would recommend that you further investigate the problem you are having with numerical underflow and more generally with precision. Possibly the code was originally developed on a system where 32bit reals were not the default, but probably a 60 or 64 bit real. You should check if 32bit reals give the required precision for the calculations being performed. If this is not the case, I would recommend that you adopt 64 bit reals (real*8) for all your variables and constants. (Named parameters are a good way to convert constants). I have found that explicitly declaring all variables and the use of IMPLICIT NONE helps when implementing this change. I would not resort to using complier options, as this is not a robust form of documentation. Les has also commented on the coding style. To me, the use of COMMON in this way is unusual and an alternative approach of common in INCLUDE files or preferably converted to MODULE would be more robust. (I have never seen so many named common in a single continued statement, although it must be legal) I hope these comment help. . John