Вычислительные шейдеры в OpenGL

Как известно, современные GPU являются мощными массивно-параллельными вычислительными устройствами, быстродействие которых заметно превосходит быстродействие современных CPU. Для исользования этих возможностей GPU в неграфических задачах существуют специалтизированные API - CUDA и OpenCL.

Это мощные и гибкие API, дающие в руки программиста доступ к вычислитльным возможностям GPU, но в некотрых случаях их использование является overkill-ом. Часто возникает задача выполнения определенных (и не очень сложных) вычислений, напрямую связанных с рендерингом данных. И было бы очень удоюбно, если бы эти вычисления можно было выполнить на GPU средствами самого OpenGL.

Именно для этих целей и было введено расширение ARB_compute_shaders, входящее в состав OpenGL 4.3. Это расширение вводит в OpenGL новый тип шейдеров - вычислительные (compute). Шейдеры этого типа не входят в традиционый конвейер рендерига, у них нет предопределеных входных значений и предопределенных выодных значений, они никак не взаимодействуют с другими шейдерами.

При этом вычислительне шейдеры пишутся на стандартном языке для написания шейдеров (GLSL), они имеют доступ к uniform-переменным и блокам, текстурам и изображениям, атомарным ссчетчикам и SSBO (см. расширение ARB_shader_storage_buffer_object).

ДЛя выполнения вычислоительных шейдеров служат специальные команды, вводимые этим расширением - glDispatchCompute и glDispatchComputeIndirect.

void glDispatchCompute         ( GLuint numGroupsX, GLuint numGroupsY, 
                                 GLuint numGroupsZ );
void glDispatchComputeIndirect ( GLintptr offset );

Отличие команды glDispatchComputeIndirect от glDispatchCompute заключается в том, что для нее три целочисленных параметра берутся из VBO, выбранного текущим для цели GL_DISPATCH_INDIRECT_BUFFER, начиная со смещения offset.

Команды glDispatchCompute* запускают вычислительный шкйдер из текущей программы. Шейдер запускается на большом количестве нитей, каждая такая нить называется work item. Все нити группируются в рабочие группы (work group). Это практически полностью совпадает с тем, что есть в CUDA и OpenCL. При этом нити, входящие в одну рабочую группу, могут обмениваться данными через разделяемую (shared) память и использовать барьерную синхронизацию (через функцию barrier).

Каждая рабочая группа представляет собой 1,2 или 3-мерный массив из нитей, при этом как сам GPU, так и реализация OpenGL могут иметь свои ограничения на максимальный размер группы по каждому измерению и максимальное число нитей в группе.

int maxX, maxY, maxZ, maxItemsPerGroup;

glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_SIZE, 0, &maxX );
glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_SIZE, 1, &maxY );
glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_SIZE, 2, &maxZ );
glGetIntegerv   ( GL_MAX_COMPUTE_WORK_GROUP_INVOCATIONS, &maxItemsPerGroup );

При запуске нитей все блоки должны иметь одну и ту же топологию. Сами группы образуют трехмерный массив размерами этого массива будут numGroupsX, numGroupsY и numGroupsZ (т.е. параметры вызова glDispatchCompute).

Также есть ограничения на максимальное число групп по каждому измерению, их можно узнать при помощи следующего фрагмента кода:

int maxGroupX, maxGroupY, maxGroupZ;

glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_COUNT, 0, &maxGroupX );
glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_COUNT, 1, &maxGroupY );
glGetIntegeri_v ( GL_MAX_COMPUTE_WORK_GROUP_COUNT, 2, &maxGroupZ );

Обратите внимане, что группы выполняются независимо друг от друга, в то время как нити внутри одной группы выполняются параллельно.

Если количество групп задается в команде glDispatchCompute*, то количество нитей внутри группы задается непосредственно в теле самого шейдера при помощи директивы layout:

layout ( local_size_x = 32, local_size_y = 16 ), local_size_z = 2 ) in;

Если один из трех размеров не был указан, то он считается равным единице. Если заданные размеры группы не соответствуют максимальным ограничениям, то такой шейдер просто не скопмилируется.

Внутри шейдера можно через встроенные (доступные только для чтения) переменные получить информацию как о количестве запущенных групп, так и о топологии самой группы и о месте текущей нити (номер ее блока и ее номер внутри блока). Ниже приводится списко доступных вычислительному шейдеру встроеных переменных.

uvec3 gl_NumWorkGroups;
uvec3 gl_WorkGroupSize;
uvec3 gl_WorkGroupID;
uvec3 gl_LocalInvocationID;
uvec3 gl_GlobalInvocationID;
uint  gl_LocalInvocationIndex;

