CUDA C использует флоп одинарной точности на двойных

Проблема

Во время проекта в CUDA C я столкнулся с неожиданным поведением в отношении операций с плавающей запятой одинарной и двойной точности. В проекте я сначала заполняю массив числами в ядре, а в другом ядре я делаю некоторые вычисления с этими числами. Все переменные и массивы имеют двойную точность, поэтому я бы не ожидал, что произойдет какая-либо операция с плавающей запятой одинарной точности. Однако, если я анализирую исполняемый файл программы с помощью NVPROF, он показывает, что выполняются операции с одинарной точностью. Как это возможно?

Минимальный, полный и проверяемый пример

Вот самая маленькая программа, которая показывает такое поведение на моей архитектуре: (упущены утверждения и перехват ошибок). Я использую видеокарту Nvidia Tesla k40.

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

#define Nx 10
#define Ny 10
#define RANDOM double(0.236954587566)

__global__ void test(double *array, size_t pitch){
    double rho, u;
    int x = threadIdx.x + blockDim.x*blockIdx.x;
    int y = threadIdx.y + blockDim.y*blockIdx.y;
    int idx = y*(pitch/sizeof(double)) + 2*x;

    if(x < Nx && y < Ny){
        rho = array[idx]; 
        u = array[idx+1]/rho;
        array[idx] = rho*u;
    }
}

__global__ void fill(double *array, size_t pitch){
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;  
    int idx = y*(pitch/sizeof(double)) + 2*x;

    if(x < Nx || y < Ny){
        array[idx] = RANDOM*idx;
        array[idx + 1] = idx*idx*RANDOM;
    }
}

int main(int argc, char* argv[]) {
    double *d_array;
    size_t pitch;
    cudaMallocPitch((void **) &d_array, &pitch, 2*Nx*sizeof(double), Ny);

    dim3 threadDistribution = dim3(8,8);
    dim3 blockDistribution = dim3( (Nx + threadDistribution.x - 1) / (threadDistribution.x), (Ny + threadDistribution.y) / (threadDistribution.y));
    fill <<< blockDistribution, threadDistribution >>> (d_array, pitch);
    cudaDeviceSynchronize();
    test <<< blockDistribution, threadDistribution >>> (d_array, pitch);

    return 0;
}

Вывод NVPROF (отредактирован, чтобы сделать его более читабельным, если вам нужен полный вывод, просто спросите в комментариях):

....
Device "Tesla K40c (0)"
Kernel: test(double*, unsigned long)
      Metric Name             Min         Max         Avg
      flop_count_sp           198         198         198
      flop_count_sp_add         0           0           0
      flop_count_sp_mul         0           0           0
      flop_count_sp_fma        99          99          99
      flop_count_sp_special   102         102         102
      flop_count_dp          1214        1214        1214
      flop_count_dp_add         0           0           0
      flop_count_dp_mul       204         204         204
      flop_count_dp_fma       505         505         505

Что я нашел до сих пор

Я обнаружил, что если я удалю деление в строке 16:

u = array[idx+1]/rho;
==>
u = array[idx+1];

результат соответствует ожидаемому: ноль операций с одинарной точностью и ровно 100 операций с двойной точностью. Кто-нибудь знает, почему деление заставляет программу использовать флоп с одинарной точностью и в 10 раз больше операций с плавающей запятой с двойной точностью? Я также пытался использовать встроенные функции (__ddiv_rn), но это не решило проблему.

Спасибо заранее!

Изменить - рабочее решение

Хотя я до сих пор не понял, почему он использует одинарную точность, я нашел «решение» этой проблемы благодаря @EOF. Замена деления на умножение обратным значением ро сработала:

u = array[idx+1]/rho;
==>
u = array[idx+1]*__drcp_rn(rho);

