Double precision in compute shader (full code)

Question

I send double precision data to my compute shader using an ssbo but somehow looses precision somewhere in the shader. I posted this problem here but has not yet received a response.

I read a position which contains 4 doubles that I use in calculations. If I in the beginning of the shader write the position to the buffer the values are correct. If I write to the buffer after my calculations I get a precision loss. I never change the values in the position. In my test I can even change "true" to "variable == variable" and loose precision only on that. Some of the code is shown in the post i linked to.

I had a similar problem with my constants which was solved with the LF suffix. Is there a similar fix for variables?

I have tried to isolate the problem, here is my code:

My .csh-file:

#version 430

struct aVec4 {
  double x, y, z, w;
};

struct FT {
  aVec4 force;
  aVec4 torque;
};

struct Triangle {
  aVec4 position;
  aVec4 normal;
  aVec4 depth;
};

layout( std430, binding = 0 ) buffer abc { Triangle bs[]; };
layout( std430, binding = 1 ) buffer def { FT forces[]; };
layout( local_size_x = 1 ) in;

dvec3 toVec3( in aVec4 vector ) {
  dvec3 v = dvec3( vector.x, vector.y, vector.z );
  return v;
}

aVec4 function1( inout Triangle triangles[2], inout int numTriangles ) {
  Triangle anotherTriangle = bs[gl_GlobalInvocationID.x];
  int numBelow = 1;

  if ( anotherTriangle.depth.x < 0.0LF ) {
    numBelow++;
  }

  if ( numBelow == numBelow ) {
    numTriangles = 1;

    aVec4 result = anotherTriangle.position;
    return result;
  }
}

FT function2( in Triangle triangles[2], in int numTriangles ) {
  FT ft;
  dvec3 force = dvec3( 0.0LF, 0.0LF, 0.0LF );

  for ( int i = 0; i < numTriangles; i++ ){
    dvec3 normal  = toVec3( triangles[i].normal );
    force  += - normal;
  }

  ft.force = triangles[0].position;
  return ft;
}

void main()
{
  Triangle triangles[2];
  triangles[0] = bs[gl_GlobalInvocationID.x];
  int numTriangles = 1;

  aVec4 result = function1( triangles, numTriangles );
  FT ft = function2( triangles, numTriangles );
  ft.torque = result;

  forces[gl_GlobalInvocationID.x] = ft;
}

My .cpp-file:

#include <GL/glew.h>
#include <GL/freeglut.h>

#include <vector>

#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>

#include <cmath>

class Vec3
{
public:
  Vec3()
    : x(0)
    , y(0)
    , z(0)
    , w(0)
  {
  }

  Vec3(double x, double y, double z)
    : x(x)
    , y(y)
    , z(z)
    , w(0)
  {
  }

  double x, y, z, w;
};


std::ostream& operator<<(std::ostream& stream, const Vec3& v)
{
  stream << std::setprecision(20) << std::setw(20) << std::scientific;
  stream << "Vec3(" <<
    std::setw(23) << v.x << ", " <<
    std::setw(23) << v.y << ", " <<
    std::setw(23) << v.z << ") (" <<
    std::setw(23) << v.w << ")";
  return stream;
}


GLuint loadShader( std::string fileName )
{
  GLuint program = glCreateProgram( );
  GLuint shader = glCreateShader( GL_COMPUTE_SHADER );

  std::ifstream shaderStream(fileName.c_str(), std::ios::in );
  if (!shaderStream.is_open())
  {
    std::cout << "Could not open shader source file '" << fileName << "'." << std::endl;
    exit(1);
  }

  std::string src;
  std::string line = "";
  while ( getline( shaderStream, line ) )
    src += "\n" + line;

  char const* srcPointer = src.c_str( );
  glShaderSource( shader, 1, &srcPointer, NULL );
  glCompileShader( shader );

  int rValue;
  glGetShaderiv( shader, GL_COMPILE_STATUS, &rValue );
  if (!rValue)
  {
    std::cout << "Unable to compile compute shader" << std::endl;
    GLchar log[ 10240 ];
    GLsizei length;
    glGetShaderInfoLog( shader, 10239, &length, log );
    fprintf( stderr, "Compiler log:\n%s\n", log );
    exit( 40 );
  }

  glAttachShader(program, shader);
  glLinkProgram(program);

  glGetProgramiv( program, GL_LINK_STATUS, &rValue );
  if (!rValue)
  {
    std::cout << "Unable to link compute shader program" << std::endl;
    GLchar log[ 10240 ];
    GLsizei length;
    glGetProgramInfoLog( program, 10239, &length, log );
    fprintf( stderr, "Linker log:\n%s\n", log );
    exit( 41 );
  }

  glDeleteShader( shader );

  return program;
}

struct ShaderData
{
  Vec3 position;
  Vec3 normal;
  Vec3 depth;
};



int main( int argc, char **argv ) {

  glutInit( &argc, argv );
  glutCreateWindow( "TEST1" ); //TODO: Why?

  glewInit( );

  std::string shaderFileName = "tutorial_computeShaderBugSearch.csh";
  GLuint program = loadShader( shaderFileName );

  GLuint ssbo;
  GLuint result;

  glGenBuffers( 1, &ssbo );
  glGenBuffers( 1, &result );

  glUseProgram( program );

  std::vector< ShaderData > buf;
  buf.resize( 1 );

  ShaderData shaderData;
  shaderData.position = Vec3( double(1) / double(3), 0, 0 );
  shaderData.normal = Vec3( 1, 1, 1 );
  shaderData.depth.x = 1;
  buf[ 0 ] = shaderData;

  glBindBuffer( GL_SHADER_STORAGE_BUFFER, ssbo );
  glBufferData( GL_SHADER_STORAGE_BUFFER, buf.size( ) * sizeof(ShaderData), buf.data( ), GL_STATIC_DRAW );

  std::vector< Vec3 > forces;
  forces.resize( 2 ); // 1 for the forces + 1 for the torques

  glBindBuffer( GL_SHADER_STORAGE_BUFFER, result );
  glBufferData( GL_SHADER_STORAGE_BUFFER, forces.size( ) * sizeof(Vec3), forces.data( ), GL_STATIC_DRAW );

  glUseProgram( program );
  glBindBufferBase( GL_SHADER_STORAGE_BUFFER, 0, ssbo );
  glBindBufferBase( GL_SHADER_STORAGE_BUFFER, 1, result );
  glDispatchCompute( 1, 1, 1 );
  glMemoryBarrier( GL_VERTEX_ATTRIB_ARRAY_BARRIER_BIT );

  glBindBuffer( GL_SHADER_STORAGE_BUFFER, result );
  Vec3* ptr = (Vec3*)glMapBuffer( GL_SHADER_STORAGE_BUFFER, GL_READ_ONLY );

  Vec3 a = ptr[ 0 ];
  Vec3 b = ptr[ 1 ];

  std::cout << "a: " << std::setprecision(20) << a.x << std::endl;
  std::cout << "b: " << std::setprecision(20) << b.x << std::endl;
  std::cout << "diff: " << std::abs(a.x-b.x) << std::endl;
}

This produces the output:

a: 0.3333333432674408
b: 0.33333333333333331
diff: 9.9341074810688212e-099