В переменной gl_NumWorkGroups содержится общее число запущенных групп по всем трем измерениям и оно совпадает с соответствующим парамтеров вызова glDispatchCompute*. В переменной gl_WorkGroupSize содержится размер рабочей группы по кажжому измерению и оно совпадает с параметрами local_size_* из директивы layout в теле шейдера.

В переменной gl_WorkGroupID содержится тремхерный индекс группы, к которой принадлежит текущая нить. В переменной gl_LocalInvocationID содержится тремерный индекс текущей нити внутри своей группы.

В переменной gl_GlobalInvocationID хранится глобальны трехмерный индекс текущей нити среди всех запущенных нитей. Он вычисляется как gl_WorkGroupID*gl_WorkGroupSize+gl_LocalInvocationID.

Довольно часто нужен не трезмерный индекс нити, а одномерный индекс - фактически просто номер нити среди всех запущенных нитей группы. Его можно получить из переменной gl_LocalInvocationIndex.

Кроме того каждая группа может получить в свое распоряжение некоторы объем разделяемой (shared) памяти, доступной сразу всем нитям группы. При этом разные группы поулчают неперекрывающиеся блоки разделяемой памяти. Ниже приводится пример описания массива в разделяемой памяти.

shared float buf [gl_WorkGroupSize.x+2][gl_WorkGroupSize.y+2];

Максимальный объем разделемой памяти, доступной группе можно получить следующим образом:

int maxSharedSize;

glGetIntegerv ( GL_MAX_COMPUTE_SHARED_MEMORY_SIZE, &maxSharedSize );

При работе с разделяемой памятью возникает необходимость в барьерной синхронизации (так как сразу много нитей могут одновременно изменять содержимое этой памяти). Ее обеспечивает встрлоенная в GLSL функция barrier. Вызов этой функции гарантирует, что или все нити блока еще не дошли до этого вызова или же они все уже прошли за него - если при этом вызове есть нити, не дошедшие до этого вызова, то остальные нити блока будут ожидать, пока те не дойдут до вызова barrier.

Кроме того вводится целый ряд специальных функций для синхронизации доступа к памяти.

void memoryBarrier ();
void memoryBarrierAtomicCounter ();
void memoryBarrierBuffer ();
void memoryBarrierImage ();
void memoryBarrierShared ();
void groupMemoryBarrier ();

Функции memoryBarrier* исполльзуются для упорядочения чтения и записи в соответствующую память разными нитями. При вызове такой функции происходит ожидание завершения всех команд чтения и записи в переменные соответсвтующего типа. По возвращении из этих функций данные, записанные одной нитю, будут гарантированно видны остальным нитям.

Например, если без использования этих функций, одна нить проивзодит две записи в coherent-переменные, то другая нить может увидеть результат второй записи до того,как увидит результат первой.

