clean up duplicate spline functions, modernize syntax a bit

camb/_compilers.py

@@ -26,10 +26,10 @@ def get_ifort_version():
     return call_command("ifort -v")
 
 
-def get_gfortran_version():
-    ver = call_command("gfortran -dumpversion")
+def get_gfortran_version(command='gfortran'):
+    ver = call_command(command + " -dumpversion")
     if ver and '.' not in ver:
-        ver = call_command("gfortran -dumpfullversion")
+        ver = call_command(command + " -dumpfullversion")
     return ver
 
 
@@ -46,24 +46,34 @@ def check_gfortran(version=gfortran_min, msg=False, retry=False):
     else:
         ok = False
     if not ok and is_windows and not retry:
-        mingw = os.path.join(os.environ["ProgramFiles"], 'mingw-w64')
+        newpath = None
         best_version = gfortran_min
-        if os.path.isdir(mingw):
-            # look for mingw installation
-            dirs = [name for name in os.listdir(mingw) if
-                    gfortran_bits in name and os.path.isdir(os.path.join(mingw, name))]
-            path = None
-            for i, x in enumerate(dirs):
-                ver = x.split('-')[1]
-                if parse_version(best_version) <= parse_version(ver):
-                    best_version = ver
-                    path = x
-            if path:
-                newpath = os.path.join(mingw, path, 'mingw64', 'bin')
-                if os.path.exists(os.path.join(newpath, 'gfortran.exe')):
-                    if not compiler_environ["PATH"].startswith(newpath):
-                        compiler_environ["PATH"] = newpath + ';' + compiler_environ["PATH"]
-                    return check_gfortran(version, msg, retry=True)
+        for root in (os.path.expanduser('~'), os.environ["ProgramFiles"]):
+            mingw = os.path.join(root, 'mingw-w64')
+            if not os.path.isdir(mingw):
+                mingw = os.path.join(root, 'mingw64')
+
+            if os.path.isdir(mingw):
+                # look for mingw installation
+                dirs = [name for name in os.listdir(mingw) if
+                        gfortran_bits in name and os.path.isdir(os.path.join(mingw, name))]
+                for i, x in enumerate(dirs):
+                    if '.' in x:
+                        ver = x.split('-')[1]
+                        if parse_version(best_version) <= parse_version(ver):
+                            best_version = ver
+                            newpath = os.path.join(mingw, x, 'mingw64', 'bin')
+                bin = os.path.join(mingw, 'bin')
+                if os.path.exists(bin):
+                    ver = get_gfortran_version('"' + os.path.join(bin, 'gfortran') + '"')
+                    if ver and parse_version(best_version) <= parse_version(ver):
+                        best_version = ver
+                        newpath = bin
+        if newpath:
+            if os.path.exists(os.path.join(newpath, 'gfortran.exe')):
+                if not compiler_environ["PATH"].startswith(newpath):
+                    compiler_environ["PATH"] = newpath + ';' + compiler_environ["PATH"]
+                return check_gfortran(version, msg, retry=True)
     if ok and is_windows:
         version_str = str(subprocess.check_output("gfortran -dumpmachine", shell=True, env=compiler_environ))
         ok = gfortran_bits in version_str

fortran/DarkAge21cm.f90

@@ -1,5 +1,6 @@
     module DarkAge21cm
     use constants
+    use splines
     !Functions for collision rates.
     implicit none
     contains
@@ -30,7 +31,7 @@ function kappa_HH_21cm(T, deriv)
         allocate(logRates(n_table),logTemps(n_table),ddlogRates(n_table))
         logRates = log(real(rates,dl)*0.01**3)
         logTemps = log(real(Temps,dl))
-        call spline(logTemps,logRates,n_table,1d30,1d30,ddlogRates)
+        call spline_def(logTemps,logRates,n_table,ddlogRates)
     end if
 
     if (T<=Temps(1)) then
@@ -135,7 +136,7 @@ function kappa_pH_21cm(T, deriv) ! from astro-ph/0702487
         allocate(logRates(n_table),logTemps(n_table),ddlogRates(n_table))
         logRates = log(real(rates,dl)*factor)
         logTemps = log(real(Temps,dl))
