Something's buggy with array assignment .=?

General Usage Performance
1 
function testnc(
    init::Dict,eroot::Dict;
    regID::AbstractString="GLB", timeID::Union{Integer,Vector}=0,
    gres::Real=0
)

    emod,epar,ereg,etime = erainitialize(
        init;
        modID="dsfc",parID="p_sfc",regID=regID,timeID=timeID,
        gres=gres
    ); ehr = hrindy(emod); nlon,nlat = ereg["size"];

    ps1  = zeros(Float32,nlon,nlat)
    dvec = collect(Date(etime["Begin"],1):Month(1):Date(etime["End"],12));


    for dtii in dvec

        nhr = ehr * daysinmonth(dtii);
        pds,pvar = erarawread(emod,epar,ereg,eroot,dtii); p = pvar[:]*1; close(pds)
        @info "Complete Dataset   Extraction | p = 1000hPa : $(sum(p.==100000))"
        psum1 = 0; psum2 = 0

        for it = 1 : nhr

            ps1 .= p[:,:,it]*1; psum1 += sum(ps1.==100000)
            ps2  = p[:,:,it]*1; psum2 += sum(ps2.==100000)
            ps3  = ps2 .- ps1;  psum3 += sum(ps3.!=0)

        end

        @info "Time-by-Time Alloc Extraction | p = 1000hPa : $(psum1)"
        @info "Time-by-Time Slice Extraction | p = 1000hPa : $(psum2)"
        @info "Comparing Extraction Methods  | p1 =/= p2   : $(psum3)"

    end

end

init,eroot = erastartup(aID=2,dID=1,path="/n/kuangdss01/lab/",welcome=false)
testnc(init,eroot,regID="SEA",timeID=1980)

returns this (I’m doing this month by month so the loop repeats 12 times over the year):

[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31572272
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 808
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 29538934
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31538320
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 30557518
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31576102
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 30557518
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31576102
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31576102
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 30557518
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31576102
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 30557518
[ Info: Complete Dataset   Extraction | ps = 1000hPa : 0
[ Info: Time-by-Time Alloc Extraction | p1 = 1000hPa : 0
[ Info: Time-by-Time Slice Extraction | p2 = 1000hPa : 0
[ Info: Comparing Extraction Methods  | p1 =/= p2   : 31576102

I would basically like to ask why ps1 and ps2 are different even though they basically are extracting from the same 3D array p

julia> versioninfo()
Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: Linux (x86_64-redhat-linux)
  CPU: Intel(R) Xeon(R) Platinum 8268 CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_PATH = /n/helmod/apps/centos7/Core/julia/1.3.1-fasrc01/bin
  JULIA_INCLUDE = /n/helmod/apps/centos7/Core/julia/1.3.1-fasrc01/include
  JULIA_DEPOT_PATH = /n/holyscratch01/kuang_lab/nwong/.julia/:/n/holyscratch01/kuang_lab/nwong/.julia/:
  JULIA_LIB = /n/helmod/apps/centos7/Core/julia/1.3.1-fasrc01/lib
  JULIA_HOME = /n/helmod/apps/centos7/Core/julia/1.3.1-fasrc01
2

Have you tried to check approximate equality instead? \approx

3

Approximate equality gives the same result. However, I did a check on the absolute error, and the absolute error is very small (about ~0.007 of an average of 100000), and the maximum error is around 2-3 or thereabouts.

This is still very weird. It’s probably okay for most of my results but this is something I need to investigate further.

4

Actually I figured out why.

This is a Float64 to Float32 conversion error. The array was Float64 and then allocated to a Float32!

Ahhhhhhhhhhhhhhh.

5

Pretty much why I don’t use anything other than Float64 (for floats). Unless you are sending that to a GPU, probably not worth your time thinking about it.

6

Certainly not worth thinking about if its premature optimization, but even on CPU, Float32 computation can often be faster, particularly (obviously) for memory-bound calculations.

7

If the code is SIMD, you can fit twice as many Float32 as Float64 into the CPU’s vector units.
But I haven’t found it to be worth the loss in accuracy; I encountered many issues when trying it before, such as matrices that were comfortably positive define in Float64 being singular with Float32.