[gdal-dev] GDALDatasetRasterIO - different results between i386 and amd64

Even Rouault even.rouault at spatialys.com
Sun Aug 16 15:11:45 PDT 2015


On Sunday 16 August 2015 23:14:17 Johan Van de Wauw wrote:
> Hello,
> 
> While investigating why one of the unit tests of rasterio fails on i386
> [1], I think that the actual cause is a different result from gdal when
> using GDALDatasetRasterIO.
> 
> When opening the test file [2] I get different results when using i386.
> I could reproduce this behaviour when using only C code (attached).
> i386:
> 0 8 5 7 0 0 0 0
> 0 6 61 15 27 15 24 128
> 0 20 152 23 15 19 28 0
> 0 18 255 21 129 24 32 0
> 9 7 14 16 19 18 36 0
> 6 27 43 207 38 31 73 0
> 0 0 0 0 74 23 0 0
> 
> AMD64 and other architectures:
> 0 8 5 7 0 0 0 0
> 0 6 61 15 27 15 24 128
> 0 20 152 23 15 19 28 0
> **0 17 255 25 255 22 32 0 **
> 9 7 14 16 19 18 36 0
> 6 27 43 207 38 31 73 0
> 0 0 0 0 74 23 0 0
> 
> 
> Note that the difference only occurs on i386: on other debian
> architectures[3,] (including 32 bit like mips, arm) the same result was
> found.

Johan,

I can also reproduce. However when using -g, or -fpmath=sse -msse2 (to avoid 
i387 fpu), I get the same result as AMD64.

When adding the following debug traces:

Index: gcore/rasterio.cpp
===================================================================
--- gcore/rasterio.cpp	(révision 29662)
+++ gcore/rasterio.cpp	(copie de travail)
@@ -595,6 +595,7 @@
             dfSrcY = (iBufYOff+0.5) * dfSrcYInc + dfYOff;
             dfSrcX = 0.5 * dfSrcXInc + dfXOff;
             iSrcY = (int) dfSrcY;
+            printf("%d %d %.18g %.18g\n", iBufYOff, iSrcY, dfSrcY, dfSrcY - 
359);
 
             iBufOffset = (GPtrDiff_t)iBufYOff * (GPtrDiff_t)nLineSpace;


I get in the AMD64 case:
0 51 51.2857142857142847 -307.714285714285722
1 153 153.857142857142861 -205.142857142857139
2 256 256.428571428571445 -102.571428571428584
3 358 359 -7.10542735760100186e-15
4 461 461.571428571428555 102.571428571428555
5 564 564.14285714285711 205.142857142857139
6 666 666.714285714285666 307.714285714285722

and in the default i386 case:
0 51 51.2857142857142847 -307.714285714285722
1 153 153.857142857142861 -205.142857142857139
2 256 256.428571428571445 -102.571428571428555
3 359 359 0
4 461 461.571428571428555 102.571428571428555
5 564 564.14285714285711 205.14285714285711
6 666 666.714285714285666 307.714285714285666


Note the difference of rounding of 359 to 358 vs 359 . This must be due to i387 
computations done on 80 bits being more precise and actually correctly 
rounding, whereas AMD64 computations are done on 64 bits.

I've tried the following code (test_round_error.c) that tries to replicate the 
computation

#include <stdlib.h>
#include <stdio.h>

int main(int argc, char* argv[])
{
    int nSrcYSize = atoi(argv[1]);
    int nBufYSize = atoi(argv[2]);
    double dfYOff = atof(argv[3]);
    int nSrcYOff = atoi(argv[4]);
    double dfSrcYSize = nSrcYSize;
    double dfSrcYInc = dfSrcYSize / (double)nBufYSize;
    double d = (nSrcYOff + 0.5) * dfSrcYInc + dfYOff;
    printf("%d\n", (int)d);
    return 0;
}

and "./test_rounding_error 718 7 0.0 3" always prints 359 on AMD64 or i386, in 
-O2 or in -O0, so there's some mystery involved...

I'm not sure what could be done. Adding a small ellipson in the truncation 
might help for that case, but it could potentially move the issue to other 
cases.

Changing very slightly the test case should hopefully hide that annoyance.

Even

> 
> Kind Regards,
> Johan
> 
> [1] https://github.com/mapbox/rasterio/issues/229
> [2] https://github.com/mapbox/rasterio/blob/master/tests/data/RGB.byte.tif
> [3] https://buildd.debian.org/status/package.php?p=rasterio&suite=unstable

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com


More information about the gdal-dev mailing list