-        call spline(logTemps,logRates,n_table,1d30,1d30,ddlogRates)
+        call spline_def(logTemps,logRates,n_table,ddlogRates)
     end if
 
     if (T<=Temps(1)) then

fortran/DarkEnergyQuintessence.f90

@@ -22,6 +22,7 @@ module Quintessence
     use results
     use constants
     use classes
+    use Interpolation
     implicit none
     private
 

fortran/SecondOrderPK.f90

@@ -176,6 +176,7 @@ end subroutine TSecondOrderPK_GetNonLinRatios_all
     subroutine GetRatios(this,CAMB_Pk, DoVel)
     !Fill the CAMB_Pk%nonlin_scaling array with sqrt(non-linear power/linear power)
     !for each redshift and wavenumber
+    use splines
     class(TSecondOrderPK) :: this
     type(MatterPowerData), target :: CAMB_Pk
     logical, intent(in):: DoVel
@@ -212,7 +213,7 @@ subroutine GetRatios(this,CAMB_Pk, DoVel)
 
             allocate(dPdLogK(CAMB_PK%num_k))
             !Get first derivative needed for series expansion at low r
-            call spline_deriv(CAMB_Pk%log_kh, CAMB_Pk%matpower(1,it), CAMB_Pk%ddmat(1,it), dPdLogK,CAMB_PK%num_k)
+            call spline_deriv(CAMB_Pk%log_kh, CAMB_Pk%matpower(:,it), CAMB_Pk%ddmat(:,it), dPdLogK,CAMB_PK%num_k)
 
             do i=1, CAMB_PK%num_k
 

fortran/SeparableBispectrum.f90

@@ -20,6 +20,7 @@ module Bispectrum
     use SpherBessels
     use constants
     use MpiUtils
+    use splines
     implicit none
 
     integer, parameter :: max_bispectrum_deltas = 5, max_bispectrum_fields=3
@@ -78,15 +79,15 @@ subroutine InitBesselDerivs(CTrans)
     do i=1, CTrans%ls%nl
         !Spline agrees well
         !  call spline_deriv(BessRanges%points,ajl(1,i),ajlpr(1,i),dJl(1,i),BessRanges%npoints)
-        !  call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))
+        !  call spline_def(BessRanges%points,dJl(1,i),BessRanges%npoints,dddJl(1,i))
 
         l1 = CTrans%ls%l(i)
         do j=1, BessRanges%npoints
             call BJL(l1-1,BessRanges%points(j),Jm)
             call BJL(l1+1,BessRanges%points(j),Jp)
             dJl(j,i)=  ( l1*Jm - (l1+1)*Jp)/(2*l1+1)
         end do
-        call spline(BessRanges%points,dJl(1,i),BessRanges%npoints,spl_large,spl_large,dddJl(1,i))
+        call spline_def(BessRanges%points,dJl(:,i),BessRanges%npoints,dddJl(:,i))
 
     end do
 

fortran/bessels.f90

@@ -11,6 +11,7 @@ module SpherBessels
     use results
     use RangeUtils
     use MpiUtils
+    use splines
     implicit none
     private
 
@@ -115,8 +116,8 @@ subroutine GenerateBessels(lSamp, CP)
         end do
 
         !     get the interpolation matrix for bessel functions
-        call spline(BessRanges%points,ajl(1,j),num_xx,spl_large,spl_large,ajlpr(1,j))
-        call spline(BessRanges%points,ajlpr(1,j),num_xx,spl_large,spl_large,ddajlpr(1,j))
+        call spline_def(BessRanges%points,ajl(:,j),num_xx,ajlpr(:,j))
+        call spline_def(BessRanges%points,ajlpr(:,j),num_xx,ddajlpr(:,j))
     end do
     !$OMP END PARALLEL DO
 

fortran/camb_python.f90

@@ -8,6 +8,7 @@ module handles
     use DarkEnergyPPF
     use ObjectLists
     use classes
