diff --git a/arrayfire/blas.py b/arrayfire/blas.py index f0e9dfdc6..1f1cb1359 100644 --- a/arrayfire/blas.py +++ b/arrayfire/blas.py @@ -202,3 +202,109 @@ def dot(lhs, rhs, lhs_opts=MATPROP.NONE, rhs_opts=MATPROP.NONE, return_scalar = safe_call(backend.get().af_dot(c_pointer(out.arr), lhs.arr, rhs.arr, lhs_opts.value, rhs_opts.value)) return out + +def gemm(lhs, rhs, alpha=1.0, beta=0.0, lhs_opts=MATPROP.NONE, rhs_opts=MATPROP.NONE, C=None): + """ + BLAS general matrix multiply (GEMM) of two af_array objects. + + This provides a general interface to the BLAS level 3 general matrix multiply (GEMM), which is generally defined as: + + C = α ∗ opA(A) opB(B)+ β∗C + + where α (alpha) and β (beta) are both scalars; A and B are the matrix multiply operands; + and opA and opB are noop (if AF_MAT_NONE) or transpose (if AF_MAT_TRANS) operations + on A or B before the actual GEMM operation. + Batched GEMM is supported if at least either A or B have more than two dimensions + (see af::matmul for more details on broadcasting). + However, only one alpha and one beta can be used for all of the batched matrix operands. + + Parameters + ---------- + + lhs : af.Array + A 2 dimensional, real or complex arrayfire array. + + rhs : af.Array + A 2 dimensional, real or complex arrayfire array. + + alpha : scalar + + beta : scalar + + lhs_opts: optional: af.MATPROP. default: af.MATPROP.NONE. + Can be one of + - af.MATPROP.NONE - If no op should be done on `lhs`. + - af.MATPROP.TRANS - If `lhs` has to be transposed before multiplying. + - af.MATPROP.CTRANS - If `lhs` has to be hermitian transposed before multiplying. + + rhs_opts: optional: af.MATPROP. default: af.MATPROP.NONE. + Can be one of + - af.MATPROP.NONE - If no op should be done on `rhs`. + - af.MATPROP.TRANS - If `rhs` has to be transposed before multiplying. + - af.MATPROP.CTRANS - If `rhs` has to be hermitian transposed before multiplying. + + Returns + ------- + + out : af.Array + Output of the matrix multiplication on `lhs` and `rhs`. + + Note + ----- + + - The data types of `lhs` and `rhs` should be the same. + - Batches are not supported. + + """ + if C is None: + out = Array() + else: + out = C + + ltype = lhs.dtype() + + if ltype == Dtype.f32: + aptr = c_cast(c_pointer(c_float_t(alpha)),c_void_ptr_t) + bptr = c_cast(c_pointer(c_float_t(beta)), c_void_ptr_t) + elif ltype == Dtype.c32: + if isinstance(alpha, af_cfloat_t): + aptr = c_cast(c_pointer(alpha), c_void_ptr_t) + elif isinstance(alpha, tuple): + aptr = c_cast(c_pointer(af_cfloat_t(alpha[0], alpha[1])), c_void_ptr_t) + else: + aptr = c_cast(c_pointer(af_cfloat_t(alpha)), c_void_ptr_t) + + if isinstance(beta, af_cfloat_t): + bptr = c_cast(c_pointer(beta), c_void_ptr_t) + elif isinstance(beta, tuple): + bptr = c_cast(c_pointer(af_cfloat_t(beta[0], beta[1])), c_void_ptr_t) + else: + bptr = c_cast(c_pointer(af_cfloat_t(beta)), c_void_ptr_t) + + elif ltype == Dtype.f64: + aptr = c_cast(c_pointer(c_double_t(alpha)),c_void_ptr_t) + bptr = c_cast(c_pointer(c_double_t(beta)), c_void_ptr_t) + elif ltype == Dtype.c64: + if isinstance(alpha, af_cdouble_t): + aptr = c_cast(c_pointer(alpha), c_void_ptr_t) + elif isinstance(alpha, tuple): + aptr = c_cast(c_pointer(af_cdouble_t(alpha[0], alpha[1])), c_void_ptr_t) + else: + aptr = c_cast(c_pointer(af_cdouble_t(alpha)), c_void_ptr_t) + + if isinstance(beta, af_cdouble_t): + bptr = c_cast(c_pointer(beta), c_void_ptr_t) + elif isinstance(beta, tuple): + bptr = c_cast(c_pointer(af_cdouble_t(beta[0], beta[1])), c_void_ptr_t) + else: + bptr = c_cast(c_pointer(af_cdouble_t(beta)), c_void_ptr_t) + elif ltype == Dtype.f16: + raise TypeError("fp16 currently unsupported gemm() input type") + else: + raise TypeError("unsupported input type") + + + safe_call(backend.get().af_gemm(c_pointer(out.arr), + lhs_opts.value, rhs_opts.value, + aptr, lhs.arr, rhs.arr, bptr)) + return out diff --git a/arrayfire/library.py b/arrayfire/library.py index 863c6746f..009a8e122 100644 --- a/arrayfire/library.py +++ b/arrayfire/library.py @@ -31,6 +31,13 @@ c_void_ptr_t = ct.c_void_p c_char_ptr_t = ct.c_char_p c_size_t = ct.c_size_t +c_cast = ct.cast + +class af_cfloat_t(ct.Structure): + _fields_ = [("real", ct.c_float), ("imag", ct.c_float)] + +class af_cdouble_t(ct.Structure): + _fields_ = [("real", ct.c_double), ("imag", ct.c_double)] AF_VER_MAJOR = '3' diff --git a/arrayfire/signal.py b/arrayfire/signal.py index 1fa6c424d..35e0fba87 100644 --- a/arrayfire/signal.py +++ b/arrayfire/signal.py @@ -27,7 +27,7 @@ def _scale_pos_axis1(y_curr, y_orig): dy = y_orig[0, 1, 0, 0] - y0 return((y_curr - y0) / dy) -def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None): +def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None, output = None): """ Interpolate along a single dimension.Interpolation is performed along axis 0 of the input array. @@ -51,6 +51,10 @@ def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None): xp : af.Array The x-coordinates of the input data points + output: None or af.Array + Optional preallocated output array. If it is a sub-array of an existing af_array, + only the corresponding portion of the af_array will be overwritten + Returns ------- @@ -65,20 +69,86 @@ def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None): where N is the length of the first dimension of `signal`. """ - output = Array() + if output is None: + output = Array() + + if(xp is not None): + pos0 = _scale_pos_axis0(x, xp) + else: + pos0 = x + + safe_call(backend.get().af_approx1(c_pointer(output.arr), signal.arr, pos0.arr, + method.value, c_float_t(off_grid))) - if(xp is not None): - pos0 = _scale_pos_axis0(x, xp) else: - pos0 = x + if(xp is not None): + pos0 = _scale_pos_axis0(x, xp) + else: + pos0 = x + safe_call(backend.get().af_approx1_v2(c_pointer(output.arr), signal.arr, pos0.arr, + method.value, c_float_t(off_grid))) + return output + + +def approx1_uniform(signal, x, interp_dim, idx_start, idx_step, method=INTERP.LINEAR, off_grid=0.0, output = None): + """ + Interpolation on one dimensional signals along specified dimension. + + af_approx1_uniform() accepts the dimension to perform the interpolation along the input. + It also accepts start and step values which define the uniform range of corresponding indices. + + Parameters + ---------- + + signal: af.Array + Input signal array (signal = f(x)) + + x: af.Array + The x-coordinates of the interpolation points. The interpolation + function is queried at these set of points. + + interp_dim: scalar + is the dimension to perform interpolation across. + + idx_start: scalar + is the first index value along interp_dim. + + idx_step: scalar + is the uniform spacing value between subsequent indices along interp_dim. - safe_call(backend.get().af_approx1(c_pointer(output.arr), signal.arr, pos0.arr, - method.value, c_float_t(off_grid))) + method: optional: af.INTERP. default: af.INTERP.LINEAR. + Interpolation method. + + off_grid: optional: scalar. default: 0.0. + The value used for positions outside the range. + + output: None or af.Array + Optional preallocated output array. If it is a sub-array of an existing af_array, + only the corresponding portion of the af_array will be overwritten + + Returns + ------- + + output: af.Array + Values calculated at interpolation points. + + """ + + if output is None: + output = Array() + + safe_call(backend.get().af_approx1_uniform(c_pointer(output.arr), signal.arr, x.arr, + c_dim_t(interp_dim), c_double_t(idx_start), c_double_t(idx_step), + method.value, c_float_t(off_grid))) + else: + safe_call(backend.get().af_approx1_uniform_v2(c_pointer(output.arr), signal.arr, x.arr, + c_dim_t(interp_dim), c_double_t(idx_start), c_double_t(idx_step), + method.value, c_float_t(off_grid))) return output + def approx2(signal, x, y, - method=INTERP.LINEAR, off_grid=0.0, xp = None, yp = None - ): + method=INTERP.LINEAR, off_grid=0.0, xp = None, yp = None, output = None): """ Interpolate along a two dimension.Interpolation is performed along axes 0 and 1 of the input array. @@ -112,6 +182,10 @@ def approx2(signal, x, y, The y-coordinates of the input data points. The convention followed is that the y-coordinates vary along axis 1 + output: None or af.Array + Optional preallocated output array. If it is a sub-array of an existing af_array, + only the corresponding portion of the af_array will be overwritten + Returns ------- @@ -127,22 +201,114 @@ def approx2(signal, x, y, and N is the length of the second dimension of `signal`. """ - output = Array() - - if(xp is not None): - pos0 = _scale_pos_axis0(x, xp) + if output is None: + output = Array() + + if(xp is not None): + pos0 = _scale_pos_axis0(x, xp) + else: + pos0 = x + + if(yp is not None): + pos1 = _scale_pos_axis1(y, yp) + else: + pos1 = y + + safe_call(backend.get().af_approx2(c_pointer(output.arr), signal.arr, + pos0.arr, pos1.arr, method.value, c_float_t(off_grid))) else: - pos0 = x + if(xp is not None): + pos0 = _scale_pos_axis0(x, xp) + else: + pos0 = x + + if(yp is not None): + pos1 = _scale_pos_axis1(y, yp) + else: + pos1 = y + + safe_call(backend.get().af_approx2_v2(c_pointer(output.arr), signal.arr, + pos0.arr, pos1.arr, method.value, c_float_t(off_grid))) + + return output + +def approx2_uniform(signal, pos0, interp_dim0, idx_start0, idx_step0, pos1, interp_dim1, idx_start1, idx_step1, + method=INTERP.LINEAR, off_grid=0.0, output = None): + """ + Interpolate along two uniformly spaced dimensions of the input array. + af_approx2_uniform() accepts two dimensions to perform the interpolation along the input. + It also accepts start and step values which define the uniform range of corresponding indices. + + Parameters + ---------- + + signal: af.Array + Input signal array (signal = f(x, y)) + + pos0 : af.Array + positions of the interpolation points along interp_dim0. + + interp_dim0: scalar + is the first dimension to perform interpolation across. + + idx_start0: scalar + is the first index value along interp_dim0. + + idx_step0: scalar + is the uniform spacing value between subsequent indices along interp_dim0. + + pos1 : af.Array + positions of the interpolation points along interp_dim1. + + interp_dim1: scalar + is the second dimension to perform interpolation across. + + idx_start1: scalar + is the first index value along interp_dim1. + + idx_step1: scalar + is the uniform spacing value between subsequent indices along interp_dim1. + + method: optional: af.INTERP. default: af.INTERP.LINEAR. + Interpolation method. + + off_grid: optional: scalar. default: 0.0. + The value used for positions outside the range. + + output: None or af.Array + Optional preallocated output array. If it is a sub-array of an existing af_array, + only the corresponding portion of the af_array will be overwritten + + Returns + ------- + + output: af.Array + Values calculated at interpolation points. + + Note + ----- + This holds applicable when x_input/y_input isn't provided: - if(yp is not None): - pos1 = _scale_pos_axis1(y, yp) + The initial measurements are assumed to have taken place at equal steps between [(0,0) - [M - 1, N - 1]] + where M is the length of the first dimension of `signal`, + and N is the length of the second dimension of `signal`. + """ + + if output is None: + output = Array() + safe_call(backend.get().af_approx2_uniform(c_pointer(output.arr), signal.arr, + pos0.arr, c_dim_t(interp_dim0), c_double_t(idx_start0), c_double_t(idx_step0), + pos1.arr, c_dim_t(interp_dim1), c_double_t(idx_start1), c_double_t(idx_step1), + method.value, c_float_t(off_grid))) else: - pos1 = y + safe_call(backend.get().af_approx2_uniform_v2(c_pointer(output.arr), signal.arr, + pos0.arr, c_dim_t(interp_dim0), c_double_t(idx_start0), c_double_t(idx_step0), + pos1.arr, c_dim_t(interp_dim1), c_double_t(idx_start1), c_double_t(idx_step1), + method.value, c_float_t(off_grid))) - safe_call(backend.get().af_approx2(c_pointer(output.arr), signal.arr, - pos0.arr, pos1.arr, method.value, c_float_t(off_grid))) return output + def fft(signal, dim0 = None , scale = None): """ Fast Fourier Transform: 1D