precision cuda single-precision
person Frank    schedule 16.01.2017    source источник
comment
Я позволю кому-то другому ответить более авторитетно, но устройства CUDA не имеют аппаратной поддержки деления с плавающей запятой. Вместо этого деление вызовет подпрограмму, сгенерированную компилятором. Вполне вероятно, что указанная подпрограмма использует внутреннюю математику одинарной точности для реализации алгоритма деления.   -  person Jason R    schedule 16.01.2017
comment
Я ожидаю, что устройство будет использовать какую-то обратную аппроксимацию, требующую итеративного уточнения.   -  person EOF    schedule 16.01.2017
comment
@EOF: Да, и делать первые несколько итераций в SP было бы вполне разумно.   -  person Martin Bonner supports Monica    schedule 16.01.2017
comment
Поскольку ваш код написан на C и вопрос касается языка C, почему тег C++? C и C++ - разные языки. Например, в C++ у вас может быть std::vector<double>, которого нет в C.   -  person Thomas Matthews    schedule 16.01.2017
comment
@thomas Matthews: код не написан на C и скомпилирован с помощью компилятора C++. А математическая библиотека устройства реализована как шаблонные перегрузки стандартной математической библиотеки C++. Для вас это может не выглядеть как C++, но уж точно не C.   -  person talonmies    schedule 16.01.2017
comment
@talonmies: Как код C++? Включает в себя stdio.h, stdlib.h и math.h библиотеки C. В C++ системные включения не имеют расширений.   -  person Thomas Matthews    schedule 16.01.2017
comment
@thomas Matthews: код, который вы видите, проходит через синтаксический анализатор, который генерирует C++. Компилятор C++ компилирует код. Библиотеки CUDA, которые автоматически импортируются, являются библиотеками шаблонов C++. Этот код не C, даже если вы думаете, что он похож на C   -  person talonmies    schedule 16.01.2017
comment
@ThomasMatthews: я обновил теги, спасибо за уведомление. Я хотел написать код CUDA C. Я скомпилировал с помощью nvcc с указанным gcc. Math.h также доступен для C, а не только для C++.   -  person Frank    schedule 16.01.2017
comment
@ЭОФ. Я тоже об этом думал, но документации по этому поводу не нашел. У вас есть идеи, где я могу найти больше информации по этой теме? Кроме того, будет ли вывод иметь двойную точность?   -  person Frank    schedule 16.01.2017
comment
@Frank Какую версию CUDA вы использовали для компиляции этого кода? Какие ключи компилятора вы передали, в частности, какую целевую архитектуру вы указали? Используя CUDA 7.5 для создания кода, как опубликовано, с -arch=sm_35, а затем дизассемблированного с cuobjdump --dump-sass, я вижу только две инструкции MUFU.RCP64H как часть вызываемой подпрограммы деления с двойной точностью, которые, возможно, считаются flop_count_sp_special. Я не вижу никаких инструкций, которые, как мне кажется, можно классифицировать как flop_count_sp_fma.   -  person njuffa    schedule 16.01.2017
comment
@Frank: На самом деле я мало что знаю о gpgpu, но процессоры, как правило, имеют инструкции для быстрого нахождения приблизительных обратных величин (и приблизительных обратных квадратных корней) и с большей пропускной способностью, чем точные инструкции деления / квадратного корня. Примером программного обеспечения является быстрый обратный квадратный корень из славы Quake. Приближения могут быть улучшены чем-то вроде итеративного уточнения рядов Ньютона-Рафсона или Тейлора. Вы можете, вероятно, в конце концов получить double-точность, но для этого потребуется довольно много инструкций.   -  person EOF    schedule 16.01.2017
comment
@njuffa: я использовал CUDA 8.0 и архитектуру sm_35. Я не знаком с дизассемблированием и MUFU.RCP64H, но я обязательно изучу его. Спасибо!   -  person Frank    schedule 16.01.2017
comment
@EOF Спасибо за информацию! Как вы можете видеть в посте, умножение на обратное не требует операций с одинарной точностью. Я думаю, вы правы в том, как GPU вычисляет деление.   -  person Frank    schedule 16.01.2017
comment
@Frank Чтобы разобрать, просто используйте cuobjdump --dump-sass [executable name]. Инструкции одинарной точности обычно начинаются с F, например. FFMA, FMUL; не всегда, напр. I2F.F32.* — преобразование целого числа в одинарную точность. MUFU.RCPH64 — это инструкция обратной аппроксимации, реализованная в multiфункциональном блоке. Возможно, подпрограмма деления двойной точности изменилась в CUDA 8.0 (предположительно лучше оптимизирована). Ищите первую инструкцию CAL в дизассемблере, чтобы найти вызов деления; внутри этой подпрограммы может быть другой код CAL для медленного пути   -  person njuffa    schedule 16.01.2017
comment
@njuffa MUFU.RCP64H, по-видимому, считается в профилировщике как Операции с плавающей запятой (Single Precision Special). В коде есть еще одна FFMA с одинарной точностью, которая, по-видимому, используется в старшем слове итеративной коррекции как более быстрое условное выражение по сравнению с 64-битной математикой.   -  person tera    schedule 16.01.2017
comment
Использование FFMA с одним аргументом, являющимся нулевым регистром, кажется мне странным, но поскольку все эффекты на флаги/предикаты недокументированы, я подумал, что лучше предоставить вам объяснение кода.   -  person tera    schedule 16.01.2017
comment
@tera Я подозревал, что что-то подобное могло быть введено в качестве оптимизации, но на данный момент у меня нет доступа к CUDA 8.0. Я бы посоветовал кому-нибудь опубликовать ответ, показывающий соответствующую разборку процедуры деления с двойной точностью из CUDA 8.0, выделяя инструкции с плавающей запятой с одинарной точностью. Вопрос не в том, зачем нужны инструкции, а в том, откуда они берутся. Ответ: Изнутри кода деления двойной точности.   -  person njuffa    schedule 16.01.2017