+    use Interpolation
     implicit none
 
     Type c_MatterTransferData
@@ -284,10 +285,9 @@ end subroutine CAMBdata_GetSigma8
 
     subroutine CAMBdata_GetSigmaRArray(State, sigma, R, nR, z_ix, nz, var1, var2)
     Type(CAMBdata) :: State
+    integer, intent(in) :: nR, nz  
     real(dl) :: sigma(nR,nz)
-    integer :: nR
     real(dl) :: R(nR)
-    integer nz
     integer var1, var2, z_ix(nz)
     integer i, ix
 
@@ -303,8 +303,8 @@ end subroutine CAMBdata_GetSigmaRArray
 
     subroutine CAMBdata_GetMatterTransferks(State, nk, ks)
     Type(CAMBdata) :: State
-    real(dl) ks(nk)
     integer nk
+    real(dl) ks(nk)
 
     if (nk/=0) then
         ks =  State%MT%TransferData(Transfer_kh,:,1) * (State%CP%H0/100)
@@ -677,6 +677,7 @@ end function CAMB_TimeEvolution
 
 
     subroutine GetBackgroundThermalEvolution(this, ntimes, times, outputs)
+    use splines
     Type(CAMBdata) :: this
     integer, intent(in) :: ntimes
     real(dl), intent(in) :: times(ntimes)

fortran/cmbmain.f90

@@ -754,7 +754,8 @@ subroutine GetSourceMem
         deallocate(ThisSources%LinearSrc)
     allocate(ThisSources%LinearSrc(ThisSources%Evolve_q%npoints,&
         ThisSources%SourceNum,State%TimeSteps%npoints), source=0._dl, stat=err)
-    if (err/=0) call GlobalError('Sources requires too much memory to allocate', error_unsupported_params)                                                                               
+    if (err/=0) call GlobalError('Sources requires too much memory to allocate', &
+        error_unsupported_params)                                                                               
 
     end subroutine GetSourceMem
 
@@ -1181,8 +1182,7 @@ subroutine MakeNonlinearSources
             !Do not use an associate for scaling. It does not work.
             scaling = State%CAMB_Pk%nonlin_ratio(ik,1:State%num_transfer_redshifts)
             if (all(abs(scaling-1) < 5e-4)) cycle
-            call spline(State%Transfer_Times(1), scaling(1), State%num_transfer_redshifts,&
-                spl_large, spl_large, ddScaling(1))
+            call spline_def(State%Transfer_Times, scaling, State%num_transfer_redshifts,ddScaling)
 
             tf_lo=1
             tf_hi=tf_lo+1
@@ -1222,8 +1222,8 @@ subroutine InitSourceInterpolation
     !$OMP PARALLEL DO DEFAULT(SHARED), SCHEDULE(STATIC), PRIVATE(i,j)
     do  i=1,State%TimeSteps%npoints
         do j=1, ThisSources%SourceNum
-            call spline(ThisSources%Evolve_q%points,ScaledSrc(1,j,i), &
-                ThisSources%Evolve_q%npoints,spl_large, spl_large, ddScaledSrc(1,j,i))
+            call spline_def(ThisSources%Evolve_q%points,ScaledSrc(:,j,i), &
+                ThisSources%Evolve_q%npoints, ddScaledSrc(:,j,i))
         end do
     end do
     !$OMP END PARALLEL DO
@@ -1374,8 +1374,8 @@ subroutine InterpolateSources(IV)
 
     if (.not.State%flat) then
         do i=1, ThisSources%SourceNum
-            call spline(State%TimeSteps%points,IV%Source_q(1,i),State%TimeSteps%npoints,&
-                spl_large,spl_large,IV%ddSource_q(1,i))
+            call spline_def(State%TimeSteps%points,IV%Source_q(:,i),State%TimeSteps%npoints,&
+                IV%ddSource_q(:,i))
         end do
     end if
 

fortran/config.f90

@@ -41,9 +41,6 @@ module config
 
     real(dl), parameter :: tol=1.0d-4 !Base tolerance for perturbation integrations
 
