<div dir="ltr"> <div class="gmail_extra"> <div class="gmail_quote">On Tue, Apr 12, 2016 at 1:05 AM, Samuel Iglesias Gonsálvez <<a href="mailto:siglesias@igalia.com" target="_blank">siglesias@igalia.com</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">From: Connor Abbott <<a href="mailto:connor.w.abbott@intel.com" target="_blank">connor.w.abbott@intel.com</a>> v2: Move to compiler/nir (Iago) Signed-off-by: Iago Toral Quiroga <<a href="mailto:itoral@igalia.com" target="_blank">itoral@igalia.com</a>> --- src/compiler/Makefile.sources | 1 + src/compiler/nir/nir.h | 7 + src/compiler/nir/nir_lower_double_ops.c | 387 ++++++++++++++++++++++++++++++++ 3 files changed, 395 insertions(+) create mode 100644 src/compiler/nir/nir_lower_double_ops.c diff --git a/src/compiler/Makefile.sources b/src/compiler/Makefile.sources index 6f09abf..db7ca3b 100644 --- a/src/compiler/Makefile.sources +++ b/src/compiler/Makefile.sources @@ -187,6 +187,7 @@ NIR_FILES = \ nir/nir_lower_alu_to_scalar.c \ nir/nir_lower_atomics.c \ nir/nir_lower_clip.c \ + nir/nir_lower_double_ops.c \ nir/nir_lower_double_packing.c \ nir/nir_lower_global_vars_to_local.c \ nir/nir_lower_gs_intrinsics.c \ diff --git a/src/compiler/nir/nir.h b/src/compiler/nir/nir.h index ebac750..434d92b 100644 --- a/src/compiler/nir/nir.h +++ b/src/compiler/nir/nir.h @@ -2282,6 +2282,13 @@ void nir_lower_to_source_mods(nir_shader *shader); bool nir_lower_gs_intrinsics(nir_shader *shader); +typedef enum { + nir_lower_drcp = (1 << 0), + nir_lower_dsqrt = (1 << 1), + nir_lower_drsq = (1 << 2), +} nir_lower_doubles_options; + +void nir_lower_doubles(nir_shader *shader, nir_lower_doubles_options options); void nir_lower_double_pack(nir_shader *shader); bool nir_normalize_cubemap_coords(nir_shader *shader); diff --git a/src/compiler/nir/nir_lower_double_ops.c b/src/compiler/nir/nir_lower_double_ops.c new file mode 100644 index 0000000..4cd153c --- /dev/null +++ b/src/compiler/nir/nir_lower_double_ops.c @@ -0,0 +1,387 @@ +/* + * Copyright © 2015 Intel Corporation + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice (including the next + * paragraph) shall be included in all copies or substantial portions of the + * Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. + * + */ + +#include "nir.h" +#include "nir_builder.h" +#include "c99_math.h" + +/* + * Lowers some unsupported double operations, using only: + * + * - pack/unpackDouble2x32 + * - conversion to/from single-precision + * - double add, mul, and fma + * - conditional select + * - 32-bit integer and floating point arithmetic + */ + +/* Creates a double with the exponent bits set to a given integer value */ +static nir_ssa_def * +set_exponent(nir_builder *b, nir_ssa_def *src, nir_ssa_def *exp) +{ + /* Split into bits 0-31 and 32-63 */ + nir_ssa_def *lo = nir_unpack_double_2x32_split_x(b, src); + nir_ssa_def *hi = nir_unpack_double_2x32_split_y(b, src); + + /* The exponent is bits 52-62, or 20-30 of the high word, so set those bits + * to 1023 </blockquote><div> </div><div>We're not setting them to 1023, we're setting it to exp. </div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> + */ + nir_ssa_def *new_hi = nir_bfi(b, nir_imm_uint(b, 0x7ff00000), + exp, hi); </blockquote><div> </div><div>Does this really need a line-wrap? It looks shorter than the comment above. </div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> + /* recombine */ + return nir_pack_double_2x32_split(b, lo, new_hi); +} + +static nir_ssa_def * +get_exponent(nir_builder *b, nir_ssa_def *src) +{ + /* get bits 32-63 */ + nir_ssa_def *hi = nir_unpack_double_2x32_split_y(b, src); + + /* extract bits 20-30 of the high word */ + return nir_ubitfield_extract(b, hi, nir_imm_int(b, 20), nir_imm_int(b, 11)); +} + +/* Return infinity with the sign of the given source which is +/-0 */ + +static nir_ssa_def * +get_signed_inf(nir_builder *b, nir_ssa_def *zero) +{ + nir_ssa_def *zero_split = nir_unpack_double_2x32(b, zero); + nir_ssa_def *zero_hi = nir_swizzle(b, zero_split, (unsigned[]) {1}, 1, false); + + /* The bit pattern for infinity is 0x7ff0000000000000, where the sign bit + * is the highest bit. Only the sign bit can be non-zero in the passed in + * source. So we essentially need to OR the infinity and the zero, except + * the low 32 bits are always 0 so we can construct the correct high 32 + * bits and then pack it together with zero low 32 bits. + */ + nir_ssa_def *inf_hi = nir_ior(b, nir_imm_uint(b, 0x7ff00000), zero_hi); + nir_ssa_def *inf_split = nir_vec2(b, nir_imm_int(b, 0), inf_hi); + return nir_pack_double_2x32(b, inf_split); </blockquote><div> </div><div>Just make this pack_double(b, nir_vec2(b, )). No need for the temporary. </div><div>Other than that </div><div>Reviewed-by: Jason Ekstrand <<a href="mailto:jason@jlekstrand.net">jason@jlekstrand.net</a>> </div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> +} + +/* + * Generates the correctly-signed infinity if the source was zero, and flushes + * the result to 0 if the source was infinity or the calculated exponent was + * too small to be representable. + */ + +static nir_ssa_def * +fix_inv_result(nir_builder *b, nir_ssa_def *res, nir_ssa_def *src, + nir_ssa_def *exp) +{ + /* If the exponent is too small or the original input was infinity/NaN, + * force the result to 0 (flush denorms) to avoid the work of handling + * denorms properly. Note that this doesn't preserve positive/negative + * zeros, but GLSL doesn't require it. + */ + res = nir_bcsel(b, nir_ior(b, nir_ige(b, nir_imm_int(b, 0), exp), + nir_feq(b, nir_fabs(b, src), + nir_imm_double(b, INFINITY))), + nir_imm_double(b, 0.0f), res); + + /* If the original input was 0, generate the correctly-signed infinity */ + res = nir_bcsel(b, nir_fne(b, src, nir_imm_double(b, 0.0f)), + res, get_signed_inf(b, src)); </blockquote><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> + + return res; + +} + +static nir_ssa_def * +lower_rcp(nir_builder *b, nir_ssa_def *src) +{ + /* normalize the input to avoid range issues */ + nir_ssa_def *src_norm = set_exponent(b, src, nir_imm_int(b, 1023)); + + /* cast to float, do an rcp, and then cast back to get an approximate + * result + */ + nir_ssa_def *ra = nir_f2d(b, nir_frcp(b, nir_d2f(b, src_norm))); + + /* Fixup the exponent of the result - note that we check if this is too + * small below. + */ + nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra), + nir_isub(b, get_exponent(b, src), + nir_imm_int(b, 1023))); + + ra = set_exponent(b, ra, new_exp); + + /* Do a few Newton-Raphson steps to improve precision. + * + * Each step doubles the precision, and we started off with around 24 bits, + * so we only need to do 2 steps to get to full precision. The step is: + * + * x_new = x * (2 - x*src) + * + * But we can re-arrange this to improve precision by using another fused + * multiply-add: + * + * x_new = x + x * (1 - x*src) + * + * See <a href="https://en.wikipedia.org/wiki/Division_algorithm" rel="noreferrer" target="_blank">https://en.wikipedia.org/wiki/Division_algorithm</a> for more details. + */ + + ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra); + ra = nir_ffma(b, ra, nir_ffma(b, ra, src, nir_imm_double(b, -1)), ra); + + return fix_inv_result(b, ra, src, new_exp); +} + +static nir_ssa_def * +lower_sqrt_rsq(nir_builder *b, nir_ssa_def *src, bool sqrt) +{ + /* We want to compute: + * + * 1/sqrt(m * 2^e) + * + * When the exponent is even, this is equivalent to: + * + * 1/sqrt(m) * 2^(-e/2) + * + * and then the exponent is odd, this is equal to: + * + * 1/sqrt(m * 2) * 2^(-(e - 1)/2) + * + * where the m * 2 is absorbed into the exponent. So we want the exponent + * inside the square root to be 1 if e is odd and 0 if e is even, and we + * want to subtract off e/2 from the final exponent, rounded to negative + * infinity. We can do the former by first computing the unbiased exponent, + * and then AND'ing it with 1 to get 0 or 1, and we can do the latter by + * shifting right by 1. + */ + + nir_ssa_def *unbiased_exp = nir_isub(b, get_exponent(b, src), + nir_imm_int(b, 1023)); + nir_ssa_def *even = nir_iand(b, unbiased_exp, nir_imm_int(b, 1)); + nir_ssa_def *half = nir_ishr(b, unbiased_exp, nir_imm_int(b, 1)); </blockquote><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> + + nir_ssa_def *src_norm = set_exponent(b, src, + nir_iadd(b, nir_imm_int(b, 1023), + even)); + + nir_ssa_def *ra = nir_f2d(b, nir_frsq(b, nir_d2f(b, src_norm))); + nir_ssa_def *new_exp = nir_isub(b, get_exponent(b, ra), half); + ra = set_exponent(b, ra, new_exp); + + /* + * The following implements an iterative algorithm that's very similar + * between sqrt and rsqrt. We start with an iteration of Goldschmit's + * algorithm, which looks like: + * + * a = the source + * y_0 = initial (single-precision) rsqrt estimate + * + * h_0 = .5 * y_0 + * g_0 = a * y_0 + * r_0 = .5 - h_0 * g_0 + * g_1 = g_0 * r_0 + g_0 + * h_1 = h_0 * r_0 + h_0 + * + * Now g_1 ~= sqrt(a), and h_1 ~= 1/(2 * sqrt(a)). We could continue + * applying another round of Goldschmit, but since we would never refer + * back to a (the original source), we would add too much rounding error. + * So instead, we do one last round of Newton-Raphson, which has better + * rounding characteristics, to get the final rounding correct. This is + * split into two cases: + * + * 1. sqrt + * + * Normally, doing a round of Newton-Raphson for sqrt involves taking a + * reciprocal of the original estimate, which is slow since it isn't + * supported in HW. But we can take advantage of the fact that we already + * computed a good estimate of 1/(2 * g_1) by rearranging it like so: + * + * g_2 = .5 * (g_1 + a / g_1) + * = g_1 + .5 * (a / g_1 - g_1) + * = g_1 + (.5 / g_1) * (a - g_1^2) + * = g_1 + h_1 * (a - g_1^2) + * + * The second term represents the error, and by splitting it out we can get + * better precision by computing it as part of a fused multiply-add. Since + * both Newton-Raphson and Goldschmit approximately double the precision of + * the result, these two steps should be enough. + * + * 2. rsqrt + * + * First off, note that the first round of the Goldschmit algorithm is + * really just a Newton-Raphson step in disguise: + * + * h_1 = h_0 * (.5 - h_0 * g_0) + h_0 + * = h_0 * (1.5 - h_0 * g_0) + * = h_0 * (1.5 - .5 * a * y_0^2) + * = (.5 * y_0) * (1.5 - .5 * a * y_0^2) + * + * which is the standard formula multiplied by .5. Unlike in the sqrt case, + * we don't need the inverse to do a Newton-Raphson step; we just need h_1, + * so we can skip the calculation of g_1. Instead, we simply do another + * Newton-Raphson step: + * + * y_1 = 2 * h_1 + * r_1 = .5 - h_1 * y_1 * a + * y_2 = y_1 * r_1 + y_1 + * + * Where the difference from Goldschmit is that we calculate y_1 * a + * instead of using g_1. Doing it this way should be as fast as computing + * y_1 up front instead of h_1, and it lets us share the code for the + * initial Goldschmit step with the sqrt case. + * + * Putting it together, the computations are: + * + * h_0 = .5 * y_0 + * g_0 = a * y_0 + * r_0 = .5 - h_0 * g_0 + * h_1 = h_0 * r_0 + h_0 + * if sqrt: + * g_1 = g_0 * r_0 + g_0 + * r_1 = a - g_1 * g_1 + * g_2 = h_1 * r_1 + g_1 + * else: + * y_1 = 2 * h_1 + * r_1 = .5 - y_1 * (h_1 * a) + * y_2 = y_1 * r_1 + y_1 + * + * For more on the ideas behind this, see "Software Division and Square + * Root Using Goldschmit's Algorithms" by Markstein and the Wikipedia page + * on square roots + * (<a href="https://en.wikipedia.org/wiki/Methods_of_computing_square_roots" rel="noreferrer" target="_blank">https://en.wikipedia.org/wiki/Methods_of_computing_square_roots</a>). + */ + + nir_ssa_def *one_half = nir_imm_double(b, 0.5); + nir_ssa_def *h_0 = nir_fmul(b, one_half, ra); + nir_ssa_def *g_0 = nir_fmul(b, src, ra); + nir_ssa_def *r_0 = nir_ffma(b, nir_fneg(b, h_0), g_0, one_half); + nir_ssa_def *h_1 = nir_ffma(b, h_0, r_0, h_0); + nir_ssa_def *res; + if (sqrt) { + nir_ssa_def *g_1 = nir_ffma(b, g_0, r_0, g_0); + nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, g_1), g_1, src); + res = nir_ffma(b, h_1, r_1, g_1); + } else { + nir_ssa_def *y_1 = nir_fmul(b, nir_imm_double(b, 2.0), h_1); + nir_ssa_def *r_1 = nir_ffma(b, nir_fneg(b, y_1), nir_fmul(b, h_1, src), + one_half); + res = nir_ffma(b, y_1, r_1, y_1); + } + + if (sqrt) { + /* Here, the special cases we need to handle are + * 0 -> 0 and + * +inf -> +inf + */ + res = nir_bcsel(b, nir_ior(b, nir_feq(b, src, nir_imm_double(b, 0.0)), + nir_feq(b, src, nir_imm_double(b, INFINITY))), + src, res); + } else { + res = fix_inv_result(b, res, src, new_exp); + } + + return res; +} + +static void +lower_doubles_instr(nir_alu_instr *instr, nir_lower_doubles_options options) +{ + assert(instr->dest.dest.is_ssa); + if (instr->dest.dest.ssa.bit_size != 64) + return; + + switch (instr->op) { + case nir_op_frcp: + if (!(options & nir_lower_drcp)) + return; + break; + + case nir_op_fsqrt: + if (!(options & nir_lower_dsqrt)) + return; + break; + + case nir_op_frsq: + if (!(options & nir_lower_drsq)) + return; + break; + + default: + return; + } + + nir_builder bld; + nir_builder_init(&bld, nir_cf_node_get_function(&instr->instr.block->cf_node)); + bld.cursor = nir_before_instr(&instr->instr); + + nir_ssa_def *src = nir_fmov_alu(&bld, instr->src[0], + instr->dest.dest.ssa.num_components); + + nir_ssa_def *result; + + switch (instr->op) { + case nir_op_frcp: + result = lower_rcp(&bld, src); + break; + case nir_op_fsqrt: + result = lower_sqrt_rsq(&bld, src, true); + break; + case nir_op_frsq: + result = lower_sqrt_rsq(&bld, src, false); + break; + default: + unreachable("unhandled opcode"); + } + + nir_ssa_def_rewrite_uses(&instr->dest.dest.ssa, nir_src_for_ssa(result)); + nir_instr_remove(&instr->instr); +} + +static bool +lower_doubles_block(nir_block *block, void *ctx) +{ + nir_lower_doubles_options options = *((nir_lower_doubles_options *) ctx); + + nir_foreach_instr_safe(block, instr) { + if (instr->type != nir_instr_type_alu) + continue; + + lower_doubles_instr(nir_instr_as_alu(instr), options); + } + + return true; +} + +static void +lower_doubles_impl(nir_function_impl *impl, nir_lower_doubles_options options) +{ + nir_foreach_block(impl, lower_doubles_block, &options); +} + +void +nir_lower_doubles(nir_shader *shader, nir_lower_doubles_options options) +{ + nir_foreach_function(shader, function) { + if (function->impl) + lower_doubles_impl(function->impl, options); + } +} -- 2.5.0 _______________________________________________ mesa-dev mailing list <a href="mailto:mesa-dev@lists.freedesktop.org" target="_blank">mesa-dev@lists.freedesktop.org</a> <a href="https://lists.freedesktop.org/mailman/listinfo/mesa-dev" rel="noreferrer" target="_blank">https://lists.freedesktop.org/mailman/listinfo/mesa-dev</a> </blockquote></div> </div></div>