|
| 1 | +# SPDX-FileCopyrightText: 2022 - 2023 Intel Corporation |
| 2 | +# |
| 3 | +# SPDX-License-Identifier: Apache-2.0 |
| 4 | +"""Numba-dpex implementation for gaussian elimination.""" |
| 5 | + |
| 6 | +import dpctl |
| 7 | +import numba_dpex |
| 8 | + |
| 9 | + |
| 10 | +@numba_dpex.kernel() |
| 11 | +def gaussian_kernel_1(m, a, size, t): |
| 12 | + """Find the multiplier matrix. |
| 13 | +
|
| 14 | + Args: |
| 15 | + m: multiplier matrix. |
| 16 | + a: input matrix. |
| 17 | + size: sizew of matrix. |
| 18 | + t: current iteration. |
| 19 | + """ |
| 20 | + if ( |
| 21 | + numba_dpex.get_local_id(2) |
| 22 | + + numba_dpex.get_group_id(2) * numba_dpex.get_local_size(2) |
| 23 | + >= size - 1 - t |
| 24 | + ): |
| 25 | + return |
| 26 | + |
| 27 | + m[ |
| 28 | + size |
| 29 | + * ( |
| 30 | + numba_dpex.get_local_size(2) * numba_dpex.get_group_id(2) |
| 31 | + + numba_dpex.get_local_id(2) |
| 32 | + + t |
| 33 | + + 1 |
| 34 | + ) |
| 35 | + + t |
| 36 | + ] = ( |
| 37 | + a[ |
| 38 | + size |
| 39 | + * ( |
| 40 | + numba_dpex.get_local_size(2) * numba_dpex.get_group_id(2) |
| 41 | + + numba_dpex.get_local_id(2) |
| 42 | + + t |
| 43 | + + 1 |
| 44 | + ) |
| 45 | + + t |
| 46 | + ] |
| 47 | + / a[size * t + t] |
| 48 | + ) |
| 49 | + |
| 50 | + |
| 51 | +@numba_dpex.kernel() |
| 52 | +def gaussian_kernel_2(m, a, b, size, t): |
| 53 | + """Perform Gaussian elimination using gaussian operations for a iteration. |
| 54 | +
|
| 55 | + Args: |
| 56 | + m: multiplier matrix. |
| 57 | + a: input matrix. |
| 58 | + b: column matrix. |
| 59 | + size: size of matrices. |
| 60 | + t: current iteration. |
| 61 | + """ |
| 62 | + if ( |
| 63 | + numba_dpex.get_local_id(2) |
| 64 | + + numba_dpex.get_group_id(2) * numba_dpex.get_local_size(2) |
| 65 | + >= size - 1 - t |
| 66 | + ): |
| 67 | + return |
| 68 | + |
| 69 | + if ( |
| 70 | + numba_dpex.get_local_id(1) |
| 71 | + + numba_dpex.get_group_id(1) * numba_dpex.get_local_size(1) |
| 72 | + >= size - t |
| 73 | + ): |
| 74 | + return |
| 75 | + |
| 76 | + xidx = numba_dpex.get_group_id(2) * numba_dpex.get_local_size( |
| 77 | + 2 |
| 78 | + ) + numba_dpex.get_local_id(2) |
| 79 | + yidx = numba_dpex.get_group_id(1) * numba_dpex.get_local_size( |
| 80 | + 1 |
| 81 | + ) + numba_dpex.get_local_id(1) |
| 82 | + |
| 83 | + a[size * (xidx + 1 + t) + (yidx + t)] -= ( |
| 84 | + m[size * (xidx + 1 + t) + t] * a[size * t + (yidx + t)] |
| 85 | + ) |
| 86 | + if yidx == 0: |
| 87 | + b[xidx + 1 + t] -= m[size * (xidx + 1 + t) + (yidx + t)] * b[t] |
| 88 | + |
| 89 | + |
| 90 | +def gaussian(a, b, m, size, block_sizeXY, result): |
| 91 | + """Perform Gaussian elimination using gaussian operations. |
| 92 | +
|
| 93 | + Args: |
| 94 | + a: input matrix. |
| 95 | + b: column matrix. |
| 96 | + m: multiplier matrix. |
| 97 | + size: size of matrices. |
| 98 | + block_sizeXY: grid size. |
| 99 | + result: result matrix. |
| 100 | + """ |
| 101 | + device = dpctl.SyclDevice() |
| 102 | + block_size = device.max_work_group_size |
| 103 | + grid_size = int((size / block_size) + 0 if not (size % block_size) else 1) |
| 104 | + |
| 105 | + blocksize2d = block_sizeXY |
| 106 | + gridsize2d = int( |
| 107 | + (size / blocksize2d) + (0 if not (size % blocksize2d) else 1) |
| 108 | + ) |
| 109 | + |
| 110 | + global_range = numba_dpex.Range(1, 1, grid_size * block_size) |
| 111 | + local_range = numba_dpex.Range(1, 1, block_size) |
| 112 | + |
| 113 | + dim_blockXY = numba_dpex.Range(1, blocksize2d, blocksize2d) |
| 114 | + dim_gridXY = numba_dpex.Range( |
| 115 | + 1, gridsize2d * blocksize2d, gridsize2d * blocksize2d |
| 116 | + ) |
| 117 | + |
| 118 | + for t in range(size - 1): |
| 119 | + gaussian_kernel_1[numba_dpex.NdRange(global_range, local_range)]( |
| 120 | + m, a, size, t |
| 121 | + ) |
| 122 | + |
| 123 | + gaussian_kernel_2[numba_dpex.NdRange(dim_gridXY, dim_blockXY)]( |
| 124 | + m, a, b, size, t |
| 125 | + ) |
| 126 | + |
| 127 | + for i in range(size): |
| 128 | + result[size - i - 1] = b[size - i - 1] |
| 129 | + for j in range(i): |
| 130 | + result[size - i - 1] -= ( |
| 131 | + a[size * (size - i - 1) + (size - j - 1)] * result[size - j - 1] |
| 132 | + ) |
| 133 | + result[size - i - 1] = ( |
| 134 | + result[size - i - 1] / a[size * (size - i - 1) + (size - i - 1)] |
| 135 | + ) |
0 commit comments