-    !     used as parameter for spline - tells it to use 'natural' end values
-    real(dl), parameter :: spl_large=1.e40_dl
-
     character(LEN=1024) :: highL_unlensed_cl_template = 'HighLExtrapTemplate_lenspotentialCls.dat'
     !fiducial high-accuracy high-L C_L used for making small cosmology-independent numerical corrections
     !to lensing and C_L interpolation. Ideally close to models of interest, but dependence is weak.

fortran/lensing.f90

@@ -39,6 +39,7 @@ module lensing
     use Precision
     use results
     use constants, only : const_pi, const_twopi, const_fourpi
+    use splines
     implicit none
     integer, parameter :: lensing_method_curv_corr=1,lensing_method_flat_corr=2, &
         lensing_method_harmonic=3
@@ -462,7 +463,7 @@ subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
             if (.false.) then
                 if (abs(sum(corrcontribs(1:jmax,1)))>1e-11) print *,i,sum(corrcontribs(1:jmax,1))
                 do j=1,4
-                    call spline(xl,corrcontribs(1,j),jmax,1d30,1d30,ddcontribs(1,j))
+                    call spline_def(xl,corrcontribs(:,j),jmax,ddcontribs(:,j))
                 end do
                 corr=0
                 llo=1
@@ -1157,7 +1158,7 @@ subroutine GetBessels(MaxArg)
             Bess(i,:)=Bessel_jn(0, maxbessel,x(i))
         end do
         do ix=0,maxbessel
-            call spline(x,Bess(1,ix),max_bes_ix,spl_large,spl_large,ddBess(1,ix))
+            call spline_def(x,Bess(:,ix),max_bes_ix,ddBess(:,ix))
         end do
 
         deallocate(x)

fortran/massive_neutrinos.f90

@@ -108,6 +108,7 @@ subroutine TNuPerturbations_init(this,Accuracy)
     end subroutine TNuPerturbations_init
 
     subroutine ThermalNuBackground_init(this)
+    use splines
     class(TThermalNuBackground) :: this
     !  Initialize interpolation tables for massive neutrinos.
     !  Use cubic splines interpolation of log rhonu and pnu vs. log a*m.
@@ -141,6 +142,7 @@ end subroutine ThermalNuBackground_init
     subroutine nuRhoPres(am,rhonu,pnu)
     !  Compute the density and pressure of one eigenstate of massive neutrinos,
     !  in units of the mean density of one flavor of massless neutrinos.
+    use splines
     real(dl),  parameter :: qmax=30._dl
     integer, parameter :: nq=100
     real(dl) dum1(nq+1),dum2(nq+1)

fortran/recfast.f90

@@ -501,7 +501,7 @@ subroutine TRecfast_init(this,State, WantTSpin)
     implicit none
     class(TRecfast), target :: this
     class(TCAMBdata), target :: State
-    real(dl) :: Trad,Tmat,Tspin,d0hi,d0lo
+    real(dl) :: Trad,Tmat,Tspin
     integer :: I
     Type(RecombinationData), pointer :: Calc
     logical, intent(in), optional :: WantTSpin
@@ -701,12 +701,10 @@ subroutine TRecfast_init(this,State, WantTSpin)
             !          write (*,'(5E15.5)') zend, Trad, Tmat, Tspin, x
         end do
 
-        d0hi=1.0d40
-        d0lo=1.0d40
-        call spline(Calc%zrec,Calc%xrec,nz,d0lo,d0hi,Calc%dxrec)
-        call spline(Calc%zrec,Calc%tmrec,nz,d0lo,d0hi,Calc%dtmrec)
+        call spline_def(Calc%zrec,Calc%xrec,nz,Calc%dxrec)
+        call spline_def(Calc%zrec,Calc%tmrec,nz,Calc%dtmrec)
         if (Calc%doTspin) then
-            call spline(Calc%zrec,Calc%tsrec,nz,d0lo,d0hi,Calc%dtsrec)
+            call spline_def(Calc%zrec,Calc%tsrec,nz,Calc%dtsrec)
         end if
     class default
         call MpiStop('Wrong state type')