claunia/cuetools.net
1
0
Fork 0
mirror of https://github.com/claunia/cuetools.net.git synced 2025-12-16 10:04:24 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
cuetools.net/CUETools.Codecs.FLACCL/flac.cl
Wolfgang Stöggl d4dd402961 Bump copyright year to 2021 
- Substitute occurrences of "-2020" with "-2021" using:
  git grep -I -l -e '-2020' -- ':(exclude)*.bak' | xargs \
  sed -b -i -e 's/-2020/-2021/g'
2021-01-14 02:18:32 +01:00

2293 lines
76 KiB
Common Lisp
Raw Permalink Blame History

/**
* CUETools.FLACCL: FLAC audio encoder using OpenCL
* Copyright (c) 2010-2021 Gregory S. Chudov
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#ifndef _FLACCL_KERNEL_H_
#define _FLACCL_KERNEL_H_
#if defined(__WinterPark__) || defined(__BeaverCreek__) || defined(__Turks__) || defined(__Caicos__) || defined(__Tahiti__) || defined(__Pitcairn__) || defined(__Capeverde__)
#define AMD
#elif defined(__Cayman__) || defined(__Barts__) || defined(__Cypress__) || defined(__Juniper__) || defined(__Redwood__) || defined(__Cedar__)
#define AMD
#elif defined(__ATI_RV770__) || defined(__ATI_RV730__) || defined(__ATI_RV710__)
#define AMD
#endif
#define VENDOR_ID_INTEL 0x8086
#define VENDOR_ID_NVIDIA 0x10DE
#define VENDOR_ID_ATIAMD 0x1002
#ifndef FLACCL_CPU
#if VENDOR_ID == VENDOR_ID_INTEL
#define WARP_SIZE 16
#else
#define WARP_SIZE 32
#endif
#endif
#if defined(HAVE_cl_khr_fp64) || defined(HAVE_cl_amd_fp64)
#define HAVE_DOUBLE
#define ZEROD 0.0
//#define FAST_DOUBLE
#else
#define double float
#define double4 float4
#define ZEROD 0.0f
#endif
#if defined(HAVE_DOUBLE) && defined(FAST_DOUBLE)
#define fastdouble double
#define fastdouble4 double4
#define ZEROFD 0.0
#else
#define fastdouble float
#define fastdouble4 float4
#define ZEROFD 0.0f
#endif
#if BITS_PER_SAMPLE > 16
#define MAX_RICE_PARAM 30
#define RICE_PARAM_BITS 5
#else
#define MAX_RICE_PARAM 14
#define RICE_PARAM_BITS 4
#endif
typedef enum
{
Constant = 0,
Verbatim = 1,
Fixed = 8,
LPC = 32
} SubframeType;
typedef struct
{
int residualOrder; // <= 32
int samplesOffs;
int shift;
int cbits;
int size;
int type;
int obits;
int blocksize;
int coding_method;
int channel;
int residualOffs;
int wbits;
int abits;
int porder;
int headerLen;
int encodingOffset;
} FLACCLSubframeData;
typedef struct
{
FLACCLSubframeData data;
int coefs[32]; // fixme: should be short?
} FLACCLSubframeTask;
#if 0
__kernel void clWindowRectangle(__global float* window, int windowOffset)
{
window[get_global_id(0)] = 1.0f;
}
__kernel void clWindowFlattop(__global float* window, int windowOffset)
{
float p = M_PI_F * get_global_id(0) / (get_global_size(0) - 1);
window[get_global_id(0)] = 1.0f
- 1.93f * cos(2 * p)
+ 1.29f * cos(4 * p)
- 0.388f * cos(6 * p)
+ 0.0322f * cos(8 * p);
}
__kernel void clWindowTukey(__global float* window, int windowOffset, float p)
{
int tid = get_global_id(0);
int Np = (int)(p / 2.0f * get_global_size(0)) - 1;
int Np2 = tid - (get_global_size(0) - Np - 1) + Np;
int n = select(max(Np, Np2), tid, tid <= Np);
window[tid] = 0.5f - 0.5f * cos(M_PI_F * n / Np);
}
#endif
#if BITS_PER_SAMPLE > 16
__kernel void clStereoDecorr(
__global int *samples,
__global unsigned char *src,
int offset
)
{
int pos = get_global_id(0);
int bpos = pos * 6;
int x = (((int)src[bpos] << 8) | ((int)src[bpos+1] << 16) | ((int)src[bpos+2] << 24)) >> 8;
int y = (((int)src[bpos+3] << 8) | ((int)src[bpos+4] << 16) | ((int)src[bpos+5] << 24)) >> 8;
samples[pos] = x;
samples[1 * offset + pos] = y;
samples[2 * offset + pos] = (x + y) >> 1;
samples[3 * offset + pos] = x - y;
}
__kernel void clChannelDecorr2(
__global int *samples,
__global unsigned char *src,
int offset
)
{
int pos = get_global_id(0);
int bpos = pos * 6;
samples[pos] = (((int)src[bpos] << 8) | ((int)src[bpos+1] << 16) | ((int)src[bpos+2] << 24)) >> 8;
samples[offset + pos] = (((int)src[bpos+3] << 8) | ((int)src[bpos+4] << 16) | ((int)src[bpos+5] << 24)) >> 8;
}
__kernel void clChannelDecorrX(
__global int *samples,
__global unsigned char *src,
int offset
)
{
int pos = get_global_id(0);
for (int ch = 0; ch < MAX_CHANNELS; ch++)
{
int bpos = 3 * (pos * MAX_CHANNELS + ch);
samples[offset * ch + pos] = (((int)src[bpos] << 8) | ((int)src[bpos+1] << 16) | ((int)src[bpos+2] << 24)) >> 8;
}
}
#else
__kernel void clStereoDecorr(
__global int4 *samples,
__global int4 *src,
int offset
)
{
int pos = get_global_id(0);
int4 s = src[pos];
int4 x = (s << 16) >> 16;
int4 y = s >> 16;
samples[pos] = x;
samples[1 * offset + pos] = y;
samples[2 * offset + pos] = (x + y) >> 1;
samples[3 * offset + pos] = x - y;
}
__kernel void clChannelDecorr2(
__global int4 *samples,
__global int4 *src,
int offset
)
{
int pos = get_global_id(0);
int4 s = src[pos];
samples[pos] = (s << 16) >> 16;
samples[offset + pos] = s >> 16;
}
__kernel void clChannelDecorrX(
__global int *samples,
__global short *src,
int offset
)
{
int pos = get_global_id(0);
for (int ch = 0; ch < MAX_CHANNELS; ch++)
{
int bpos = pos * MAX_CHANNELS + ch;
samples[offset * ch + pos] = src[bpos];
}
}
#endif
//__kernel void clChannelDecorr(
// int *samples,
// short *src,
// int offset
//)
//{
// int pos = get_global_id(0);
// if (pos < offset)
// samples[get_group_id(1) * offset + pos] = src[pos * get_num_groups(1) + get_group_id(1)];
//}
#define __ffs(a) (32 - clz(a & (-a)))
//#define __ffs(a) (33 - clz(~a & (a - 1)))
#ifdef FLACCL_CPU
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clFindWastedBits(
__global FLACCLSubframeTask *tasks,
__global int *samples,
int tasksPerChannel
)
{
__global FLACCLSubframeTask* ptask = &tasks[get_group_id(0) * tasksPerChannel];
int w = 0, a = 0;
for (int pos = 0; pos < ptask->data.blocksize; pos ++)
{
int smp = samples[ptask->data.samplesOffs + pos];
w |= smp;
a |= smp ^ (smp >> 31);
}
w = max(0,__ffs(w) - 1);
a = 32 - clz(a) - w;
for (int i = 0; i < tasksPerChannel; i++)
{
ptask[i].data.wbits = w;
ptask[i].data.abits = a;
//ptask[i].data.size = ptask[i].data.obits * ptask[i].data.blocksize;
}
}
#else
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clFindWastedBits(
__global FLACCLSubframeTask *tasks,
__global int *samples,
int tasksPerChannel
)
{
__local int abits[GROUP_SIZE];
__local int wbits[GROUP_SIZE];
__local FLACCLSubframeData task;
int tid = get_local_id(0);
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[get_group_id(0) * tasksPerChannel].data))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
int w = 0, a = 0;
for (int pos = tid; pos < task.blocksize; pos += GROUP_SIZE)
{
int smp = samples[task.samplesOffs + pos];
w |= smp;
a |= smp ^ (smp >> 31);
}
wbits[tid] = w;
abits[tid] = a;
barrier(CLK_LOCAL_MEM_FENCE);
for (int s = GROUP_SIZE / 2; s > 0; s >>= 1)
{
if (tid < s)
{
wbits[tid] |= wbits[tid + s];
abits[tid] |= abits[tid + s];
}
barrier(CLK_LOCAL_MEM_FENCE);
}
w = max(0,__ffs(wbits[0]) - 1);
a = 32 - clz(abits[0]) - w;
if (tid < tasksPerChannel)
{
int i = get_group_id(0) * tasksPerChannel + tid;
tasks[i].data.wbits = w;
tasks[i].data.abits = a;
//tasks[i].data.size = tasks[i].data.obits * tasks[i].data.blocksize;
}
}
#endif
#ifdef FLACCL_CPU
#define TEMPBLOCK 512
#define STORE_AC(ro, val) if (ro <= MAX_ORDER) pout[ro] = val;
#define STORE_AC4(ro, val) STORE_AC(ro*4+0, val##ro.x) STORE_AC(ro*4+1, val##ro.y) STORE_AC(ro*4+2, val##ro.z) STORE_AC(ro*4+3, val##ro.w)
// get_num_groups(0) == number of tasks
// get_num_groups(1) == number of windows
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clComputeAutocor(
__global float *output,
__global const int *samples,
__global const float *window,
__global FLACCLSubframeTask *tasks,
const int taskCount // tasks per block
)
{
FLACCLSubframeData task = tasks[get_group_id(0) * taskCount].data;
int len = task.blocksize;
int windowOffs = get_group_id(1) * len;
float data[TEMPBLOCK + MAX_ORDER + 3];
double4 ac0 = 0.0, ac1 = 0.0, ac2 = 0.0, ac3 = 0.0, ac4 = 0.0, ac5 = 0.0, ac6 = 0.0, ac7 = 0.0, ac8 = 0.0;
for (int pos = 0; pos < len; pos += TEMPBLOCK)
{
for (int tid = 0; tid < TEMPBLOCK + MAX_ORDER + 3; tid++)
data[tid] = tid < len - pos ? samples[task.samplesOffs + pos + tid] * window[windowOffs + pos + tid] : 0.0f;
for (int j = 0; j < TEMPBLOCK;)
{
float4 temp0 = 0.0f, temp1 = 0.0f, temp2 = 0.0f, temp3 = 0.0f, temp4 = 0.0f, temp5 = 0.0f, temp6 = 0.0f, temp7 = 0.0f, temp8 = 0.0f;
for (int k = 0; k < 32; k++)
{
float d0 = data[j];
temp0 += d0 * vload4(0, &data[j]);
temp1 += d0 * vload4(1, &data[j]);
#if MAX_ORDER >= 8
temp2 += d0 * vload4(2, &data[j]);
#if MAX_ORDER >= 12
temp3 += d0 * vload4(3, &data[j]);
#if MAX_ORDER >= 16
temp4 += d0 * vload4(4, &data[j]);
temp5 += d0 * vload4(5, &data[j]);
temp6 += d0 * vload4(6, &data[j]);
temp7 += d0 * vload4(7, &data[j]);
temp8 += d0 * vload4(8, &data[j]);
#endif
#endif
#endif
j++;
}
ac0 += convert_double4(temp0);
ac1 += convert_double4(temp1);
#if MAX_ORDER >= 8
ac2 += convert_double4(temp2);
#if MAX_ORDER >= 12
ac3 += convert_double4(temp3);
#if MAX_ORDER >= 16
ac4 += convert_double4(temp4);
ac5 += convert_double4(temp5);
ac6 += convert_double4(temp6);
ac7 += convert_double4(temp7);
ac8 += convert_double4(temp8);
#endif
#endif
#endif
}
}
__global float * pout = &output[(get_group_id(0) * get_num_groups(1) + get_group_id(1)) * (MAX_ORDER + 1)];
STORE_AC4(0, ac) STORE_AC4(1, ac) STORE_AC4(2, ac) STORE_AC4(3, ac)
STORE_AC4(4, ac) STORE_AC4(5, ac) STORE_AC4(6, ac) STORE_AC4(7, ac)
STORE_AC4(8, ac)
}
#else
// get_num_groups(0) == number of tasks
// get_num_groups(1) == number of windows
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clComputeAutocor(
__global float *output,
__global const int *samples,
__global const float *window,
__global FLACCLSubframeTask *tasks,
const int taskCount // tasks per block
)
{
__local fastdouble data[GROUP_SIZE * 2];
__local FLACCLSubframeData task;
const int tid = get_local_id(0);
// fetch task data
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(tasks + taskCount * get_group_id(0)))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
int bs = task.blocksize;
data[tid] = ZEROFD;
const int THREADS_FOR_ORDERS = MAX_ORDER < 8 ? 8 : MAX_ORDER < 16 ? 16 : MAX_ORDER < 32 ? 32 : 64;
int lag = tid & (THREADS_FOR_ORDERS - 1);
int tid1 = tid + GROUP_SIZE - lag;
int pos = 0;
const __global float * wptr = &window[get_group_id(1) * bs];
// const __global int * sptr = &samples[task.samplesOffs];
double corr = ZEROD;
for (pos = 0; pos + GROUP_SIZE - 1 < bs; pos += GROUP_SIZE)
{
int off = pos + tid;
// fetch samples
fastdouble nextData = samples[task.samplesOffs + off] * wptr[off];
data[tid + GROUP_SIZE] = nextData;
barrier(CLK_LOCAL_MEM_FENCE);
fastdouble4 tmp = ZEROFD;
for (int i = 0; i < THREADS_FOR_ORDERS / 4; i++)
tmp += vload4(i, &data[tid1 - lag]) * vload4(i, &data[tid1]);
corr += (tmp.x + tmp.y) + (tmp.w + tmp.z);
barrier(CLK_LOCAL_MEM_FENCE);
data[tid] = nextData;
}
{
int off = pos + tid;
// fetch samples
double nextData = off < bs ? samples[task.samplesOffs + off] * wptr[off] : ZEROD;
data[tid + GROUP_SIZE] = nextData;
barrier(CLK_LOCAL_MEM_FENCE);
fastdouble4 tmp = ZEROFD;
for (int i = 0; i < THREADS_FOR_ORDERS / 4; i++)
tmp += vload4(i, &data[tid1 - lag]) * vload4(i, &data[tid1]);
corr += (tmp.x + tmp.y) + (tmp.w + tmp.z);
}
data[tid] = corr;
barrier(CLK_LOCAL_MEM_FENCE);
for (int i = GROUP_SIZE / 2; i >= THREADS_FOR_ORDERS; i >>= 1)
{
if (tid < i)
data[tid] += data[tid + i];
barrier(CLK_LOCAL_MEM_FENCE);
}
if (tid <= MAX_ORDER)
output[(get_group_id(0) * get_num_groups(1) + get_group_id(1)) * (MAX_ORDER + 1) + tid] = data[tid];
}
#endif
#ifdef FLACCL_CPU
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clComputeLPC(
__global float *pautoc,
__global float *lpcs,
int windowCount
)
{
int lpcOffs = (get_group_id(0) + get_group_id(1) * windowCount) * (MAX_ORDER + 1) * 32;
int autocOffs = (get_group_id(0) + get_group_id(1) * get_num_groups(0)) * (MAX_ORDER + 1);
volatile double ldr[32];
volatile double gen0[32];
volatile double gen1[32];
volatile double err[32];
__global float* autoc = pautoc + autocOffs;
for (int i = 0; i < MAX_ORDER; i++)
{
gen0[i] = gen1[i] = autoc[i + 1];
ldr[i] = 0.0;
}
// Compute LPC using Schur and Levinson-Durbin recursion
double error = autoc[0];
for (int order = 0; order < MAX_ORDER; order++)
{
// Schur recursion
double reff = -gen1[0] / error;
//error += gen1[0] * reff; // Equivalent to error *= (1 - reff * reff);
error *= (1 - reff * reff);
for (int j = 0; j < MAX_ORDER - 1 - order; j++)
{
gen1[j] = gen1[j + 1] + reff * gen0[j];
gen0[j] = gen1[j + 1] * reff + gen0[j];
}
err[order] = error;
// Levinson-Durbin recursion
ldr[order] = reff;
for (int j = 0; j < order / 2; j++)
{
double tmp = ldr[j];
ldr[j] += reff * ldr[order - 1 - j];
ldr[order - 1 - j] += reff * tmp;
}
if (0 != (order & 1))
ldr[order / 2] += ldr[order / 2] * reff;
// Output coeffs
for (int j = 0; j <= order; j++)
lpcs[lpcOffs + order * 32 + j] = -ldr[order - j];
}
// Output prediction error estimates
for (int j = 0; j < MAX_ORDER; j++)
lpcs[lpcOffs + MAX_ORDER * 32 + j] = err[j];
}
#else
__kernel __attribute__((reqd_work_group_size(32, 1, 1)))
void clComputeLPC(
__global float *autoc,
__global float *lpcs,
int windowCount
)
{
__local struct {
volatile double ldr[32];
volatile double gen1[32];
volatile float error[32];
volatile float autoc[33];
} shared;
const int tid = get_local_id(0);// + get_local_id(1) * 32;
int autocOffs = (get_group_id(0) + get_group_id(1) * get_num_groups(0)) * (MAX_ORDER + 1);
int lpcOffs = autocOffs * 32;
shared.autoc[get_local_id(0)] = get_local_id(0) <= MAX_ORDER ? autoc[autocOffs + get_local_id(0)] : 0;
if (get_local_id(0) + get_local_size(0) <= MAX_ORDER)
shared.autoc[get_local_id(0) + get_local_size(0)] = autoc[autocOffs + get_local_id(0) + get_local_size(0)];
barrier(CLK_LOCAL_MEM_FENCE);
// Compute LPC using Schur and Levinson-Durbin recursion
double gen0 = shared.gen1[get_local_id(0)] = shared.autoc[get_local_id(0)+1];
shared.ldr[get_local_id(0)] = ZEROD;
double error = shared.autoc[0];
#ifdef DEBUGPRINT1
int magic = autocOffs == 0; // shared.autoc[0] == 177286873088.0f;
if (magic && get_local_id(0) <= MAX_ORDER)
printf("autoc[%d] == %f\n", get_local_id(0), shared.autoc[get_local_id(0)]);
#endif
barrier(CLK_LOCAL_MEM_FENCE);
for (int order = 0; order < MAX_ORDER; order++)
{
// Schur recursion
double reff = -shared.gen1[0] / error;
//error += shared.gen1[0] * reff; // Equivalent to error *= (1 - reff * reff);
error *= (1 - reff * reff);
double gen1;
if (get_local_id(0) < MAX_ORDER - 1 - order)
{
gen1 = shared.gen1[get_local_id(0) + 1] + reff * gen0;
gen0 += shared.gen1[get_local_id(0) + 1] * reff;
}
barrier(CLK_LOCAL_MEM_FENCE);
if (get_local_id(0) < MAX_ORDER - 1 - order)
shared.gen1[get_local_id(0)] = gen1;
#ifdef DEBUGPRINT1
if (magic && get_local_id(0) == 0)
printf("order == %d, reff == %f, error = %f\n", order, reff, error);
if (magic && get_local_id(0) <= MAX_ORDER)
printf("gen[%d] == %f, %f\n", get_local_id(0), gen0, gen1);
#endif
// Store prediction error
if (get_local_id(0) == 0)
shared.error[order] = error;
// Levinson-Durbin recursion
double ldr = shared.ldr[get_local_id(0)];
barrier(CLK_LOCAL_MEM_FENCE);
if (get_local_id(0) < order)
shared.ldr[order - 1 - get_local_id(0)] += reff * ldr;
if (get_local_id(0) == order)
shared.ldr[get_local_id(0)] += reff;
barrier(CLK_LOCAL_MEM_FENCE);
// Output coeffs
if (get_local_id(0) <= order)
lpcs[lpcOffs + order * 32 + get_local_id(0)] = -shared.ldr[order - get_local_id(0)];
//if (get_local_id(0) <= order + 1 && fabs(-shared.ldr[0]) > 3000)
// printf("coef[%d] == %f, autoc == %f, error == %f\n", get_local_id(0), -shared.ldr[order - get_local_id(0)], shared.autoc[get_local_id(0)], shared.error[get_local_id(0)]);
}
barrier(CLK_LOCAL_MEM_FENCE);
#ifdef DEBUGPRINT1
if (magic && get_local_id(0) < MAX_ORDER)
printf("error[%d] == %f\n", get_local_id(0), shared.error[get_local_id(0)]);
#endif
// Output prediction error estimates
if (get_local_id(0) < MAX_ORDER)
lpcs[lpcOffs + MAX_ORDER * 32 + get_local_id(0)] = shared.error[get_local_id(0)];
}
#endif
#ifdef FLACCL_CPU
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clQuantizeLPC(
__global FLACCLSubframeTask *tasks,
__global float*lpcs,
int taskCount, // tasks per block
int taskCountLPC, // tasks per set of coeffs (<= 32)
int minprecision,
int precisions
)
{
int bs = tasks[get_group_id(1) * taskCount].data.blocksize;
int abits = tasks[get_group_id(1) * taskCount].data.abits;
int obits = tasks[get_group_id(1) * taskCount].data.obits;
int lpcOffs = (get_group_id(0) + get_group_id(1) * get_num_groups(0)) * (MAX_ORDER + 1) * 32;
float error[MAX_ORDER];
int best_orders[MAX_ORDER];
int best8 = 0;
// Load prediction error estimates based on Akaike's Criteria
for (int tid = 0; tid < MAX_ORDER; tid++)
{
error[tid] = bs * log(1.0f + lpcs[lpcOffs + MAX_ORDER * 32 + tid]) + tid * 4.12f * log((float)bs);
best_orders[tid] = tid;
if (tid < 8 && error[tid] < error[best8])
best8 = tid;
}
#if 0
for (int i = best8 + 1; i < MAX_ORDER; i++)
error[i] += 20.5f * log((float)bs);
#endif
// Select best orders
for (int i = 0; i < MAX_ORDER && i < taskCountLPC; i++)
{
for (int j = i + 1; j < MAX_ORDER; j++)
{
if (error[best_orders[j]] < error[best_orders[i]])
{
int tmp = best_orders[j];
best_orders[j] = best_orders[i];
best_orders[i] = tmp;
}
}
}
// Quantization
for (int i = 0; i < taskCountLPC; i ++)
{
int order = best_orders[i >> precisions];
int tmpi = 0;
for (int tid = 0; tid <= order; tid ++)
{
float lpc = lpcs[lpcOffs + order * 32 + tid];
// get 15 bits of each coeff
int c = convert_int_rte(lpc * (1 << 15));
// remove sign bits
tmpi |= c ^ (c >> 31);
}
// choose precision
//int cbits = max(3, min(10, 5 + (abits >> 1))); // - convert_int_rte(shared.PE[order - 1])
#if BITS_PER_SAMPLE > 16
int cbits = max(3, min(15 - minprecision + (i - ((i >> precisions) << precisions)) - (bs <= 2304) - (bs <= 1152) - (bs <= 576), abits));
#else
int cbits = max(3, min(min(13 - minprecision + (i - ((i >> precisions) << precisions)) - (bs <= 2304) - (bs <= 1152) - (bs <= 576), abits), clz(order) + 1 - obits));
#endif
// calculate shift based on precision and number of leading zeroes in coeffs
int shift = max(0,min(15, clz(tmpi) - 18 + cbits));
int taskNo = get_group_id(1) * taskCount + get_group_id(0) * taskCountLPC + i;
tmpi = 0;
for (int tid = 0; tid <= order; tid ++)
{
float lpc = lpcs[lpcOffs + order * 32 + tid];
// quantize coeffs with given shift
int c = convert_int_rte(clamp(lpc * (1 << shift), (float)((-1 << (cbits - 1)) + 1), (float)((1 << (cbits - 1)) - 1)));
// remove sign bits
tmpi |= c ^ (c >> 31);
tasks[taskNo].coefs[tid] = c;
}
// calculate actual number of bits (+1 for sign)
cbits = 1 + 32 - clz(tmpi);
// output shift, cbits, ro
tasks[taskNo].data.shift = shift;
tasks[taskNo].data.cbits = cbits;
tasks[taskNo].data.residualOrder = order + 1;
}
}
#else
__kernel __attribute__((reqd_work_group_size(32, 1, 1)))
void clQuantizeLPC(
__global FLACCLSubframeTask *tasks,
__global float*lpcs,
int taskCount, // tasks per block
int taskCountLPC, // tasks per set of coeffs (<= 32)
int minprecision,
int precisions
)
{
__local struct {
FLACCLSubframeData task;
volatile int index[64];
volatile float error[64];
volatile int maxcoef[32];
// volatile int best8;
} shared;
const int tid = get_local_id(0);
// fetch task data
if (tid < sizeof(shared.task) / sizeof(int))
((__local int*)&shared.task)[tid] = ((__global int*)(tasks + get_group_id(1) * taskCount))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
const int lpcOffs = (get_group_id(0) + get_group_id(1) * get_num_groups(0)) * (MAX_ORDER + 1) * 32;
// Select best orders based on Akaike's Criteria
shared.index[tid] = min(MAX_ORDER - 1, tid);
shared.error[tid] = shared.task.blocksize * 64 + tid;
shared.index[32 + tid] = MAX_ORDER - 1;
shared.error[32 + tid] = shared.task.blocksize * 64 + tid + 32;
shared.maxcoef[tid] = 0;
// Load prediction error estimates
if (tid < MAX_ORDER)
shared.error[tid] = shared.task.blocksize * log(lpcs[lpcOffs + MAX_ORDER * 32 + tid]) + tid * 4.12f * log((float)shared.task.blocksize);
//shared.error[get_local_id(0)] = shared.task.blocksize * log(lpcs[lpcOffs + MAX_ORDER * 32 + get_local_id(0)]) + get_local_id(0) * 0.30f * (shared.task.abits + 1) * log(shared.task.blocksize);
#if 0
if (tid == 0)
{
int b8 = 0;
for (int i = 1; i < 8; i++)
if (shared.error[i] < shared.error[b8])
b8 = i;
shared.best8 = b8;
}
shared.error[tid] += select(0.0f, 20.5f * log((float)shared.task.blocksize), tid > shared.best8);
#endif
barrier(CLK_LOCAL_MEM_FENCE);
// Sort using bitonic sort
for(int size = 2; size < 64; size <<= 1){
//Bitonic merge
int ddd = (tid & (size / 2)) == 0;
for(int stride = size / 2; stride > 0; stride >>= 1){
int pos = 2 * tid - (tid & (stride - 1));
float e0 = shared.error[pos];
float e1 = shared.error[pos + stride];
int i0 = shared.index[pos];
int i1 = shared.index[pos + stride];
barrier(CLK_LOCAL_MEM_FENCE);
if ((e0 >= e1) == ddd)
{
shared.error[pos] = e1;
shared.error[pos + stride] = e0;
shared.index[pos] = i1;
shared.index[pos + stride] = i0;
}
barrier(CLK_LOCAL_MEM_FENCE);
}
}
//ddd == dir for the last bitonic merge step
{
for(int stride = 32; stride > 0; stride >>= 1){
//barrier(CLK_LOCAL_MEM_FENCE);
int pos = 2 * tid - (tid & (stride - 1));
float e0 = shared.error[pos];
float e1 = shared.error[pos + stride];
int i0 = shared.index[pos];
int i1 = shared.index[pos + stride];
barrier(CLK_LOCAL_MEM_FENCE);
if (e0 >= e1)
{
shared.error[pos] = e1;
shared.error[pos + stride] = e0;
shared.index[pos] = i1;
shared.index[pos + stride] = i0;
}
barrier(CLK_LOCAL_MEM_FENCE);
}
}
//shared.index[tid] = MAX_ORDER - 1;
//barrier(CLK_LOCAL_MEM_FENCE);
// Quantization
for (int i = 0; i < taskCountLPC; i ++)
{
int order = shared.index[i >> precisions];
float lpc = tid <= order ? lpcs[lpcOffs + order * 32 + tid] : 0.0f;
// get 15 bits of each coeff
int coef = convert_int_rte(lpc * (1 << 15));
// remove sign bits
atomic_or(shared.maxcoef + i, coef ^ (coef >> 31));
barrier(CLK_LOCAL_MEM_FENCE);
int cbits = min(51 - 2 * clz(shared.task.blocksize), shared.task.abits) - minprecision + (i - ((i >> precisions) << precisions));
#if BITS_PER_SAMPLE <= 16
// Limit cbits so that 32-bit arithmetics will be enough when calculating residual
// (1 << (obits - 1)) * ((1 << (cbits - 1)) - 1) * (order + 1) < (1 << 31)
// (1 << (cbits - 1)) - 1 < (1 << (32 - obits)) / (order + 1)
// (1 << (cbits - 1)) <= (1 << (32 - obits)) / (order + 1)
// (1 << (cbits - 1)) <= (1 << (32 - obits - (32 - clz(order))) <= (1 << (32 - obits)) / (order + 1)
// (1 << (cbits - 1)) <= (1 << (clz(order) - obits))
// cbits - 1 <= clz(order) - obits
// cbits <= clz(order) - obits + 1
cbits = min(cbits, clz(order) + 1 - shared.task.obits);
#endif
cbits = clamp(cbits, 3, 15);
// Calculate shift based on precision and number of leading zeroes in coeffs.
// We know that if shifted by 15, coefs require
// 33 - clz(shared.maxcoef[i]) bits;
// So to get the desired cbits, we need to shift coefs by
// 15 + cbits - (33 - clz(shared.maxcoef[i]));
int shift = clamp(clz(shared.maxcoef[i]) - 18 + cbits, 0, 15);
int lim = (1 << (cbits - 1)) - 1;
coef = clamp(convert_int_rte(lpc * (1 << shift)), -lim, lim);
// output shift, cbits and output coeffs
int taskNo = get_group_id(1) * taskCount + get_group_id(0) * taskCountLPC + i;
if (tid == 0)
tasks[taskNo].data.shift = shift;
if (tid == 0)
tasks[taskNo].data.cbits = cbits;
if (tid == 0)
tasks[taskNo].data.residualOrder = order + 1;
if (tid <= order)
tasks[taskNo].coefs[tid] = coef;
}
}
#endif
#ifdef FLACCL_CPU
#define TEMPBLOCK1 TEMPBLOCK
__kernel __attribute__(( vec_type_hint (int4))) __attribute__((reqd_work_group_size(1, 1, 1)))
void clEstimateResidual(
__global int*samples,
__global int*selectedTasks,
__global FLACCLSubframeTask *tasks
)
{
int selectedTask = selectedTasks[get_group_id(0)];
FLACCLSubframeTask task = tasks[selectedTask];
int ro = task.data.residualOrder;
int bs = task.data.blocksize;
#define ERPARTS (MAX_BLOCKSIZE >> 6)
float len[ERPARTS]; // blocksize / 64!!!!
__global int *data = &samples[task.data.samplesOffs];
for (int i = 0; i < ERPARTS; i++)
len[i] = 0.0f;
if (ro <= 4)
{
float fcoef[4];
for (int tid = 0; tid < 4; tid++)
fcoef[tid] = tid + ro - 4 < 0 ? 0.0f : - ((float) task.coefs[tid + ro - 4]) / (1 << task.data.shift);
float4 fc0 = vload4(0, &fcoef[0]);
float fdata[4];
for (int pos = 0; pos < 4; pos++)
fdata[pos] = pos + ro - 4 < 0 ? 0.0f : (float)(data[pos + ro - 4] >> task.data.wbits);
float4 fd0 = vload4(0, &fdata[0]);
for (int pos = ro; pos < bs; pos ++)
{
float4 sum4 = fc0 * fd0;
float2 sum2 = sum4.s01 + sum4.s23;
fd0 = fd0.s1230;
fd0.s3 = (float)(data[pos] >> task.data.wbits);
len[pos >> 6] += fabs(fd0.s3 + (sum2.x + sum2.y));
}
}
else if (ro <= 8)
{
float fcoef[8];
for (int tid = 0; tid < 8; tid++)
fcoef[tid] = tid + ro - 8 < 0 ? 0.0f : - ((float) task.coefs[tid + ro - 8]) / (1 << task.data.shift);
float8 fc0 = vload8(0, &fcoef[0]);
float fdata[8];
for (int pos = 0; pos < 8; pos++)
fdata[pos] = pos + ro - 8 < 0 ? 0.0f : (float)(data[pos + ro - 8] >> task.data.wbits);
float8 fd0 = vload8(0, &fdata[0]);
for (int pos = ro; pos < bs; pos ++)
{
float8 sum8 = fc0 * fd0;
float4 sum4 = sum8.s0123 + sum8.s4567;
float2 sum2 = sum4.s01 + sum4.s23;
fd0 = fd0.s12345670;
fd0.s7 = (float)(data[pos] >> task.data.wbits);
len[pos >> 6] += fabs(fd0.s7 + (sum2.x + sum2.y));
}
}
else if (ro <= 12)
{
float fcoef[12];
for (int tid = 0; tid < 12; tid++)
fcoef[tid] = tid + ro - 12 >= 0 ? - ((float) task.coefs[tid + ro - 12]) / (1 << task.data.shift) : 0.0f;
float4 fc0 = vload4(0, &fcoef[0]);
float4 fc1 = vload4(1, &fcoef[0]);
float4 fc2 = vload4(2, &fcoef[0]);
float fdata[12];
for (int pos = 0; pos < 12; pos++)
fdata[pos] = pos + ro - 12 < 0 ? 0.0f : (float)(data[pos + ro - 12] >> task.data.wbits);
float4 fd0 = vload4(0, &fdata[0]);
float4 fd1 = vload4(1, &fdata[0]);
float4 fd2 = vload4(2, &fdata[0]);
for (int pos = ro; pos < bs; pos ++)
{
float4 sum4 = fc0 * fd0 + fc1 * fd1 + fc2 * fd2;
float2 sum2 = sum4.s01 + sum4.s23;
fd0 = fd0.s1230;
fd1 = fd1.s1230;
fd2 = fd2.s1230;
fd0.s3 = fd1.s3;
fd1.s3 = fd2.s3;
fd2.s3 = (float)(data[pos] >> task.data.wbits);
len[pos >> 6] += fabs(fd2.s3 + (sum2.x + sum2.y));
}
}
else
{
float fcoef[32];
for (int tid = 0; tid < 32; tid++)
fcoef[tid] = tid < MAX_ORDER && tid + ro - MAX_ORDER >= 0 ? - ((float) task.coefs[tid + ro - MAX_ORDER]) / (1 << task.data.shift) : 0.0f;
float4 fc0 = vload4(0, &fcoef[0]);
float4 fc1 = vload4(1, &fcoef[0]);
float4 fc2 = vload4(2, &fcoef[0]);
float fdata[MAX_ORDER + TEMPBLOCK1 + 32];
for (int pos = 0; pos < MAX_ORDER; pos++)
fdata[pos] = 0.0f;
for (int pos = MAX_ORDER + TEMPBLOCK1; pos < MAX_ORDER + TEMPBLOCK1 + 32; pos++)
fdata[pos] = 0.0f;
for (int bpos = 0; bpos < bs; bpos += TEMPBLOCK1)
{
int end = min(bpos + TEMPBLOCK1, bs);
for (int pos = max(bpos - ro, 0); pos < max(bpos, ro); pos++)
fdata[MAX_ORDER + pos - bpos] = (float)(data[pos] >> task.data.wbits);
for (int pos = max(bpos, ro); pos < end; pos ++)
{
float next = (float)(data[pos] >> task.data.wbits);
float * dptr = fdata + pos - bpos;
dptr[MAX_ORDER] = next;
float4 sum
= fc0 * vload4(0, dptr)
+ fc1 * vload4(1, dptr)
#if MAX_ORDER > 8
+ fc2 * vload4(2, dptr)
#if MAX_ORDER > 12
+ vload4(3, &fcoef[0]) * vload4(3, dptr)
#if MAX_ORDER > 16
+ vload4(4, &fcoef[0]) * vload4(4, dptr)
+ vload4(5, &fcoef[0]) * vload4(5, dptr)
+ vload4(6, &fcoef[0]) * vload4(6, dptr)
+ vload4(7, &fcoef[0]) * vload4(7, dptr)
#endif
#endif
#endif
;
next += sum.x + sum.y + sum.z + sum.w;
len[pos >> 6] += fabs(next);
}
}
}
int total = 0;
for (int i = 0; i < ERPARTS; i++)
{
int res = convert_int_sat_rte(len[i] * 2);
int k = clamp(31 - clz(res) - 6, 0, MAX_RICE_PARAM); // 25 - clz(res) == clz(64) - clz(res) == log2(res / 64)
total += (k << 6) + (res >> k);
}
int partLen = min(0x7ffffff, total) + (bs - ro);
int obits = task.data.obits - task.data.wbits;
tasks[selectedTask].data.size = min(obits * bs,
task.data.type == Fixed ? ro * obits + 6 + RICE_PARAM_BITS + partLen :
task.data.type == LPC ? ro * obits + 4 + 5 + ro * task.data.cbits + 6 + RICE_PARAM_BITS/* << porder */ + partLen :
task.data.type == Constant ? obits * select(1, bs, partLen != bs - ro) :
obits * bs);
}
#else
#define ESTPARTLOG 5
__kernel /*__attribute__(( vec_type_hint (int4)))*/ __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clEstimateResidual(
__global int*samples,
__global int*selectedTasks,
__global FLACCLSubframeTask *tasks
)
{
__local float data[GROUP_SIZE * 2 + 32];
#if !defined(AMD)
__local volatile uint idata[GROUP_SIZE + 16];
#endif
__local FLACCLSubframeTask task;
__local uint psum[MAX_BLOCKSIZE >> ESTPARTLOG];
__local float fcoef[32];
__local int selectedTask;
if (get_local_id(0) == 0)
selectedTask = selectedTasks[get_group_id(0)];
barrier(CLK_LOCAL_MEM_FENCE);
const int tid = get_local_id(0);
if (tid < sizeof(task)/sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[selectedTask]))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
int ro = task.data.residualOrder;
int bs = task.data.blocksize;
if (tid < 32)
//fcoef[tid] = select(0.0f, - ((float) task.coefs[tid]) / (1 << task.data.shift), tid < ro);
fcoef[tid] = tid < MAX_ORDER && tid + ro - MAX_ORDER >= 0 ? - ((float) task.coefs[tid + ro - MAX_ORDER]) / (1 << task.data.shift) : 0.0f;
for (int offs = tid; offs < (MAX_BLOCKSIZE >> ESTPARTLOG); offs += GROUP_SIZE)
psum[offs] = 0;
data[tid] = 0.0f;
// need to initialize "extra" data, because NaNs can produce weird results even when multiplied by zero extra coefs
if (tid < 32)
data[GROUP_SIZE * 2 + tid] = 0.0f;
barrier(CLK_LOCAL_MEM_FENCE);
float4 fc0 = vload4(0, &fcoef[0]);
float4 fc1 = vload4(1, &fcoef[0]);
#if MAX_ORDER > 8
float4 fc2 = vload4(2, &fcoef[0]);
#endif
__global int * rptr = &samples[task.data.samplesOffs];
int wb = task.data.wbits;
int pos;
for (pos = 0; pos + GROUP_SIZE - 1 < bs; pos += GROUP_SIZE)
{
// fetch samples
int offs = pos + tid;
float nextData = rptr[offs] >> wb;
data[tid + GROUP_SIZE] = nextData;
barrier(CLK_LOCAL_MEM_FENCE);
// compute residual
__local float* dptr = &data[tid + GROUP_SIZE - MAX_ORDER];
float4 sum4
= fc0 * vload4(0, dptr)
+ fc1 * vload4(1, dptr)
#if MAX_ORDER > 8
+ fc2 * vload4(2, dptr)
#if MAX_ORDER > 12
+ vload4(3, &fcoef[0]) * vload4(3, dptr)
#if MAX_ORDER > 16
+ vload4(4, &fcoef[0]) * vload4(4, dptr)
+ vload4(5, &fcoef[0]) * vload4(5, dptr)
+ vload4(6, &fcoef[0]) * vload4(6, dptr)
+ vload4(7, &fcoef[0]) * vload4(7, dptr)
#endif
#endif
#endif
;
float2 sum2 = sum4.s01 + sum4.s23;
int it = convert_int_sat_rte(nextData + (sum2.s0 + sum2.s1));
// int t = (int)(nextData + sum.x + sum.y + sum.z + sum.w);
barrier(CLK_LOCAL_MEM_FENCE);
data[tid] = nextData;
// convert to unsigned
uint t = (it << 1) ^ (it >> 31);
// ensure we're within frame bounds
t = select(0U, t, offs >= ro);
// overflow protection
t = min(t, 0x7ffffffU);
#if defined(AMD)
atomic_add(&psum[min(MAX_BLOCKSIZE - 1, offs) >> ESTPARTLOG], t);
#else
idata[tid] = t;
#if WARP_SIZE <= (1 << (ESTPARTLOG - 1))
barrier(CLK_LOCAL_MEM_FENCE);
for (int l = 1 << (ESTPARTLOG - 1); l >= WARP_SIZE; l >>= 1) {
if (!(tid & l)) idata[tid] += idata[tid + l];
barrier(CLK_LOCAL_MEM_FENCE);
}
for (int l = WARP_SIZE / 2; l > 1; l >>= 1)
idata[tid] += idata[tid + l];
#else
for (int l = 1 << (ESTPARTLOG - 1); l > 1; l >>= 1)
idata[tid] += idata[tid + l];
#endif
if ((tid & (1 << ESTPARTLOG) - 1) == 0)
psum[min(MAX_BLOCKSIZE - 1, offs) >> ESTPARTLOG] = idata[tid] + idata[tid + 1];
#endif
}
if (pos < bs)
{
// fetch samples
int offs = pos + tid;
float nextData = offs < bs ? rptr[offs] >> wb : 0.0f;
data[tid + GROUP_SIZE] = nextData;
barrier(CLK_LOCAL_MEM_FENCE);
// compute residual
__local float* dptr = &data[tid + GROUP_SIZE - MAX_ORDER];
float4 sum
= fc0 * vload4(0, dptr)
+ fc1 * vload4(1, dptr)
#if MAX_ORDER > 8
+ fc2 * vload4(2, dptr)
#if MAX_ORDER > 12
+ vload4(3, &fcoef[0]) * vload4(3, dptr)
#if MAX_ORDER > 16
+ vload4(4, &fcoef[0]) * vload4(4, dptr)
+ vload4(5, &fcoef[0]) * vload4(5, dptr)
+ vload4(6, &fcoef[0]) * vload4(6, dptr)
+ vload4(7, &fcoef[0]) * vload4(7, dptr)
#endif
#endif
#endif
;
int it = convert_int_sat_rte(nextData + sum.x + sum.y + sum.z + sum.w);
barrier(CLK_LOCAL_MEM_FENCE);
data[tid] = nextData;
// convert to unsigned
uint t = (it << 1) ^ (it >> 31);
// ensure we're within frame bounds
t = select(0U, t, offs >= ro && offs < bs);
// overflow protection
t = min(t, 0x7ffffffU);
#if defined(AMD)
atomic_add(&psum[min(MAX_BLOCKSIZE - 1, offs) >> ESTPARTLOG], t);
#else
idata[tid] = t;
#if WARP_SIZE <= (1 << (ESTPARTLOG - 1))
barrier(CLK_LOCAL_MEM_FENCE);
for (int l = 1 << (ESTPARTLOG - 1); l >= WARP_SIZE; l >>= 1) {
if (!(tid & l)) idata[tid] += idata[tid + l];
barrier(CLK_LOCAL_MEM_FENCE);
}
for (int l = WARP_SIZE / 2; l > 1; l >>= 1)
idata[tid] += idata[tid + l];
#else
for (int l = 1 << (ESTPARTLOG - 1); l > 1; l >>= 1)
idata[tid] += idata[tid + l];
#endif
if ((tid & (1 << ESTPARTLOG) - 1) == 0)
psum[min(MAX_BLOCKSIZE - 1, offs) >> ESTPARTLOG] = idata[tid] + idata[tid + 1];
#endif
}
// calculate rice partition bit length for every 32 samples
barrier(CLK_LOCAL_MEM_FENCE);
#if (MAX_BLOCKSIZE >> (ESTPARTLOG + 1)) > GROUP_SIZE
#error MAX_BLOCKSIZE is too large for this GROUP_SIZE
#endif
uint pl = tid < (MAX_BLOCKSIZE >> (ESTPARTLOG + 1)) ? psum[tid * 2] + psum[tid * 2 + 1] : 0;
barrier(CLK_LOCAL_MEM_FENCE);
// for (int pos = 0; pos < (MAX_BLOCKSIZE >> ESTPARTLOG) / 2; pos += GROUP_SIZE)
// {
//int offs = pos + tid;
//int pl = offs < (MAX_BLOCKSIZE >> ESTPARTLOG) / 2 ? psum[offs * 2] + psum[offs * 2 + 1] : 0;
////int pl = psum[offs * 2] + psum[offs * 2 + 1];
//barrier(CLK_LOCAL_MEM_FENCE);
//if (offs < (MAX_BLOCKSIZE >> ESTPARTLOG) / 2)
// psum[offs] = pl;
// }
int k = clamp(31 - (int)clz(pl) - (ESTPARTLOG + 1), 0, MAX_RICE_PARAM); // 26 - clz(res) == clz(32) - clz(res) == log2(res / 32)
if (tid < MAX_BLOCKSIZE >> (ESTPARTLOG + 1))
psum[tid] = (k << (ESTPARTLOG + 1)) + (pl >> k);
barrier(CLK_LOCAL_MEM_FENCE);
for (int l = MAX_BLOCKSIZE >> (ESTPARTLOG + 2); l > 0; l >>= 1)
{
if (tid < l)
psum[tid] += psum[tid + l];
barrier(CLK_LOCAL_MEM_FENCE);
}
if (tid == 0)
{
int pl = (int)psum[0] + (bs - ro);
int obits = task.data.obits - task.data.wbits;
int len = min(obits * task.data.blocksize,
task.data.type == Fixed ? task.data.residualOrder * obits + 6 + RICE_PARAM_BITS + pl :
task.data.type == LPC ? task.data.residualOrder * obits + 4 + 5 + task.data.residualOrder * task.data.cbits + 6 + RICE_PARAM_BITS/* << porder */ + pl :
task.data.type == Constant ? obits * select(1, task.data.blocksize, pl != task.data.blocksize - task.data.residualOrder) :
obits * task.data.blocksize);
tasks[selectedTask].data.size = len;
}
}
#endif
__kernel
void clSelectStereoTasks(
__global FLACCLSubframeTask *tasks,
__global int*selectedTasks,
__global int*selectedTasksSecondEstimate,
__global int*selectedTasksBestMethod,
int tasksWindow,
int windowCount,
int tasksToSecondEstimate,
int taskCount,
int selectedCount
)
{
int best_size[4];
int best_wind[4];
for (int ch = 0; ch < 4; ch++)
{
int first_no = selectedTasks[(get_global_id(0) * 4 + ch) * selectedCount];
int best_len = tasks[first_no].data.size;
int best_wnd = 0;
for (int i = 1; i < selectedCount; i++)
{
int task_no = selectedTasks[(get_global_id(0) * 4 + ch) * selectedCount + i];
int task_len = tasks[task_no].data.size;
int task_wnd = (task_no - first_no) / tasksWindow;
task_wnd = select(0, task_wnd, task_wnd < windowCount);
best_wnd = select(best_wnd, task_wnd, task_len < best_len);
best_len = min(task_len, best_len);
}
best_size[ch] = best_len;
best_wind[ch] = best_wnd;
}
int bitsBest = best_size[2] + best_size[3]; // MidSide
int chMask = 2 | (3 << 2);
int bits = best_size[3] + best_size[1];
chMask = select(chMask, 3 | (1 << 2), bits < bitsBest); // RightSide
bitsBest = min(bits, bitsBest);
bits = best_size[0] + best_size[3];
chMask = select(chMask, 0 | (3 << 2), bits < bitsBest); // LeftSide
bitsBest = min(bits, bitsBest);
bits = best_size[0] + best_size[1];
chMask = select(chMask, 0 | (1 << 2), bits < bitsBest); // LeftRight
bitsBest = min(bits, bitsBest);
for (int ich = 0; ich < 2; ich++)
{
int ch = select(chMask & 3, chMask >> 2, ich > 0);
int roffs = tasks[(get_global_id(0) * 4 + ich) * taskCount].data.samplesOffs;
int nonSelectedNo = 0;
for (int j = taskCount - 1; j >= 0; j--)
{
int i = select(j, (j % windowCount) * tasksWindow + (j / windowCount), j < windowCount * tasksWindow);
int no = (get_global_id(0) * 4 + ch) * taskCount + i;
selectedTasksBestMethod[(get_global_id(0) * 2 + ich) * taskCount + i] = no;
tasks[no].data.residualOffs = roffs;
if (j >= selectedCount)
tasks[no].data.size = 0x7fffffff;
if (nonSelectedNo < tasksToSecondEstimate)
if (tasksToSecondEstimate == taskCount - selectedCount || best_wind[ch] == i / tasksWindow || i >= windowCount * tasksWindow)
selectedTasksSecondEstimate[(get_global_id(0) * 2 + ich) * tasksToSecondEstimate + nonSelectedNo++] = no;
}
}
}
__kernel
void clChooseBestMethod(
__global FLACCLSubframeTask *tasks_out,
__global FLACCLSubframeTask *tasks,
__global int*selectedTasks,
int taskCount
)
{
int best_no = selectedTasks[get_global_id(0) * taskCount];
int best_len = tasks[best_no].data.size;
for (int i = 1; i < taskCount; i++)
{
int task_no = selectedTasks[get_global_id(0) * taskCount + i];
int task_len = tasks[task_no].data.size;
best_no = select(best_no, task_no, task_len < best_len);
best_len = min(best_len, task_len);
}
tasks_out[get_global_id(0)] = tasks[best_no];
}
#ifdef DO_PARTITIONS
#if BITS_PER_SAMPLE > 16
#define residual_t long
#define convert_bps_sat convert_int_sat
#else
#define residual_t int
#define convert_bps_sat
#endif
#ifdef FLACCL_CPU
inline residual_t calc_residual(__global int *ptr, int * coefs, int ro)
{
residual_t sum = 0;
for (int i = 0; i < ro; i++)
sum += (residual_t)ptr[i] * coefs[i];
//sum += upsample(mul_hi(ptr[i], coefs[i]), as_uint(ptr[i] * coefs[i]));
return sum;
}
#define ENCODE_N(cro,action) for (int pos = cro; pos < bs; pos ++) { \
residual_t t = (data[pos] - (calc_residual(data + pos - cro, task.coefs, cro) >> task.data.shift)) >> task.data.wbits; \
action; \
}
#define SWITCH_N(action) \
switch (ro) \
{ \
case 0: ENCODE_N(0, action) break; \
case 1: ENCODE_N(1, action) break; \
case 2: ENCODE_N(2, action) /*if (task.coefs[0] == -1 && task.coefs[1] == 2) ENCODE_N(2, 2 * ptr[1] - ptr[0], action) else*/ break; \
case 3: ENCODE_N(3, action) break; \
case 4: ENCODE_N(4, action) break; \
case 5: ENCODE_N(5, action) break; \
case 6: ENCODE_N(6, action) break; \
case 7: ENCODE_N(7, action) break; \
case 8: ENCODE_N(8, action) break; \
case 9: ENCODE_N(9, action) break; \
case 10: ENCODE_N(10, action) break; \
case 11: ENCODE_N(11, action) break; \
case 12: ENCODE_N(12, action) break; \
default: ENCODE_N(ro, action) \
}
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clEncodeResidual(
__global ulong *partition_lengths,
__global int *residual,
__global int *samples,
__global FLACCLSubframeTask *tasks,
int max_porder, // <= 8
int psize // == task.blocksize >> max_porder?
)
{
FLACCLSubframeTask task = tasks[get_group_id(0)];
int bs = task.data.blocksize;
int ro = task.data.residualOrder;
__global int *data = &samples[task.data.samplesOffs];
__global ulong *pl = partition_lengths + (1 << (max_porder + 1)) * get_group_id(0);
int r;
for (int p = 0; p < (1 << max_porder); p++)
pl[p] = 0UL;
__global int *rptr = residual + task.data.residualOffs;
if (psize == 16)
{
SWITCH_N((rptr[pos] = r = convert_bps_sat(t), pl[pos >> 4] += (uint)((r << 1) ^ (r >> 31))));
}
else
{
SWITCH_N((rptr[pos] = r = convert_bps_sat(t), pl[pos / psize] += (uint)((r << 1) ^ (r >> 31))));
}
}
#else
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clEncodeResidual(
__global int *partition_lengths,
__global int *output,
__global int *samples,
__global FLACCLSubframeTask *tasks,
int max_porder, // <= 8
int psize // == task.blocksize >> max_porder?
)
{
__local FLACCLSubframeTask task;
__local int data[GROUP_SIZE * 2 + MAX_ORDER];
const int tid = get_local_id(0);
if (get_local_id(0) < sizeof(task) / sizeof(int))
((__local int*)&task)[get_local_id(0)] = ((__global int*)(&tasks[get_group_id(0)]))[get_local_id(0)];
barrier(CLK_LOCAL_MEM_FENCE);
int bs = task.data.blocksize;
int ro = task.data.residualOrder;
if (tid < 32 && tid >= ro)
task.coefs[tid] = 0;
barrier(CLK_LOCAL_MEM_FENCE);
int4 cptr0 = vload4(0, &task.coefs[0]);
int4 cptr1 = vload4(1, &task.coefs[0]);
#if MAX_ORDER > 8
int4 cptr2 = vload4(2, &task.coefs[0]);
#endif
// We tweaked coeffs so that (task.cbits + task.obits + log2i(ro) <= 32)
// when BITS_PER_SAMPLE == 16, so we don't need 64bit arithmetics.
data[tid] = 0;
for (int pos = 0; pos < bs; pos += GROUP_SIZE)
{
// fetch samples
int off = pos + tid;
int nextData = off < bs ? samples[task.data.samplesOffs + off] >> task.data.wbits : 0;
data[tid + GROUP_SIZE] = nextData;
barrier(CLK_LOCAL_MEM_FENCE);
// compute residual
__local int* dptr = &data[tid + GROUP_SIZE - ro];
#if BITS_PER_SAMPLE > 16
long4 sum
= upsample(mul_hi(cptr0, vload4(0, dptr)), as_uint4(cptr0 * vload4(0, dptr)))
+ upsample(mul_hi(cptr1, vload4(1, dptr)), as_uint4(cptr1 * vload4(1, dptr)))
#if MAX_ORDER > 8
+ upsample(mul_hi(cptr2, vload4(2, dptr)), as_uint4(cptr2 * vload4(2, dptr)))
#if MAX_ORDER > 12
+ upsample(mul_hi(vload4(3, &task.coefs[0]), vload4(3, dptr)), as_uint4(vload4(3, &task.coefs[0]) * vload4(3, dptr)))
#if MAX_ORDER > 16
+ upsample(mul_hi(vload4(4, &task.coefs[0]), vload4(4, dptr)), as_uint4(vload4(4, &task.coefs[0]) * vload4(4, dptr)))
+ upsample(mul_hi(vload4(5, &task.coefs[0]), vload4(5, dptr)), as_uint4(vload4(5, &task.coefs[0]) * vload4(5, dptr)))
+ upsample(mul_hi(vload4(6, &task.coefs[0]), vload4(6, dptr)), as_uint4(vload4(6, &task.coefs[0]) * vload4(6, dptr)))
+ upsample(mul_hi(vload4(7, &task.coefs[0]), vload4(7, dptr)), as_uint4(vload4(7, &task.coefs[0]) * vload4(7, dptr)))
#endif
#endif
#endif
#else
int4 sum
= cptr0 * vload4(0, dptr)
+ cptr1 * vload4(1, dptr)
#if MAX_ORDER > 8
+ cptr2 * vload4(2, dptr)
#if MAX_ORDER > 12
+ vload4(3, &task.coefs[0]) * vload4(3, dptr)
#if MAX_ORDER > 16
+ vload4(4, &task.coefs[0]) * vload4(4, dptr)
+ vload4(5, &task.coefs[0]) * vload4(5, dptr)
+ vload4(6, &task.coefs[0]) * vload4(6, dptr)
+ vload4(7, &task.coefs[0]) * vload4(7, dptr)
#endif
#endif
#endif
#endif
;
if (off >= ro && off < bs)
output[task.data.residualOffs + off] = convert_bps_sat(nextData - ((sum.x + sum.y + sum.z + sum.w) >> task.data.shift));
barrier(CLK_LOCAL_MEM_FENCE);
data[tid] = nextData;
}
}
#endif
#ifndef FLACCL_CPU
// get_group_id(0) == partition index / (GROUP_SIZE / 16)
// get_group_id(1) == task index
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clCalcPartition(
__global int *partition_lengths,
__global int *residual,
__global FLACCLSubframeTask *tasks,
int max_porder, // <= 8
int psize // == task.blocksize >> max_porder?
)
{
__local uint pl[(GROUP_SIZE / 16)][MAX_RICE_PARAM + 1];
__local FLACCLSubframeData task;
const int tid = get_local_id(0);
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[get_group_id(1)]))[tid];
if (tid < (GROUP_SIZE / 16))
{
for (int k = 0; k <= MAX_RICE_PARAM; k++)
pl[tid][k] = 0;
}
barrier(CLK_LOCAL_MEM_FENCE);
int start = get_group_id(0) * psize * (GROUP_SIZE / 16);
int end = min(start + psize * (GROUP_SIZE / 16), task.blocksize);
for (int offs = start + tid; offs < end; offs += GROUP_SIZE)
{
// fetch residual
int s = (offs >= task.residualOrder && offs < end) ? residual[task.residualOffs + offs] : 0;
// convert to unsigned
uint t = (s << 1) ^ (s >> 31);
// calc number of unary bits for each residual sample with each rice parameter
int part = (offs - start) / psize;
// we must ensure that psize * (t >> k) doesn't overflow;
uint lim = 0x7fffffffU / (uint)psize;
for (int k = 0; k <= MAX_RICE_PARAM; k++)
atomic_add(&pl[part][k], min(lim, t >> k));
//pl[part][k] += s >> k;
}
barrier(CLK_LOCAL_MEM_FENCE);
int part = get_group_id(0) * (GROUP_SIZE / 16) + tid;
if (tid < (GROUP_SIZE / 16) && part < (1 << max_porder))
{
for (int k = 0; k <= MAX_RICE_PARAM; k++)
{
// output length
const int pos = ((MAX_RICE_PARAM + 1) << (max_porder + 1)) * get_group_id(1) + (k << (max_porder + 1));
uint plen = pl[tid][k];
partition_lengths[pos + part] = min(0x007fffffU, plen) + (uint)(psize - select(0, task.residualOrder, part == 0)) * (k + 1);
// if (get_group_id(1) == 0)
//printf("pl[%d][%d] == %d\n", k, part, min(0x7fffff, pl[k][tid]) + (psize - task.residualOrder * (part == 0)) * (k + 1));
}
}
}
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clCalcPartition16(
__global unsigned int *partition_lengths,
__global int *residual,
__global FLACCLSubframeTask *tasks,
int max_porder // <= 8
)
{
__local FLACCLSubframeData task;
__local unsigned int res[GROUP_SIZE];
__local unsigned int pl[GROUP_SIZE >> 4][MAX_RICE_PARAM + 1];
const int tid = get_local_id(0);
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[get_group_id(0)]))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
int bs = task.blocksize;
int ro = task.residualOrder;
barrier(CLK_LOCAL_MEM_FENCE);
for (int pos = 0; pos < bs; pos += GROUP_SIZE)
{
int offs = pos + tid;
// fetch residual
int s = (offs >= ro && offs < bs) ? residual[task.residualOffs + offs] : 0;
// convert to unsigned
res[tid] = (s << 1) ^ (s >> 31);
barrier(CLK_LOCAL_MEM_FENCE);
// we must ensure that psize * (t >> k) doesn't overflow;
uint4 lim = 0x07ffffffU;
int x = tid >> 4;
__local uint * chunk = &res[x << 4];
for (int k0 = 0; k0 <= MAX_RICE_PARAM; k0 += 16)
{
// calc number of unary bits for each group of 16 residual samples
// with each rice parameter.
int k = k0 + (tid & 15);
uint4 rsum
= min(lim, vload4(0,chunk) >> k)
+ min(lim, vload4(1,chunk) >> k)
+ min(lim, vload4(2,chunk) >> k)
+ min(lim, vload4(3,chunk) >> k)
;
uint rs = rsum.x + rsum.y + rsum.z + rsum.w;
// We can safely limit length here to 0x007fffffU, not causing length
// mismatch, because any such length would cause Verbatim frame anyway.
// And this limit protects us from overflows when calculating larger
// partitions, as we can have a maximum of 2^8 partitions, resulting
// in maximum partition length of 0x7fffffffU + change.
if (k <= MAX_RICE_PARAM) pl[x][k] = min(0x007fffffU, rs) + (uint)(16 - select(0, ro, offs < 16)) * (k + 1);
}
barrier(CLK_LOCAL_MEM_FENCE);
for (int k0 = 0; k0 <= MAX_RICE_PARAM; k0 += 16)
{
int k1 = k0 + (tid >> (GROUP_SIZE_LOG - 4)), x1 = tid & ((1 << (GROUP_SIZE_LOG - 4)) - 1);
if (k1 <= MAX_RICE_PARAM && (pos >> 4) + x1 < (1 << max_porder))
partition_lengths[((MAX_RICE_PARAM + 1) << (max_porder + 1)) * get_group_id(0) + (k1 << (max_porder + 1)) + (pos >> 4) + x1] = pl[x1][k1];
}
}
}
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clCalcPartition32(
__global unsigned int *partition_lengths,
__global int *residual,
__global FLACCLSubframeTask *tasks,
int max_porder // <= 8
)
{
__local FLACCLSubframeData task;
__local unsigned int res[GROUP_SIZE];
__local unsigned int pl[GROUP_SIZE >> 5][32];
const int tid = get_local_id(0);
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[get_group_id(0)]))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
int bs = task.blocksize;
int ro = task.residualOrder;
barrier(CLK_LOCAL_MEM_FENCE);
for (int pos = 0; pos < bs; pos += GROUP_SIZE)
{
int offs = pos + tid;
// fetch residual
int s = (offs >= ro && offs < bs) ? residual[task.residualOffs + offs] : 0;
// convert to unsigned
res[tid] = (s << 1) ^ (s >> 31);
barrier(CLK_LOCAL_MEM_FENCE);
// we must ensure that psize * (t >> k) doesn't overflow;
uint4 lim = 0x03ffffffU;
int x = tid >> 5;
__local uint * chunk = &res[x << 5];
// calc number of unary bits for each group of 32 residual samples
// with each rice parameter.
int k = tid & 31;
uint4 rsum
= min(lim, vload4(0,chunk) >> k)
+ min(lim, vload4(1,chunk) >> k)
+ min(lim, vload4(2,chunk) >> k)
+ min(lim, vload4(3,chunk) >> k)
+ min(lim, vload4(4,chunk) >> k)
+ min(lim, vload4(5,chunk) >> k)
+ min(lim, vload4(6,chunk) >> k)
+ min(lim, vload4(7,chunk) >> k)
;
uint rs = rsum.x + rsum.y + rsum.z + rsum.w;
// We can safely limit length here to 0x007fffffU, not causing length
// mismatch, because any such length would cause Verbatim frame anyway.
// And this limit protects us from overflows when calculating larger
// partitions, as we can have a maximum of 2^8 partitions, resulting
// in maximum partition length of 0x7fffffffU + change.
if (k <= MAX_RICE_PARAM) pl[x][k] = min(0x007fffffU, rs) + (uint)(32 - select(0, ro, offs < 32)) * (k + 1);
barrier(CLK_LOCAL_MEM_FENCE);
int k1 = (tid >> (GROUP_SIZE_LOG - 5)), x1 = tid & ((1 << (GROUP_SIZE_LOG - 5)) - 1);
if (k1 <= MAX_RICE_PARAM && (pos >> 5) + x1 < (1 << max_porder))
partition_lengths[((MAX_RICE_PARAM + 1) << (max_porder + 1)) * get_group_id(0) + (k1 << (max_porder + 1)) + (pos >> 5) + x1] = pl[x1][k1];
}
}
#endif
#ifdef FLACCL_CPU
// Sums partition lengths for a certain k == get_group_id(0)
// get_group_id(0) == k
// get_group_id(1) == task index
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clSumPartition(
__global ulong* partition_lengths,
int max_porder
)
{
if (get_group_id(0) != 0) // ignore k != 0
return;
__global ulong * sums = partition_lengths + (1 << (max_porder + 1)) * get_group_id(1);
for (int i = max_porder - 1; i >= 0; i--)
{
for (int j = 0; j < (1 << i); j++)
{
sums[(2 << i) + j] = sums[2 * j] + sums[2 * j + 1];
// if (get_group_id(1) == 0)
//printf("[%d][%d]: %d + %d == %d\n", i, j, sums[2 * j], sums[2 * j + 1], sums[2 * j] + sums[2 * j + 1]);
}
sums += 2 << i;
}
}
#else
// Sums partition lengths for a certain k == get_group_id(0)
// Requires 128 threads
// get_group_id(0) == k
// get_group_id(1) == task index
__kernel __attribute__((reqd_work_group_size(128, 1, 1)))
void clSumPartition(
__global uint* partition_lengths,
int max_porder
)
{
__local uint data[256]; // max_porder <= 8, data length <= 1 << 9.
const int pos = ((MAX_RICE_PARAM + 1) << (max_porder + 1)) * get_group_id(1) + (get_group_id(0) << (max_porder + 1));
// fetch partition lengths
uint2 pl = get_local_id(0) * 2 < (1 << max_porder) ? vload2(get_local_id(0),&partition_lengths[pos]) : 0;
data[get_local_id(0)] = pl.x + pl.y;
barrier(CLK_LOCAL_MEM_FENCE);
int in_pos = (get_local_id(0) << 1);
int out_pos = (1 << (max_porder - 1)) + get_local_id(0);
for (int bs = 1 << (max_porder - 2); bs > 0; bs >>= 1)
{
if (get_local_id(0) < bs) data[out_pos] = data[in_pos] + data[in_pos + 1];
in_pos += bs << 1;
out_pos += bs;
barrier(CLK_LOCAL_MEM_FENCE);
}
if (get_local_id(0) < (1 << max_porder))
partition_lengths[pos + (1 << max_porder) + get_local_id(0)] = data[get_local_id(0)];
if (get_local_size(0) + get_local_id(0) < (1 << max_porder))
partition_lengths[pos + (1 << max_porder) + get_local_size(0) + get_local_id(0)] = data[get_local_size(0) + get_local_id(0)];
}
#endif
#ifdef FLACCL_CPU
// Finds optimal rice parameter for each partition.
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clFindRiceParameter(
__global FLACCLSubframeTask *tasks,
__global int* rice_parameters,
__global ulong* partition_lengths,
int max_porder
)
{
__global FLACCLSubframeTask* task = tasks + get_group_id(0);
const int tid = get_local_id(0);
int lim = (2 << max_porder) - 1;
//int psize = task->data.blocksize >> max_porder;
int bs = task->data.blocksize;
int ro = task->data.residualOrder;
__global ulong* ppl = &partition_lengths[get_group_id(0) << (max_porder + 1)];
__global int* prp = &rice_parameters[get_group_id(0) << (max_porder + 2)];
__global int* pol = prp + (1 << (max_porder + 1));
for (int porder = max_porder; porder >= 0; porder--)
{
int pos = (2 << max_porder) - (2 << porder);
int fin = pos + (1 << porder);
ulong pl = ppl[pos];
int ps = (bs >> porder) - ro;
int k = clamp(63 - (int)clz(pl / max(1, ps)), 0, MAX_RICE_PARAM);
int plk = ps * (k + 1) + (int)(pl >> k);
// output rice parameter
prp[pos] = k;
// output length
pol[pos] = plk;
ps = (bs >> porder);
for (int offs = pos + 1; offs < fin; offs++)
{
pl = ppl[offs];
k = clamp(63 - (int)clz(pl / ps), 0, MAX_RICE_PARAM);
plk = ps * (k + 1) + (int)(pl >> k);
// output rice parameter
prp[offs] = k;
// output length
pol[offs] = plk;
}
}
}
#else
// Finds optimal rice parameter for each partition.
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clFindRiceParameter(
__global FLACCLSubframeTask *tasks,
__global int* rice_parameters,
__global uint* partition_lengths,
int max_porder
)
{
for (int offs = get_local_id(0); offs < (2 << max_porder); offs += GROUP_SIZE)
{
const int pos = ((MAX_RICE_PARAM + 1) << (max_porder + 1)) * get_group_id(0) + offs;
uint best_l = partition_lengths[pos];
int best_k = 0;
for (int k = 1; k <= MAX_RICE_PARAM; k++)
{
uint l = partition_lengths[pos + (k << (max_porder + 1))];
best_k = select(best_k, k, l < best_l);
best_l = min(best_l, l);
}
// output rice parameter
rice_parameters[(get_group_id(0) << (max_porder + 2)) + offs] = best_k;
// output length
rice_parameters[(get_group_id(0) << (max_porder + 2)) + (1 << (max_porder + 1)) + offs] = best_l;
}
}
#endif
#ifdef FLACCL_CPU
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clFindPartitionOrder(
__global int *residual,
__global int* best_rice_parameters,
__global FLACCLSubframeTask *tasks,
__global int* rice_parameters,
int max_porder
)
{
__global FLACCLSubframeTask* task = tasks + get_group_id(0);
int partlen[9];
for (int p = 0; p < 9; p++)
partlen[p] = 0;
// fetch partition lengths
const int pos = (get_group_id(0) << (max_porder + 2)) + (2 << max_porder);
for (int porder = max_porder; porder >= 0; porder--)
{
int start = (2 << max_porder) - (2 << porder);
for (int offs = 0; offs < (1 << porder); offs ++)
partlen[porder] += rice_parameters[pos + start + offs];
}
int best_length = partlen[0] + RICE_PARAM_BITS;
int best_porder = 0;
for (int porder = 1; porder <= max_porder; porder++)
{
int length = (RICE_PARAM_BITS << porder) + partlen[porder];
best_porder = select(best_porder, porder, length < best_length);
best_length = min(best_length, length);
}
best_length = (RICE_PARAM_BITS << best_porder) + task->data.blocksize - task->data.residualOrder;
int best_psize = task->data.blocksize >> best_porder;
int start = task->data.residualOffs + task->data.residualOrder;
int fin = task->data.residualOffs + best_psize;
for (int p = 0; p < (1 << best_porder); p++)
{
int k = rice_parameters[pos - (2 << best_porder) + p];
best_length += k * (fin - start);
for (int i = start; i < fin; i++)
{
int t = residual[i];
best_length += ((t << 1) ^ (t >> 31)) >> k;
}
start = fin;
fin += best_psize;
}
int obits = task->data.obits - task->data.wbits;
task->data.porder = best_porder;
task->data.headerLen =
task->data.type == Constant ? obits :
task->data.type == Verbatim ? obits * task->data.blocksize :
task->data.type == Fixed ? task->data.residualOrder * obits + 6 :
task->data.type == LPC ? task->data.residualOrder * obits + 6 + 4 + 5 + task->data.residualOrder * task->data.cbits : 0;
task->data.size =
task->data.headerLen + ((task->data.type == Fixed || task->data.type == LPC) ? best_length : 0);
if (task->data.size >= obits * task->data.blocksize)
{
task->data.headerLen = task->data.size = obits * task->data.blocksize;
task->data.type = Verbatim;
}
for (int offs = 0; offs < (1 << best_porder); offs ++)
best_rice_parameters[(get_group_id(0) << max_porder) + offs] = rice_parameters[pos - (2 << best_porder) + offs];
}
#else
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clFindPartitionOrder(
__global int *residual,
__global int* best_rice_parameters,
__global FLACCLSubframeTask *tasks,
__global int* rice_parameters,
int max_porder
)
{
__local int partlen[16];
__local FLACCLSubframeData task;
const int pos = (get_group_id(0) << (max_porder + 2)) + (2 << max_porder);
if (get_local_id(0) < sizeof(task) / sizeof(int))
((__local int*)&task)[get_local_id(0)] = ((__global int*)(&tasks[get_group_id(0)]))[get_local_id(0)];
if (get_local_id(0) < 16)
partlen[get_local_id(0)] = 0;
barrier(CLK_LOCAL_MEM_FENCE);
// fetch partition lengths
int lim = (2 << max_porder) - 1;
for (int offs = get_local_id(0); offs < lim; offs += GROUP_SIZE)
{
int len = rice_parameters[pos + offs];
int porder = 31 - clz(lim - offs);
atomic_add(&partlen[porder], len);
}
barrier(CLK_LOCAL_MEM_FENCE);
int best_length = partlen[0] + RICE_PARAM_BITS;
int best_porder = 0;
for (int porder = 1; porder <= max_porder; porder++)
{
int length = (RICE_PARAM_BITS << porder) + partlen[porder];
best_porder = select(best_porder, porder, length < best_length);
best_length = min(best_length, length);
}
if (get_local_id(0) == 0)
{
task.porder = best_porder;
int obits = task.obits - task.wbits;
task.headerLen =
task.type == Fixed ? task.residualOrder * obits + 6 :
task.type == LPC ? task.residualOrder * obits + 6 + 4 + 5 + task.residualOrder * task.cbits :
task.type == Constant ? obits :
/* task.type == Verbatim ? */ obits * task.blocksize;
task.size = task.headerLen + select(0, best_length, task.type == Fixed || task.type == LPC);
if (task.size >= obits * task.blocksize)
{
task.headerLen = task.size = obits * task.blocksize;
task.type = Verbatim;
}
}
barrier(CLK_LOCAL_MEM_FENCE);
if (get_local_id(0) < sizeof(task) / sizeof(int))
((__global int*)(&tasks[get_group_id(0)]))[get_local_id(0)] = ((__local int*)&task)[get_local_id(0)];
for (int offs = get_local_id(0); offs < (1 << best_porder); offs += GROUP_SIZE)
best_rice_parameters[(get_group_id(0) << max_porder) + offs] = rice_parameters[pos - (2 << best_porder) + offs];
// FIXME: should be bytes?
}
#endif
#ifdef DO_RICE
#ifdef FLACCL_CPU
typedef struct BitWriter_t
{
__global unsigned int *buffer;
unsigned int bit_buf;
int bit_left;
int buf_ptr;
} BitWriter;
inline void writebits(BitWriter *bw, int bits, int v)
{
uint val = ((uint)v) & ((1 << bits) - 1);
if (bits < bw->bit_left)
{
bw->bit_buf = (bw->bit_buf << bits) | val;
bw->bit_left -= bits;
}
else
{
// if (bits >= 32) printf("\n\n\n\n-------------------------\n\n\n");
unsigned int bb = (bw->bit_buf << bw->bit_left) | (val >> (bits - bw->bit_left));
bw->buffer[bw->buf_ptr++] = (bb >> 24) | ((bb >> 8) & 0xff00) | ((bb << 8) & 0xff0000) | ((bb << 24) & 0xff000000);
bw->bit_left += (32 - bits);
bw->bit_buf = val;
// bw->bit_buf = val & ((1 << (32 - bw->bit_left)) - 1);
}
}
inline void flush(BitWriter *bw)
{
if (bw->bit_left < 32)
writebits(bw, bw->bit_left, 0);
}
#endif
inline int len_utf8(int n)
{
int bts = 31 - clz(n);
return select(8, 8 * ((bts + 4) / 5), bts > 6);
}
#ifdef FLACCL_CPU
// get_global_id(0) * channels == task index
__kernel __attribute__((reqd_work_group_size(1, 1, 1)))
void clCalcOutputOffsets(
__global int *residual,
__global int *samples,
__global FLACCLSubframeTask *tasks,
int channels,
int frameCount,
int firstFrame
)
{
int offset = 0;
for (int iFrame = 0; iFrame < frameCount; iFrame++)
{
//printf("len_utf8(%d) == %d\n", firstFrame + iFrame, len_utf8(firstFrame + iFrame));
offset += 15 + 1 + 4 + 4 + 4 + 3 + 1 + len_utf8(firstFrame + iFrame)
// + 8-16 // custom block size
// + 8-16 // custom sample rate
;
int bs = tasks[iFrame * channels].data.blocksize;
//public static readonly int[] flac_blocksizes = new int[15] { 0, 192, 576, 1152, 2304, 4608, 0, 0, 256, 512, 1024, 2048, 4096, 8192, 16384 };
offset += select(0, select(8, 16, bs >= 256), bs != 4096 && bs != 4608); // TODO: check all other standard sizes
// assert (offset % 8) == 0
offset += 8;
for (int ch = 0; ch < channels; ch++)
{
__global FLACCLSubframeTask* task = tasks + iFrame * channels + ch;
offset += 8 + task->data.wbits;
// Add 32 bits to separate frames if header is too small so they can intersect
offset += 64;
task->data.encodingOffset = offset + task->data.headerLen;
offset += task->data.size;
}
offset = (offset + 7) & ~7;
offset += 16;
}
}
#else
// get_global_id(0) * channels == task index
__kernel __attribute__((reqd_work_group_size(32, 1, 1)))
void clCalcOutputOffsets(
__global int *residual,
__global int *samples,
__global FLACCLSubframeTask *tasks,
int channels1,
int frameCount,
int firstFrame
)
{
__local FLACCLSubframeData ltasks[MAX_CHANNELS];
__local volatile int mypos[MAX_CHANNELS];
int offset = 0;
for (int iFrame = 0; iFrame < frameCount; iFrame++)
{
if (get_local_id(0) < sizeof(ltasks[0]) / sizeof(int))
for (int ch = 0; ch < MAX_CHANNELS; ch++)
((__local int*)&ltasks[ch])[get_local_id(0)] = ((__global int*)(&tasks[iFrame * MAX_CHANNELS + ch]))[get_local_id(0)];
//printf("len_utf8(%d) == %d\n", firstFrame + iFrame, len_utf8(firstFrame + iFrame));
offset += 15 + 1 + 4 + 4 + 4 + 3 + 1 + len_utf8(firstFrame + iFrame)
// + 8-16 // custom block size
// + 8-16 // custom sample rate
;
int bs = ltasks[0].blocksize;
//public static readonly int[] flac_blocksizes = new int[15] { 0, 192, 576, 1152, 2304, 4608, 0, 0, 256, 512, 1024, 2048, 4096, 8192, 16384 };
offset += select(0, select(8, 16, bs >= 256), bs != 4096 && bs != 4608); // TODO: check all other standard sizes
// assert (offset % 8) == 0
offset += 8;
if (get_local_id(0) < MAX_CHANNELS)
{
int ch = get_local_id(0);
// Add 64 bits to separate frames if header is too small so they can intersect
int mylen = 8 + ltasks[ch].wbits + 64 + ltasks[ch].size;
mypos[ch] = mylen;
for (int offset = 1; offset < WARP_SIZE && offset < MAX_CHANNELS; offset <<= 1)
if (ch >= offset) mypos[ch] += mypos[ch - offset];
mypos[ch] += offset;
tasks[iFrame * MAX_CHANNELS + ch].data.encodingOffset = mypos[ch] - ltasks[ch].size + ltasks[ch].headerLen;
}
offset = mypos[MAX_CHANNELS - 1];
offset = (offset + 7) & ~7;
offset += 16;
}
}
#endif
// get_group_id(0) == task index
__kernel __attribute__((reqd_work_group_size(GROUP_SIZE, 1, 1)))
void clRiceEncoding(
__global int *residual,
__global int *samples,
__global int* best_rice_parameters,
__global FLACCLSubframeTask *tasks,
__global unsigned int* output,
int max_porder
)
{
#ifdef FLACCL_CPU
__global FLACCLSubframeTask* task = tasks + get_group_id(0);
if (task->data.type == Fixed || task->data.type == LPC)
{
int ro = task->data.residualOrder;
int bs = task->data.blocksize;
int porder = task->data.porder;
int psize = bs >> porder;
BitWriter bw;
bw.buffer = output;
bw.buf_ptr = task->data.encodingOffset / 32;
bw.bit_left = 32 - (task->data.encodingOffset & 31);
bw.bit_buf = 0;
//if (get_group_id(0) == 0) printf("%d\n", offs);
int res_cnt = psize - ro;
// residual
int j = ro;
__global int * kptr = &best_rice_parameters[get_group_id(0) << max_porder];
for (int p = 0; p < (1 << porder); p++)
{
int k = kptr[p];
writebits(&bw, RICE_PARAM_BITS, k);
//if (get_group_id(0) == 0) printf("[%x] ", k);
//if (get_group_id(0) == 0) printf("(%x) ", bw.bit_buf);
if (p == 1) res_cnt = psize;
int cnt = min(res_cnt, bs - j);
unsigned int kexp = 1U << k;
__global int *rptr = &residual[task->data.residualOffs + j];
for (int i = 0; i < cnt; i++)
{
int iv = rptr[i];
unsigned int v = (iv << 1) ^ (iv >> 31);
// write quotient in unary
int bits = k + (v >> k) + 1;
while (bits > 31)
{
int b = min(bits - 31, 31);
if (b < bw.bit_left)
{
bw.bit_buf <<= b;
bw.bit_left -= b;
}
else
{
unsigned int bb = bw.bit_buf << bw.bit_left;
bw.bit_buf = 0;
bw.bit_left += (32 - b);
bw.buffer[bw.buf_ptr++] = as_int(as_char4(bb).wzyx);
}
bits -= b;
}
unsigned int val = (v & (kexp - 1)) | kexp;
if (bits < bw.bit_left)
{
bw.bit_buf = (bw.bit_buf << bits) | val;
bw.bit_left -= bits;
}
else
{
unsigned int bb = (bw.bit_buf << bw.bit_left) | (val >> (bits - bw.bit_left));
bw.bit_buf = val;
bw.bit_left += (32 - bits);
bw.buffer[bw.buf_ptr++] = as_int(as_char4(bb).wzyx);
}
////if (get_group_id(0) == 0) printf("%x ", v);
//writebits(&bw, (v >> k) + 1, 1);
////if (get_group_id(0) == 0) printf("(%x) ", bw.bit_buf);
//writebits(&bw, k, v);
////if (get_group_id(0) == 0) printf("(%x) ", bw.bit_buf);
}
j += cnt;
}
//if (bw.buf_ptr * 32 + 32 - bw.bit_left != task->data.encodingOffset - task->data.headerLen + task->data.size)
// printf("bit length mismatch: encodingOffset == %d, headerLen == %d, size == %d, so should be %d, but is %d\n",
// task->data.encodingOffset, task->data.headerLen, task->data.size,
// task->data.encodingOffset - task->data.headerLen + task->data.size,
// bw.buf_ptr * 32 + 32 - bw.bit_left
// );
//if (get_group_id(0) == 0) printf("\n");
flush(&bw);
}
#else
__local uint data[GROUP_SIZE];
__local volatile int mypos[GROUP_SIZE+1];
#if 0
__local int brp[256];
#endif
__local volatile int warppos[WARP_SIZE];
__local FLACCLSubframeData task;
int tid = get_local_id(0);
if (tid < sizeof(task) / sizeof(int))
((__local int*)&task)[tid] = ((__global int*)(&tasks[get_group_id(0)]))[tid];
barrier(CLK_LOCAL_MEM_FENCE);
if (task.type != Fixed && task.type != LPC)
return;
if (tid == 0)
mypos[GROUP_SIZE] = 0;
if (tid < WARP_SIZE)
warppos[tid] = 0;
#if 0
for (int offs = tid; offs < (1 << task.porder); offs ++)
brp[offs] = best_rice_parameters[(get_group_id(0) << max_porder) + offs];
#endif
data[tid] = 0;
barrier(CLK_LOCAL_MEM_FENCE);
const int bs = task.blocksize;
int start = task.encodingOffset;
int plen = bs >> task.porder;
//int plenoffs = 12 - task.porder;
uint remainder = 0U;
int pos;
for (pos = 0; pos + GROUP_SIZE - 1 < bs; pos += GROUP_SIZE)
{
int offs = pos + tid;
int iv = residual[task.residualOffs + offs];
int part = offs / plen;
//int part = offs >> plenoffs;
#if 0
int k = brp[part];
#else
int k = best_rice_parameters[(get_group_id(0) << max_porder) + part];
#endif
int pstart = offs == part * plen;
//int pstart = offs == part << plenoffs;
uint v = (iv << 1) ^ (iv >> 31);
int mylen = select(0, (int)(v >> k) + 1 + k, offs >= task.residualOrder && offs < bs) + select(0, RICE_PARAM_BITS, pstart);
mypos[tid] = mylen;
// Inclusive scan(+)
int lane = (tid & (WARP_SIZE - 1));
for (int offset = 1; offset < WARP_SIZE; offset <<= 1)
mypos[tid] += mypos[select(GROUP_SIZE, tid - offset, lane >= offset)];
int mp = mypos[tid];
if ((tid & (WARP_SIZE - 1)) == WARP_SIZE - 1)
warppos[tid/WARP_SIZE] = mp;
barrier(CLK_LOCAL_MEM_FENCE);
if (tid < GROUP_SIZE/WARP_SIZE)
{
for (int offset = 1; offset < GROUP_SIZE/WARP_SIZE; offset <<= 1)
warppos[tid] += warppos[select(GROUP_SIZE/WARP_SIZE, tid - offset, tid >= offset)];
}
barrier(CLK_LOCAL_MEM_FENCE);
mp += start + select(0, warppos[tid / WARP_SIZE - 1], tid / WARP_SIZE > 0);
int start32 = start >> 5;
start += mypos[GROUP_SIZE - 1] + warppos[GROUP_SIZE / WARP_SIZE - 2];
//if (start / 32 - start32 >= GROUP_SIZE - 3)
// tasks[get_group_id(0)].data.size = 1;
//if (tid == GROUP_SIZE - 1 && mypos[tid] > (GROUP_SIZE/2) * 32)
// printf("Oops: %d\n", mypos[tid]);
data[tid] = select(0U, remainder, tid == 0);
barrier(CLK_LOCAL_MEM_FENCE);
if (pstart)
{
int kpos = mp - mylen;
int kpos0 = (kpos >> 5) - start32;
int kpos1 = kpos & 31;
uint kval = (uint)k << (32 - RICE_PARAM_BITS);
uint kval0 = kval >> kpos1;
uint kval1 = kval << (32 - kpos1);
if (kval0) atomic_or(&data[kpos0], kval0);
if (kpos1 && kval1) atomic_or(&data[kpos0 + 1], kval1);
}
if (offs >= task.residualOrder && offs < bs)
{
int qpos = mp - k - 1;
int qpos0 = (qpos >> 5) - start32;
int qpos1 = qpos & 31;
uint qval = (1U << 31) | (v << (31 - k));
uint qval0 = qval >> qpos1;
uint qval1= qval << (32 - qpos1);
if (qval0) atomic_or(&data[qpos0], qval0);
if (qpos1 && qval1) atomic_or(&data[qpos0 + 1], qval1);
}
barrier(CLK_LOCAL_MEM_FENCE);
if ((start32 + tid) * 32 <= start)
output[start32 + tid] = as_int(as_char4(data[tid]).wzyx);
remainder = data[start / 32 - start32];
}
if (pos < bs)
{
int offs = pos + tid;
int iv = offs < bs ? residual[task.residualOffs + offs] : 0;
int part = offs / plen; // >> plenoffs;
//int k = brp[min(255, part)];
int k = offs < bs ? best_rice_parameters[(get_group_id(0) << max_porder) + part] : 0;
int pstart = offs == part * plen;
uint v = (iv << 1) ^ (iv >> 31);
int mylen = select(0, (int)(v >> k) + 1 + k, offs >= task.residualOrder && offs < bs) + select(0, RICE_PARAM_BITS, pstart);
mypos[tid] = mylen;
// Inclusive scan(+)
int lane = (tid & (WARP_SIZE - 1));
for (int offset = 1; offset < WARP_SIZE; offset <<= 1)
mypos[tid] += mypos[select(GROUP_SIZE, tid - offset, lane >= offset)];
int mp = mypos[tid];
if ((tid & (WARP_SIZE - 1)) == WARP_SIZE - 1)
warppos[tid/WARP_SIZE] = mp;
barrier(CLK_LOCAL_MEM_FENCE);
if (tid < GROUP_SIZE/WARP_SIZE)
{
for (int offset = 1; offset < GROUP_SIZE/WARP_SIZE; offset <<= 1)
warppos[tid] += warppos[select(GROUP_SIZE/WARP_SIZE, tid - offset, tid >= offset)];
}
barrier(CLK_LOCAL_MEM_FENCE);
mp += start + select(0, warppos[tid / WARP_SIZE - 1], tid / WARP_SIZE > 0);
int start32 = start >> 5;
start += mypos[GROUP_SIZE - 1] + warppos[GROUP_SIZE / WARP_SIZE - 2];
//if (tid == GROUP_SIZE - 1 && mypos[tid] > (GROUP_SIZE/2) * 32)
// printf("Oops: %d\n", mypos[tid]);
data[tid] = select(0U, remainder, tid == 0);
barrier(CLK_LOCAL_MEM_FENCE);
if (pstart)
{
int kpos = mp - mylen;
int kpos0 = (kpos >> 5) - start32;
int kpos1 = kpos & 31;
uint kval = (uint)k << (32 - RICE_PARAM_BITS);
uint kval0 = kval >> kpos1;
uint kval1 = kval << (32 - kpos1);
if (kval0) atomic_or(&data[kpos0], kval0);
if (kpos1 && kval1) atomic_or(&data[kpos0 + 1], kval1);
}
if (offs >= task.residualOrder && offs < bs)
{
int qpos = mp - k - 1;
int qpos0 = (qpos >> 5) - start32;
int qpos1 = qpos & 31;
uint qval = (1U << 31) | (v << (31 - k));
uint qval0 = qval >> qpos1;
uint qval1= qval << (32 - qpos1);
if (qval0) atomic_or(&data[qpos0], qval0);
if (qpos1 && qval1) atomic_or(&data[qpos0 + 1], qval1);
}
barrier(CLK_LOCAL_MEM_FENCE);
if ((start32 + tid) * 32 <= start)
output[start32 + tid] = as_int(as_char4(data[tid]).wzyx);
remainder = data[start / 32 - start32];
}
// if (tid == 0 && start != task.encodingOffset - task.headerLen + task.size)
//printf("size mismatch: %d != %d\n", start, task.encodingOffset - task.headerLen + task.size);
#endif
}
#endif /* DO_RICE */
#endif /* DO_PARTITIONS */
#endif
Reference in New Issue View Git Blame Copy Permalink