Ответы (1)


arrow_upward
6
arrow_downward

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

Даже встроенный __ddiv_rn() скомпилирован в эту последовательность инструкций с помощью ptxas, поэтому его использование не имеет значения.

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

В качестве примера я буду использовать следующее базовое ядро ​​деления:

__global__ void div(double x, double y, double *z) {
    *z = x / y;
}

Это скомпилировано в следующую сборку шейдера для устройства с вычислительными возможностями 3.5:

    Function : _Z3divddPd
.headerflags    @"EF_CUDA_SM35 EF_CUDA_PTX_SM(EF_CUDA_SM35)"
                                                                                          /* 0x08a0109c10801000 */
    /*0008*/                   MOV R1, c[0x0][0x44];                                      /* 0x64c03c00089c0006 */
    /*0010*/                   MOV R0, c[0x0][0x14c];                                     /* 0x64c03c00299c0002 */
    /*0018*/                   MOV32I R2, 0x1;                                            /* 0x74000000009fc00a */
    /*0020*/                   MOV R8, c[0x0][0x148];                                     /* 0x64c03c00291c0022 */
    /*0028*/                   MOV R9, c[0x0][0x14c];                                     /* 0x64c03c00299c0026 */
    /*0030*/                   MUFU.RCP64H R3, R0;                                        /* 0x84000000031c000e */
    /*0038*/                   MOV32I R0, 0x35b7333;                                      /* 0x7401adb9999fc002 */
                                                                                          /* 0x08a080a080a4a4a4 */
    /*0048*/                   DFMA R4, -R8, R2, c[0x2][0x0];                             /* 0x9b880840001c2012 */
    /*0050*/                   DFMA R4, R4, R4, R4;                                       /* 0xdb801000021c1012 */
    /*0058*/                   DFMA R4, R4, R2, R2;                                       /* 0xdb800800011c1012 */
    /*0060*/                   DMUL R6, R4, c[0x0][0x140];                                /* 0x64000000281c101a */
    /*0068*/                   FSETP.GE.AND P0, PT, R0, |c[0x0][0x144]|, PT;              /* 0x5db09c00289c001e */
    /*0070*/                   DFMA R8, -R8, R6, c[0x0][0x140];                           /* 0x9b881800281c2022 */
    /*0078*/                   MOV R2, c[0x0][0x150];                                     /* 0x64c03c002a1c000a */
                                                                                          /* 0x0880acb0a0ac8010 */
    /*0088*/                   MOV R3, c[0x0][0x154];                                     /* 0x64c03c002a9c000e */
    /*0090*/                   DFMA R4, R8, R4, R6;                                       /* 0xdb801800021c2012 */
    /*0098*/               @P0 BRA 0xb8;                                                  /* 0x120000000c00003c */
    /*00a0*/                   FFMA R0, RZ, c[0x0][0x14c], R5;                            /* 0x4c001400299ffc02 */
    /*00a8*/                   FSETP.GT.AND P0, PT, |R0|, c[0x2][0x8], PT;                /* 0x5da01c40011c021e */
    /*00b0*/               @P0 BRA 0xe8;                                                  /* 0x120000001800003c */
    /*00b8*/                   MOV R4, c[0x0][0x140];                                     /* 0x64c03c00281c0012 */
                                                                                          /* 0x08a1b810b8008010 */
    /*00c8*/                   MOV R5, c[0x0][0x144];                                     /* 0x64c03c00289c0016 */
    /*00d0*/                   MOV R7, c[0x0][0x14c];                                     /* 0x64c03c00299c001e */
    /*00d8*/                   MOV R6, c[0x0][0x148];                                     /* 0x64c03c00291c001a */
    /*00e0*/                   CAL 0xf8;                                                  /* 0x1300000008000100 */
    /*00e8*/                   ST.E.64 [R2], R4;                                          /* 0xe5800000001c0810 */
    /*00f0*/                   EXIT;                                                      /* 0x18000000001c003c */
    /*00f8*/                   LOP32I.AND R0, R7, 0x40000000;                             /* 0x20200000001c1c00 */
                                                                                          /* 0x08a08010a010b010 */
    /*0108*/                   MOV32I R15, 0x1ff00000;                                    /* 0x740ff800001fc03e */
    /*0110*/                   ISETP.LT.U32.AND P0, PT, R0, c[0x2][0xc], PT;              /* 0x5b101c40019c001e */
    /*0118*/                   MOV R8, RZ;                                                /* 0xe4c03c007f9c0022 */
    /*0120*/                   SEL R9, R15, c[0x2][0x10], !P0;                            /* 0x65002040021c3c26 */
    /*0128*/                   MOV32I R12, 0x1;                                           /* 0x74000000009fc032 */
    /*0130*/                   DMUL R10, R8, R6;                                          /* 0xe4000000031c202a */
    /*0138*/                   LOP32I.AND R0, R5, 0x7f800000;                             /* 0x203fc000001c1400 */
                                                                                          /* 0x08a0108ca01080a0 */
    /*0148*/                   MUFU.RCP64H R13, R11;                                      /* 0x84000000031c2c36 */
    /*0150*/                   DFMA R16, -R10, R12, c[0x2][0x0];                          /* 0x9b883040001c2842 */
    /*0158*/                   ISETP.LT.U32.AND P0, PT, R0, c[0x2][0x14], PT;             /* 0x5b101c40029c001e */
    /*0160*/                   MOV R14, RZ;                                               /* 0xe4c03c007f9c003a */
    /*0168*/                   DFMA R16, R16, R16, R16;                                   /* 0xdb804000081c4042 */
    /*0170*/                   SEL R15, R15, c[0x2][0x10], !P0;                           /* 0x65002040021c3c3e */
    /*0178*/                   SSY 0x3a0;                                                 /* 0x1480000110000000 */
                                                                                          /* 0x08acb4a4a4a4a480 */
    /*0188*/                   DMUL R14, R14, R4;                                         /* 0xe4000000021c383a */
    /*0190*/                   DFMA R12, R16, R12, R12;                                   /* 0xdb803000061c4032 */
    /*0198*/                   DMUL R16, R14, R12;                                        /* 0xe4000000061c3842 */
    /*01a0*/                   DFMA R10, -R10, R16, R14;                                  /* 0xdb883800081c282a */
    /*01a8*/                   DFMA R10, R10, R12, R16;                                   /* 0xdb804000061c282a */
    /*01b0*/                   DSETP.LEU.AND P0, PT, |R10|, RZ, PT;                       /* 0xdc581c007f9c2a1e */
    /*01b8*/              @!P0 BRA 0x1e0;                                                 /* 0x120000001020003c */
                                                                                          /* 0x088010b010b8acb4 */
    /*01c8*/                   DSETP.EQ.AND P0, PT, R10, RZ, PT;                          /* 0xdc101c007f9c281e */
    /*01d0*/              @!P0 BRA 0x358;                                                 /* 0x12000000c020003c */
    /*01d8*/                   DMUL.S R8, R4, R6;                                         /* 0xe4000000035c1022 */
    /*01e0*/                   ISETP.GT.U32.AND P0, PT, R0, c[0x2][0x18], PT;             /* 0x5b401c40031c001e */
    /*01e8*/                   MOV32I R0, 0x1ff00000;                                     /* 0x740ff800001fc002 */
    /*01f0*/                   MOV R14, RZ;                                               /* 0xe4c03c007f9c003a */
    /*01f8*/                   SEL R15, R0, c[0x2][0x10], !P0;                            /* 0x65002040021c003e */
                                                                                          /* 0x08b4a49c849c849c */
    /*0208*/                   DMUL R12, R10, R8;                                         /* 0xe4000000041c2832 */
    /*0210*/                   DMUL R18, R10, R14;                                        /* 0xe4000000071c284a */
    /*0218*/                   DMUL R10, R12, R14;                                        /* 0xe4000000071c302a */
    /*0220*/                   DMUL R16, R8, R18;                                         /* 0xe4000000091c2042 */
    /*0228*/                   DFMA R8, R10, R6, -R4;                                     /* 0xdb901000031c2822 */
    /*0230*/                   DFMA R12, R16, R6, -R4;                                    /* 0xdb901000031c4032 */
    /*0238*/                   DSETP.GT.AND P0, PT, |R8|, |R12|, PT;                      /* 0xdc209c00061c221e */
                                                                                          /* 0x08b010ac10b010a0 */
    /*0248*/                   SEL R9, R17, R11, P0;                                      /* 0xe5000000059c4426 */
    /*0250*/                   FSETP.GTU.AND P1, PT, |R9|, 1.469367938527859385e-39, PT;  /* 0xb5e01c00801c263d */
    /*0258*/                   MOV R11, R9;                                               /* 0xe4c03c00049c002e */
    /*0260*/                   SEL R8, R16, R10, P0;                                      /* 0xe5000000051c4022 */
    /*0268*/               @P1 NOP.S;                                                     /* 0x8580000000443c02 */
    /*0270*/                   FSETP.LT.AND P0, PT, |R5|, 1.5046327690525280102e-36, PT;  /* 0xb5881c20001c161d */
    /*0278*/                   MOV32I R0, 0x3ff00000;                                     /* 0x741ff800001fc002 */
                                                                                          /* 0x0880a48090108c10 */
    /*0288*/                   MOV R16, RZ;                                               /* 0xe4c03c007f9c0042 */
    /*0290*/                   SEL R17, R0, c[0x2][0x1c], !P0;                            /* 0x65002040039c0046 */
    /*0298*/                   LOP.OR R10, R8, 0x1;                                       /* 0xc2001000009c2029 */
    /*02a0*/                   LOP.AND R8, R8, -0x2;                                      /* 0xca0003ffff1c2021 */
    /*02a8*/                   DMUL R4, R16, R4;                                          /* 0xe4000000021c4012 */
    /*02b0*/                   DMUL R6, R16, R6;                                          /* 0xe4000000031c401a */
    /*02b8*/                   DFMA R14, R10, R6, -R4;                                    /* 0xdb901000031c283a */
                                                                                          /* 0x08b010b010a0b4a4 */
    /*02c8*/                   DFMA R12, R8, R6, -R4;                                     /* 0xdb901000031c2032 */
    /*02d0*/                   DSETP.GT.AND P0, PT, |R12|, |R14|, PT;                     /* 0xdc209c00071c321e */
    /*02d8*/                   SEL R8, R10, R8, P0;                                       /* 0xe5000000041c2822 */
    /*02e0*/                   LOP.AND R0, R8, 0x1;                                       /* 0xc2000000009c2001 */
    /*02e8*/                   IADD R11.CC, R8, -0x1;                                     /* 0xc88403ffff9c202d */
    /*02f0*/                   ISETP.EQ.U32.AND P0, PT, R0, 0x1, PT;                      /* 0xb3201c00009c001d */
    /*02f8*/                   IADD.X R0, R9, -0x1;                                       /* 0xc88043ffff9c2401 */
                                                                                          /* 0x08b4a480a010b010 */
    /*0308*/                   SEL R10, R11, R8, !P0;                                     /* 0xe5002000041c2c2a */
    /*0310*/               @P0 IADD R8.CC, R8, 0x1;                                       /* 0xc084000000802021 */
    /*0318*/                   SEL R11, R0, R9, !P0;                                      /* 0xe5002000049c002e */
    /*0320*/               @P0 IADD.X R9, R9, RZ;                                         /* 0xe08040007f802426 */
    /*0328*/                   DFMA R14, R10, R6, -R4;                                    /* 0xdb901000031c283a */
    /*0330*/                   DFMA R4, R8, R6, -R4;                                      /* 0xdb901000031c2012 */
    /*0338*/                   DSETP.GT.AND P0, PT, |R4|, |R14|, PT;                      /* 0xdc209c00071c121e */
                                                                                          /* 0x08b4acb4a010b810 */
    /*0348*/                   SEL R8, R10, R8, P0;                                       /* 0xe5000000041c2822 */
    /*0350*/                   SEL.S R9, R11, R9, P0;                                     /* 0xe500000004dc2c26 */
    /*0358*/                   MOV R8, RZ;                                                /* 0xe4c03c007f9c0022 */
    /*0360*/                   MUFU.RCP64H R9, R7;                                        /* 0x84000000031c1c26 */
    /*0368*/                   DSETP.GT.AND P0, PT, |R8|, RZ, PT;                         /* 0xdc201c007f9c221e */
    /*0370*/               @P0 BRA.U 0x398;                                               /* 0x120000001000023c */
    /*0378*/              @!P0 DSETP.NEU.AND P1, PT, |R6|, +INF , PT;                     /* 0xb4681fff80201a3d */
                                                                                          /* 0x0800b8a010ac0010 */
    /*0388*/              @!P0 SEL R9, R7, R9, P1;                                        /* 0xe500040004a01c26 */
    /*0390*/              @!P0 SEL R8, R6, RZ, P1;                                        /* 0xe50004007fa01822 */
    /*0398*/                   DMUL.S R8, R8, R4;                                         /* 0xe4000000025c2022 */
    /*03a0*/                   MOV R4, R8;                                                /* 0xe4c03c00041c0012 */
    /*03a8*/                   MOV R5, R9;                                                /* 0xe4c03c00049c0016 */
    /*03b0*/                   RET;                                                       /* 0x19000000001c003c */
    /*03b8*/                   BRA 0x3b8;                                                 /* 0x12007ffffc1c003c */

Инструкция MUFU.RCP64H обеспечивает начальное приближение обратной величины. Он работает со старшими 32 битами знаменателя (y) и обеспечивает старшие 32 бита аппроксимации двойной точности, поэтому профилировщик считает их операциями с плавающей запятой (Single Precision Special). Ниже находится еще одна инструкция одинарной точности FFMA, по-видимому, используемая как высокопроизводительная версия проверки условного оператора, где не требуется полная точность.

person tera    schedule 16.01.2017
comment
Предположительно, сравнения с одинарной точностью FSETP также учитываются для подсчета профилировщиком инструкций с одинарной точностью. - person njuffa; 16.01.2017
comment
Да, профилировщик считает 400 инструкций с плавающей запятой одинарной точности, но только 198 операций с плавающей запятой одинарной точности. - person tera; 16.01.2017