В качестве примера использование вычислительных шейдеров давайте рассмотрим анимацию системы частиц под действием сил (lde[ "черных дыр"). Движение каждой частицы будет полностью определяться результирующей силой притяжения от этих "черных дыр".

Параметры честиц будут храниться в двух буферах - в одном будут находиться координаты частицы (vec4), а в другом - скорости (vec4). Для анимации этих частиц мы будем использовать следующий вычислительный шейдер.

#version 430

layout( local_size_x = 1000 ) in;

uniform float gravity1        = 1000.0;
uniform vec3  blackHolePos1   = vec3(5,0,0);
uniform float gravity2        = 1000.0;
uniform vec3  blackHolePos2   = vec3(-5,0,0);
uniform float particleInvMass = 1.0 / 0.1;
uniform float deltaT          = 0.0005;
uniform float maxDist         = 45.0;

layout(std430, binding = 0) buffer Pos 
{
    vec4 position [];
};

layout(std430, binding = 1) buffer Vel 
{
    vec4 velocity [];
};

void main() 
{
    uint idx = gl_GlobalInvocationID.x;
    vec3 p   = position [idx].xyz;

        // Force from black hole #1
    vec3  d     = blackHolePos1 - p;
    float dist  = length(d);
    vec3  force = (gravity1 / dist) * normalize(d);

        // Force from black hole #2
    d      = blackHolePos2 - p;
    dist   = length(d);
    force += (gravity2 / dist) * normalize(d);

        // Reset particles that get too far from the attractors
    if ( dist > maxDist ) 
        position [idx] = vec4(0,0,0,1);
    else 
    {
            // Apply simple Euler integrator
        vec3 a = force * particleInvMass;

        position [idx] = vec4( p + velocity[idx].xyz * deltaT + 0.5 * a * deltaT * deltaT, 1.0 );
        velocity [idx] = vec4( velocity[idx].xyz + a * deltaT, 0.0 );
    }
}

Ниже приводится код на C++ для построения и рендеринга анимации частиц. Обратите внимание на использование функции glMemoryBarrier после glDispatchCompute. Это необходимо для того, чтобы гарантирвать что результаты работы вычислительного шейдера будут видны при рендеринге частиц.

#include    <vector>
#include    "GlutRotateWindow.h"
#include    "Program.h"
#include    "VertexArray.h"
#include    "VertexBuffer.h"
#include    "vec3.h"
#include    "mat4.h"

class   MyWindow : public GlutRotateWindow
{
    vec3         eye;
    Program      compute;           // compute shader for particle animation
    Program      render;            // shader to render particles
    VertexArray  vao;               // VAO to hold buffers settings for render
    VertexBuffer vel;               // velocities buffer
    VertexBuffer pos;               // position buffer
    size_t       n;
    size_t       numParticles;
    float        t;                 // current time in seconds
    mat4         proj;

public:
    MyWindow () : GlutRotateWindow ( 200, 200, 800, 600, "Compute shader particles" ), eye ( -0.5, 0.5, 10 )
    {
                // load shaders
        if ( !compute.loadProgram ( "particles-compute.glsl" ) )
        {
            printf ( "Error loading shader: %s\n", compute.getLog ().c_str () );
            
            exit ( 1 );
        }
            
        if ( !render.loadProgram ( "particles-render.glsl" ) )
        {
            printf ( "Error loading shader: %s\n", render.getLog ().c_str () );
            
            exit ( 1 );
        }
        
                    // set particles count and other vars
        t            = 0;
        n            = 100;
        numParticles = n * n * n;
        
                    // init buffers with particle data
        std::vector<vec4> vb;
        std::vector<vec4> pb;
        float             h = 2.0f / (n - 1);

        for ( int i = 0; i < n; i++ )
            for ( int j = 0; j < n; j++ )
                for ( int k = 0; k < n; k++ ) 
                {
                    vec4    p ( h * i - 1, h * j - 1, h * k - 1, 1 );
                    
                    pb.push_back ( p );
                    vb.push_back ( vec4::zero );
                }

                    // create VBO'a and VAO
                    // use bindBase to bind to binding points
        vel.create   ();
        vel.bindBase ( GL_SHADER_STORAGE_BUFFER, 1 );
        vel.setData  ( numParticles * sizeof ( vec4 ), vb.data (), GL_DYNAMIC_DRAW );
        vel.unbind   ();
        pos.create   ();
        pos.bindBase ( GL_SHADER_STORAGE_BUFFER, 0 );
        pos.setData  ( numParticles * sizeof ( vec4 ), pb.data (), GL_DYNAMIC_DRAW );
        vel.unbind   ();

        vao.create        ();
        vao.bind          ();
        pos.bind          ( GL_ARRAY_BUFFER );
        render.bind       ();
        render.setAttrPtr ( "pos", 4, sizeof ( vec4 ), (void *) 0 );
        render.unbind     ();
        vao.unbind        ();

        glEnable    ( GL_BLEND );
        glBlendFunc ( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
    }

    void redisplay ()
    {
        glClear ( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

                // animate particles, not using of memory barrier
        compute.bind      ();
        glDispatchCompute ( numParticles / 1000, 1, 1 );
        glMemoryBarrier   ( GL_SHADER_STORAGE_BARRIER_BIT );
        compute.unbind    ();

                // render particles
        vao.bind                ();
        render.bind             ();
        render.setUniformMatrix ( "mvp",  proj*getRotation () );
        glPointSize             ( 1.0f );
        glDrawArrays            ( GL_POINTS, 0, numParticles );
        render.unbind           ();
        vao.unbind              ();
    }
    
    void reshape ( int w, int h )
    {
        GlutWindow::reshape ( w, h );
        
        proj = perspective ( 60.0f, (float)w / (float)h, 0.1f, 100.0f ) * lookAt ( eye, vec3 :: zero, vec3 ( 0, 0, 1 ) );

        glViewport ( 0, 0, (GLsizei)w, (GLsizei)h );
    }

    void mouseWheel ( int wheel, int dir, int x, int y )
    {
        eye += 0.5 * vec3 ( dir, dir, dir );
        
                    // since eye value has changed
        reshape ( getWidth(), getHeight() );
        
        glutPostRedisplay ();
    }
    
    void    keyTyped ( unsigned char key, int modifiers, int x, int y )
    {
        if ( key == 27 || key == 'q' || key == 'Q' )    //  quit requested
            exit ( 0 );
    }
    
    void    idle ()
    {
        float now = getTime ();

        compute.bind            ();
        compute.setUniformFloat ( "deltaT", (now - t)/10 );
        compute.unbind          ();
        
        t = now;
        
        glutPostRedisplay ();
    }
};
    
int main ( int argc, char * argv [] )
{
    GlutWindow::init( argc, argv );
    
    MyWindow    win;
    
    GlutWindow::run ();
    
    return 0;
}

По этой ссылке можно скачать весь исходный код к этой статье. Также доступны для скачивания откомпилированные версии для M$ Windows.