diff --git a/docs/articles_en/documentation/openvino-ir-format/operation-sets/available-opsets/opset14.rst b/docs/articles_en/documentation/openvino-ir-format/operation-sets/available-opsets/opset14.rst index f5b6f025673..b06528d3b27 100644 --- a/docs/articles_en/documentation/openvino-ir-format/operation-sets/available-opsets/opset14.rst +++ b/docs/articles_en/documentation/openvino-ir-format/operation-sets/available-opsets/opset14.rst @@ -103,7 +103,7 @@ Table of Contents * :doc:`I420toRGB <../operation-specs/image/i420-to-rgb-8>` * :doc:`If <../operation-specs/condition/if-8>` * :doc:`Interpolate <../operation-specs/image/interpolate-11>` -* :doc:`Inverse <../operation-specs/matrix/Inverse_14>` +* :doc:`Inverse <../operation-specs/matrix/inverse-14>` * :doc:`IRDFT <../operation-specs/signals/irdft-9>` * :doc:`IsInf <../operation-specs/comparison/isinf-10>` * :doc:`IsNaN <../operation-specs/comparison/isnan-10>` diff --git a/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/Inverse_14.rst b/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/Inverse_14.rst deleted file mode 100644 index 60812117ff5..00000000000 --- a/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/Inverse_14.rst +++ /dev/null @@ -1,135 +0,0 @@ -.. {#openvino_docs_ops_matrix_Inverse_14} - -Inverse -======= - - -.. meta:: - :description: Learn about Inverse-14 - a matrix operation that computes the inverse or adjoint for one matrix or a batch of input matrice. - -**Versioned name**: *Inverse-14* - -**Category**: *Matrix* - -**Short description**: *Inverse* operation computes either the inverse or adjoint (conjugate transposes) of one or a batch of square invertible matrices. - -**Detailed description**: *Inverse* operation computes the inverse of a square matrix. The operation uses LU decomposition with partial pivoting to compute the inverses. (`LU decomposotion with partial pivoting algorithm `__) - -The inverse matrix A^(-1) of a square matrix A is defined as: - -.. math:: - - A \cdot A^{-1} = A^{-1} \cdot A = I - -where **I** is the n-dimensional identity matrix. - -The inverse matrix exists if and only if the input matrix is invertible. In that case, the inverse is unique. If the matrix is not invertible however, the operation may raise an exception or return an undefined result. - -This operation can be used to compute the adjugate matrix instead of the inverse. - -The adjugate matrix adj(A) of a square matrix A is defined as: - -.. math:: - - adj(A) = det(A) \cdot A^{-1} - -where **A^{-1}** is the matrix inverse of A, and **det(A)** is the determinant of A. - -The adjugate matrix exists if and only if the inverse matrix exists. - -**Attribute**: - -* *adjoint* - - * **Description**: Modifies the return value of the operation. If true, the operation returns the adjoint (conjugate transpose) of the input matrices instead of finding the inverse. - * **Range of values**: `true`, `false` - - * ``true`` - output adjugate matrix. - * ``false`` - output inverse matrix. - - * **Type**: `bool` - * **Default value**: `false` - * **Required**: *No* - -**Input**: - -* **1**: `input` - A tensor of shape [B1, B2, ..., Bn, ROW, COL] and type `T` representing the input square matrices. **The inner-most 2 dimensions form square matrices and must be of the same size.** B1, B2, ..., Bn represent any amount of batch dimensions (can be 0 for a single matrix input). **Required.** - -**Output**: - -* **1**: `output` - A tensor with the same type `T` as the input and same shape [B1, B2, ..., Bn, ROW, COL] as the input, representing the inverse matrices (or adjugate matrices) of the input matrices. - -**Types** - -* **T**: any supported floating-point type. - -*Example 1: 2D input matrix.* - -.. code-block:: xml - :force: - - - - - - 3 - 3 - - - - - 3 - 3 - - - - -*Example 2: 3D input tensor with one batch dimension and adjoint=true.* - -.. code-block:: xml - :force: - - - - - - 2 - 4 - 4 - - - - - 2 - 4 - 4 - - - - -*Example 3: 5D input tensor with three batch dimensions.* - -.. code-block:: xml - :force: - - - - - - 5 - 4 - 3 - 2 - 2 - - - - - 5 - 4 - 3 - 2 - 2 - - - diff --git a/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/inverse-14.rst b/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/inverse-14.rst new file mode 100644 index 00000000000..07c173abc6a --- /dev/null +++ b/docs/articles_en/documentation/openvino-ir-format/operation-sets/operation-specs/matrix/inverse-14.rst @@ -0,0 +1,200 @@ +.. {#openvino_docs_ops_matrix_Inverse_14} + +Inverse +======= + + +.. meta:: + :description: Learn about Inverse-14 - a matrix operation that computes the inverse or adjoint for one matrix or a batch of input matrice. + +**Versioned name**: *Inverse-14* + +**Category**: *Matrix* + +**Short description**: *Inverse* operation computes either the inverse or adjoint (conjugate transposes) of one or a batch of square invertible matrices. + +**Detailed description**: *Inverse* operation computes the inverse of a square matrix. + +The inverse matrix :math:`A^{-1}` of a square n-dimensional matrix :math:`A` is defined as: + +.. math:: + + A \cdot A^{-1} = A^{-1} \cdot A = I + +where :math:`I` is the n-dimensional identity matrix. + +The inverse matrix exists if and only if the input matrix is invertible (satisfies any of the properties of the *Invertible Matrix Theorem*). In that case, the inverse exists, and is unique. However, if the matrix is not invertible, the operation may raise an exception or return an undefined result. The operation may return slightly different results for the same input data on different devices due to parallelism and different data types implementations. + +This operation can be used to compute the adjugate matrix instead of the inverse. + +The adjugate matrix :math:`adj(A)` of a square matrix :math:`A` is defined as: + +.. math:: + + adj(A) = det(A) \cdot A^{-1} + +where :math:`A^{-1}` is the matrix inverse of :math:`A`, and :math:`det(A)` is the determinant of :math:`A`. + +The adjugate matrix exists if and only if the inverse matrix exists. + +**Algorithm formulation**: + +This operation uses LU decomposition with partial pivoting to compute the inverse (or adjugate) matrix. + +.. note:: + + LU decomposition decomposes matrix :math:`A` into 2 matrices :math:`L` and :math:`U`, where :math:`L` is the lower triangular matrix with all diagonal elements equal to 1 and :math:`U` is the upper triangular matrix. + +.. math:: + + A = L \cdot U + +.. note:: + + LU decomposition allows to easily obtain determinant of :math:`A`. Notice that since :math:`L`` is a lower triangular matrix with all diagonal elements equal to 1, :math:`det(L) = 1`. + +.. math:: + + det(A) = det(L) * det(U) = 1 * det(U) = det(U) + +.. note:: + + To compute the inverse of :math:`A`, given its LU decomposition, it is enough to solve a set of n simple linear equations. + A simple linear equation is of the form :math:`Ax=b`, where x and b are vectors, and x is the solution of the equation. The goal is to compute x, such that it is the i-th column of matrix :math:`A^{-1}`. To do so, let b be a vector of zeros, except for i-th spot that has a value of one (Then b is an i-th column of matrix I). + It is easy to notice that the x-vectors create the matrix :math:`A^{-1}`, and the b-vectors create the Identity matrix. + +.. math:: + + A \cdot A^{-1} = I \implies A^{-1} = [x_1 \& x_2 \& ... \& x_n], A \cdot x_i = b_i, b_i \in col_i(I) + +.. note:: + + Using the LU decomposition, the simple linear equation can be replaced by two, even simpler linear equations. + +.. math:: + + A \cdot x = b \iff L \cdot U \cdot x = b \iff L \cdot y = b, U \cdot x = y + +Algorithm pseudocode: + +1. Start with original matrix :math:`A`. If the data type of :math:`A` is not f32, convert them to f32 to avoid accumulating rounding errors. +2. Copy initial matrix into matrix :math:`U`. Initialize matrix :math:`L` to be the Identity matrix (zero matrix with all diagonal elements set to 1). +3. Perform LU decomposition with partial pivoting. + + * Repeat this step for each column in the input matrix. + * Let *c* be the index of the currently processed column. + * Find the index of the row with the highest value in a given column - *pivot*. + * If :math:`pivot \neq c`, swap the *pivot* and *c* row in :math:`L`. Repeat for :math:`U`. Note that this operation flips the sign of the determinant, so this has to be accounted for. + * Perform standard Gaussian elimination. + +4. To obtain the inverse, solve for each column of :math:`A^{-1}` as explained above. + + * Solve linear equation :math:`Ly = b` for y (forward substitution) + * Solve linear equation :math:`Ux = y` for x (backward substitution) + * Set x as the corresponding column of the output inverse matrix :math:`A^{-1}` + +5. If adjoint == true, then it is necessary to multiply :math:`A^{-1}` by its determinant. + + * As explained above, it is enough to compute :math:`det(U)`, since :math:`det(U) = det(A)`. + * :math:`det(U)` is just a multiplication of its diagonal elements. + * Account for each row swap in the LU decomposition step - for every row swap, swap the sign of the determinant. + * Multiply all elements of :math:`A^{-1}` by the determinant to obtain the adjugate matrix. + +6. Return the computed matrix. Convert it back from f32 to its original element type. + +**Attribute**: + +* *adjoint* + + * **Description**: Modifies the return value of the operation. If true, the operation returns the adjoint (conjugate transpose) of the input matrices instead of finding the inverse. + * **Range of values**: `true`, `false` + + * ``true`` - output adjugate matrix. + * ``false`` - output inverse matrix. + + * **Type**: `bool` + * **Default value**: `false` + * **Required**: *No* + +**Input**: + +* **1**: `input` - A tensor of shape [B1, B2, ..., Bn, ROW, COL] and type `T` representing the input square matrices. **The inner-most 2 dimensions form square matrices and must be of the same size.** B1, B2, ..., Bn represent any amount of batch dimensions (can be 0 for a single matrix input). **Required.** + +**Output**: + +* **1**: `output` - A tensor with the same type `T` as the input and same shape [B1, B2, ..., Bn, ROW, COL] as the input, representing the inverse matrices (or adjugate matrices) of the input matrices. + +**Types** + +* **T**: any supported floating-point type. Any type other than f32 will be converted to f32 before executing this op, and then converted back to the original input type to avoid accumulating rounding errors. + +*Example 1: 2D input matrix.* + +.. code-block:: xml + :force: + + + + + + 3 + 3 + + + + + 3 + 3 + + + + +*Example 2: 3D input tensor with one batch dimension and adjoint=true.* + +.. code-block:: xml + :force: + + + + + + 2 + 4 + 4 + + + + + 2 + 4 + 4 + + + + +*Example 3: 5D input tensor with three batch dimensions.* + +.. code-block:: xml + :force: + + + + + + 5 + 4 + 3 + 2 + 2 + + + + + 5 + 4 + 3 + 2 + 2 + + + diff --git a/src/core/reference/include/openvino/reference/inverse.hpp b/src/core/reference/include/openvino/reference/inverse.hpp index b3df294dd82..83a60d0b30e 100644 --- a/src/core/reference/include/openvino/reference/inverse.hpp +++ b/src/core/reference/include/openvino/reference/inverse.hpp @@ -9,107 +9,30 @@ #include #include "openvino/core/shape.hpp" +#include "openvino/reference/convert.hpp" namespace ov { namespace reference { -namespace inverse { namespace internal { -template -void lu_decomposition(const T* input, - std::vector>& L, - std::vector>& U, - std::vector& P, +void lu_decomposition(std::vector& input, + std::vector& L, + std::vector& U, + std::vector& P, bool& sign, size_t b, - size_t n) { - // Make L identity, U a copy of input and P a range(0, n) - const auto batch_idx = b * n * n; - for (size_t i = 0; i < n; ++i) { - P[i] = static_cast(i); - L[i][i] = T{1}; + size_t n, + size_t n_squared); - auto i_idx = i * n; - for (size_t j = 0; j < n; ++j) { - U[i][j] = input[batch_idx + i_idx + j]; - } - } - - for (size_t k = 0; k < n; ++k) { - // Partial Pivoting - auto pivot_row = k; - for (auto i = k + 1; i < n; ++i) { - if (std::abs(U[i][k]) > std::abs(U[pivot_row][k])) { - pivot_row = i; - } - } - - if (pivot_row != k) { - // Swap rows in L, U (A) and P - std::swap(U[k], U[pivot_row]); - std::swap(L[k], L[pivot_row]); - std::swap(P[k], P[pivot_row]); - sign = !sign; - } - - for (auto i = k + 1; i < n; ++i) { - L[i][k] = U[i][k] / U[k][k]; - for (auto j = k; j < n; ++j) { - U[i][j] -= L[i][k] * U[k][j]; - } - } - } -} - -template -void lu_solve(T* output, - std::vector>& L, - std::vector>& U, - std::vector& P, +void lu_solve(std::vector& output, + std::vector& L, + std::vector& U, + std::vector& P, size_t b, size_t n, - size_t column) { - std::vector B(n, T{0}); - std::vector X(n, T{0}); - std::vector Y(n, T{0}); - B[column] = T{1}; + size_t n_squared); - // Forward substitution: Ly = Pb - for (size_t i = 0; i < n; ++i) { - Y[i] = B[P[i]]; - for (size_t j = 0; j < i; ++j) { - Y[i] -= L[i][j] * Y[j]; - } - } - - // Backward substitution: Ux = y - for (int i = static_cast(n - 1); i >= 0; --i) { - X[i] = Y[i]; - for (size_t j = static_cast(i) + 1; j < n; ++j) { - X[i] -= U[i][j] * X[j]; - } - X[i] /= U[i][i]; - } - - size_t batch_idx = b * n * n; - for (size_t row = 0; row < n; ++row) { - output[batch_idx + row * n + column] = X[row]; - } -} - -template -void to_adjoint(T* output, std::vector>& U, bool sign, size_t b, size_t n) { - T determinant = sign ? T{1} : T{-1}; - - for (size_t i = 0; i < n; ++i) { - determinant *= U[i][i]; - } - - const auto batch_idx = b * n * n; - for (auto idx = batch_idx; idx < batch_idx + n * n; ++idx) { - output[idx] *= determinant; - } -} +void to_adjoint(std::vector& output, std::vector& U, bool sign, size_t b, size_t n, size_t n_squared); } // namespace internal /** @@ -124,26 +47,35 @@ void to_adjoint(T* output, std::vector>& U, bool sign, size_t b, template void inverse(const T* input, T* output, const Shape& shape, const bool adjoint) { const size_t n = shape.back(); - const auto total_elements = shape_size(shape); - const auto batch_size = total_elements / n / n; + const auto n_squared = n * n; + size_t batch_size = 1; + + for (size_t i = 0; i < shape.size() - 2; ++i) { + batch_size = batch_size * shape[i]; + } + + auto input_conv = std::vector(batch_size * n_squared); + auto output_conv = std::vector(batch_size * n_squared); + + convert(input, input_conv.data(), batch_size * n_squared); + + std::vector L(n_squared); + std::vector U(n_squared); + std::vector P(n); for (size_t b = 0; b < batch_size; ++b) { - std::vector> L(n, std::vector(n, T{0})); - std::vector> U(n, std::vector(n, T{0})); - std::vector P(n); bool sign = true; - internal::lu_decomposition(input, L, U, P, sign, b, n); + internal::lu_decomposition(input_conv, L, U, P, sign, b, n, n_squared); - for (size_t column = 0; column < n; ++column) { - internal::lu_solve(output, L, U, P, b, n, column); - } + internal::lu_solve(output_conv, L, U, P, b, n, n_squared); if (adjoint) { - internal::to_adjoint(output, U, sign, b, n); + internal::to_adjoint(output_conv, U, sign, b, n, n_squared); } } + + convert(output_conv.data(), output, batch_size * n_squared); } -} // namespace inverse } // namespace reference } // namespace ov diff --git a/src/core/reference/src/op/inverse.cpp b/src/core/reference/src/op/inverse.cpp new file mode 100644 index 00000000000..e02c4aeb996 --- /dev/null +++ b/src/core/reference/src/op/inverse.cpp @@ -0,0 +1,124 @@ +// Copyright (C) 2018-2024 Intel Corporation +// SPDX-License-Identifier: Apache-2.0 +// + +#include "openvino/reference/inverse.hpp" + +namespace ov { +namespace reference { +namespace internal { + +void lu_decomposition(std::vector& input, + std::vector& L, + std::vector& U, + std::vector& P, + bool& sign, + size_t b, + size_t n, + size_t n_squared) { + // Make L identity, U a copy of input and P a range(0, n) + const auto batch_idx = b * n_squared; + + std::fill(L.begin(), L.end(), 0.0f); + memcpy(&U[0], &input[batch_idx], sizeof(float) * n_squared); + + for (size_t i = 0; i < n; ++i) { + L[i * n + i] = 1.0f; + P[i] = i; + } + + for (size_t k = 0; k < n; ++k) { + // Partial Pivoting + auto pivot_row = k; + auto pivot_idx = pivot_row * n; + const auto k_idx = k * n; + + for (auto i = (k + 1) * n, j = k + 1; i < n_squared; i += n, ++j) { + if (std::abs(U[i + k]) > std::abs(U[pivot_idx + k])) { + pivot_row = j; + pivot_idx = pivot_row * n; + } + } + + if (pivot_row != k) { + // Swap rows in L, U (A) and P + sign = !sign; + std::swap(P[k], P[pivot_row]); + std::swap_ranges(&U[k_idx], &U[k_idx + n], &U[pivot_idx]); + std::swap_ranges(&L[k_idx], &L[k_idx + n], &L[pivot_idx]); + } + + const auto remaining_columns = n - k; + const auto remaining_rows = remaining_columns - 1; + + for (size_t i = 0; i < remaining_rows; ++i) { + const auto i_idx = (i + k + 1) * n; + L[i_idx + k] = U[i_idx + k] / U[k_idx + k]; + } + + for (size_t i = 0; i < remaining_rows * remaining_columns; ++i) { + const auto i_idx = (i / remaining_columns + k + 1) * n; + const auto j_idx = i % remaining_columns + k; + U[i_idx + j_idx] = U[i_idx + j_idx] - L[i_idx + k] * U[k_idx + j_idx]; + } + } +} + +void lu_solve(std::vector& output, + std::vector& L, + std::vector& U, + std::vector& P, + size_t b, + size_t n, + size_t n_squared) { + std::vector X(n); + std::vector Y(n); + + for (size_t column = 0; column < n; ++column) { + std::fill(X.begin(), X.end(), 0.0f); + std::fill(Y.begin(), Y.end(), 0.0f); + + // Forward substitution: Ly = Pb + for (size_t i = 0; i < n; ++i) { + if (P[i] == column) { + Y[i] = 1.0f; + } + const auto i_idx = i * n; + for (size_t j = 0; j < i; ++j) { + Y[i] -= L[i_idx + j] * Y[j]; + } + } + + // Backward substitution: Ux = y + for (size_t i = 0; i < n; ++i) { + const auto i_adj = n - i - 1; + const auto i_idx = i_adj * n; + X[i_adj] = Y[i_adj]; + for (size_t j = i_adj + 1; j < n; ++j) { + X[i_adj] = X[i_adj] - U[i_idx + j] * X[j]; + } + X[i_adj] = X[i_adj] / U[i_idx + i_adj]; + } + + const auto batch_column_idx = b * n_squared + column; + for (size_t row = 0; row < n; ++row) { + output[batch_column_idx + row * n] = X[row]; + } + } +} + +void to_adjoint(std::vector& output, std::vector& U, bool sign, size_t b, size_t n, size_t n_squared) { + auto determinant = sign ? 1.0f : -1.0f; + + for (size_t i = 0; i < n; ++i) { + determinant *= U[i * n + i]; + } + + const auto batch_idx = b * n_squared; + for (size_t i = 0; i < n_squared; ++i) { + output[batch_idx + i] *= determinant; + } +} +} // namespace internal +} // namespace reference +} // namespace ov diff --git a/src/core/tests/type_prop/inverse.cpp b/src/core/tests/type_prop/inverse.cpp index 70d18ae0b4f..1a60993a1dd 100644 --- a/src/core/tests/type_prop/inverse.cpp +++ b/src/core/tests/type_prop/inverse.cpp @@ -33,6 +33,17 @@ TEST_F(TypePropInverseV14Test, input_f64_2x2_constant) { EXPECT_EQ(op->get_output_partial_shape(0), (ov::PartialShape{2, 2})); } +TEST_F(TypePropInverseV14Test, symbol_propagation) { + auto input_shape = ov::PartialShape{2, 2}; + auto A = std::make_shared(), B = std::make_shared(); + set_shape_symbols(input_shape, {A, B}); + const auto data = std::make_shared(ov::element::f32, input_shape); + const auto op = make_op(data, false); + EXPECT_EQ(op->get_element_type(), ov::element::f32); + EXPECT_EQ(op->get_output_partial_shape(0), (input_shape)); + EXPECT_EQ(get_shape_symbols(op->get_output_partial_shape(0)), get_shape_symbols(input_shape)); +} + TEST_F(TypePropInverseV14Test, input_f32_static) { const auto data = std::make_shared(ov::element::f32, ov::Shape{10, 4, 4}); const auto op = make_op(data, false); diff --git a/src/plugins/intel_cpu/src/nodes/inverse.cpp b/src/plugins/intel_cpu/src/nodes/inverse.cpp index 04c283fc2ff..6f95482c418 100644 --- a/src/plugins/intel_cpu/src/nodes/inverse.cpp +++ b/src/plugins/intel_cpu/src/nodes/inverse.cpp @@ -11,9 +11,6 @@ #include "openvino/op/inverse.hpp" #include "utils/bfloat16.hpp" -// Parallel LU decomposition algorithm with partial pivoting -// Based on the lectures by Prof. Dr. Thomas Huckle, Parallel Numerics - namespace ov { namespace intel_cpu { namespace node { @@ -56,7 +53,7 @@ void Inverse::getSupportedDescriptors() { void Inverse::initSupportedPrimitiveDescriptors() { m_input_precision = getOriginalInputPrecisionAtPort(INPUT_PORT); - if (!one_of(m_input_precision, ov::element::f32, ov::element::f16, ov::element::bf16)) { + if (m_input_precision != ov::element::f32) { m_input_precision = ov::element::f32; } @@ -88,68 +85,56 @@ bool Inverse::created() const { } void Inverse::execute(dnnl::stream strm) { - OV_SWITCH(intel_cpu, - InverseExecute, - this, - m_input_precision, - OV_CASE(ov::element::bf16, bfloat16_t), - OV_CASE(ov::element::f16, ov::float16), - OV_CASE(ov::element::f32, float)) + inverse(); } void Inverse::executeDynamicImpl(dnnl::stream strm) { execute(strm); } -template void Inverse::inverse() { - const auto* data = getSrcDataAtPortAs(INPUT_PORT); - auto* output = getDstDataAtPortAs(OUTPUT_PORT); + const auto* data = getSrcDataAtPortAs(INPUT_PORT); + auto* output = getDstDataAtPortAs(OUTPUT_PORT); - std::vector L(m_side_squared, T{0}); - std::vector U(m_side_squared, T{0}); - std::vector P(m_side, T{0}); + std::vector L(m_side_squared); + std::vector U(m_side_squared); + std::vector P(m_side); for (size_t b = 0; b < m_batches_count; ++b) { bool sign = true; + lu_decomposition(data, L, U, P, sign, b); - for (size_t column = 0; column < m_side; ++column) { - lu_solve(output, L, U, P, b, column); - } + lu_solve(output, L, U, P, b); if (m_adjoint) { - // Multiply by det(A) = det(U) to_adjoint(output, U, sign, b); } } } -template -void Inverse::lu_decomposition(const T* data, - std::vector& L, - std::vector& U, - std::vector& P, +void Inverse::lu_decomposition(const float* data, + std::vector& L, + std::vector& U, + std::vector& P, bool& sign, size_t b) { // Make L identity, U a copy of data and P a range(0, side) - size_t batch_idx = b * m_side_squared; + const auto batch_idx = b * m_side_squared; - std::fill(L.begin(), L.end(), T{0}); - cpu_parallel_memcpy(&U[0], &data[batch_idx], sizeof(T) * m_side_squared); + std::fill(L.begin(), L.end(), 0.0f); + cpu_parallel_memcpy(&U[0], &data[batch_idx], sizeof(float) * m_side_squared); parallel_for(m_side, [&](size_t i) { - L[i * m_side + i] = T{1}; - P[i] = static_cast(i); + L[i * m_side + i] = 1.0f; + P[i] = i; }); for (size_t k = 0; k < m_side; ++k) { - size_t pivot_row = k; - size_t k_idx = k * m_side; - size_t pivot_idx = pivot_row * m_side; - - size_t remaining_columns = m_side - k; - size_t remaining_rows = remaining_columns - 1; + // Partial Pivoting + auto pivot_row = k; + auto pivot_idx = pivot_row * m_side; + const auto k_idx = k * m_side; // Find maximum value pivot - non-parallel for (size_t i = (k + 1) * m_side, j = k + 1; i < m_side_squared; i += m_side, ++j) { @@ -169,63 +154,65 @@ void Inverse::lu_decomposition(const T* data, }); } + const auto remaining_columns = m_side - k; + const auto remaining_rows = remaining_columns - 1; + parallel_for(remaining_rows, [&](size_t i) { - size_t i_idx = (i + k + 1) * m_side; + const auto i_idx = (i + k + 1) * m_side; L[i_idx + k] = U[i_idx + k] / U[k_idx + k]; }); parallel_for(remaining_rows * remaining_columns, [&](size_t i) { - size_t i_idx = (i / remaining_columns + k + 1) * m_side; - size_t j_idx = i % remaining_columns + k; - + const auto i_idx = (i / remaining_columns + k + 1) * m_side; + const auto j_idx = i % remaining_columns + k; U[i_idx + j_idx] = U[i_idx + j_idx] - L[i_idx + k] * U[k_idx + j_idx]; }); } } -template -void Inverse::lu_solve(T* output, std::vector& L, std::vector& U, std::vector& P, size_t b, size_t column) { - std::vector B(m_side, T{0}); - std::vector X(m_side, T{0}); - std::vector Y(m_side, T{0}); - B[column] = T{1}; +void Inverse::lu_solve(float* output, std::vector& L, std::vector& U, std::vector& P, size_t b) { + parallel_for(m_side, [&](size_t column) { + std::vector X(m_side, 0.0f); + std::vector Y(m_side, 0.0f); - // Forward substitution: Ly = Pb - not possible to be parallel - for (size_t i = 0; i < m_side; ++i) { - Y[i] = B[P[i]]; - size_t i_idx = i * m_side; - for (size_t j = 0; j < i; ++j) { - Y[i] = Y[i] - L[i_idx + j] * Y[j]; + // Forward substitution: Ly = Pb + for (size_t i = 0; i < m_side; ++i) { + if (P[i] == column) { + Y[i] = 1.0f; + } + const auto i_idx = i * m_side; + for (size_t j = 0; j < i; ++j) { + Y[i] = Y[i] - L[i_idx + j] * Y[j]; + } } - } - // Backward substitution: Ux = y - not possible to be parallel - for (size_t i = 0; i < m_side; ++i) { - size_t i_adj = m_side - i - 1; - size_t i_idx = i_adj * m_side; - X[i_adj] = Y[i_adj]; - for (size_t j = i_adj + 1; j < m_side; ++j) { - X[i_adj] = X[i_adj] - U[i_idx + j] * X[j]; + // Backward substitution: Ux = y + for (size_t i = 0; i < m_side; ++i) { + size_t i_adj = m_side - i - 1; + size_t i_idx = i_adj * m_side; + X[i_adj] = Y[i_adj]; + for (size_t j = i_adj + 1; j < m_side; ++j) { + X[i_adj] = X[i_adj] - U[i_idx + j] * X[j]; + } + X[i_adj] = X[i_adj] / U[i_idx + i_adj]; } - X[i_adj] = X[i_adj] / U[i_idx + i_adj]; - } - // Substitute back to get result - size_t batch_idx = b * m_side_squared; - parallel_for(m_side, [&](size_t row) { - output[batch_idx + row * m_side + column] = X[row]; + // Substitute back to get result + const auto batch_column_idx = b * m_side_squared + column; + for (size_t row = 0; row < m_side; ++row) { + output[batch_column_idx + row * m_side] = X[row]; + } }); } -template -void Inverse::to_adjoint(T* output, std::vector& U, bool sign, size_t b) { - T determinant = sign ? T{1} : T{-1}; +void Inverse::to_adjoint(float* output, std::vector& U, bool sign, size_t b) { + auto determinant = sign ? 1.0f : -1.0f; for (size_t i = 0; i < m_side; ++i) { determinant = determinant * U[i * m_side + i]; } - size_t batch_idx = b * m_side_squared; + const auto batch_idx = b * m_side_squared; parallel_for(m_side_squared, [&](size_t i) { output[batch_idx + i] = output[batch_idx + i] * determinant; }); diff --git a/src/plugins/intel_cpu/src/nodes/inverse.hpp b/src/plugins/intel_cpu/src/nodes/inverse.hpp index 223f3fd9313..c87122955ce 100644 --- a/src/plugins/intel_cpu/src/nodes/inverse.hpp +++ b/src/plugins/intel_cpu/src/nodes/inverse.hpp @@ -48,24 +48,18 @@ private: size_t m_batches_count = 0; // Helper functions - template void inverse(); - template - void lu_decomposition(const T* data, std::vector& L, std::vector& U, std::vector& P, bool& sign, size_t b); + void lu_decomposition(const float* data, + std::vector& L, + std::vector& U, + std::vector& P, + bool& sign, + size_t b); - template - void to_adjoint(T* output, std::vector& U, bool sign, size_t b); + void to_adjoint(float* output, std::vector& U, bool sign, size_t b); - template - void lu_solve(T* output, std::vector& L, std::vector& U, std::vector& P, size_t b, size_t column); - - template - struct InverseExecute { - void operator()(Inverse* node) { - node->inverse(); - } - }; + void lu_solve(float* output, std::vector& L, std::vector& U, std::vector& P, size_t b); }; } // namespace node diff --git a/src/plugins/intel_cpu/tests/functional/custom/single_layer_tests/instances/common/inverse.cpp b/src/plugins/intel_cpu/tests/functional/custom/single_layer_tests/instances/common/inverse.cpp index 59fea8f9afa..734971ab97f 100644 --- a/src/plugins/intel_cpu/tests/functional/custom/single_layer_tests/instances/common/inverse.cpp +++ b/src/plugins/intel_cpu/tests/functional/custom/single_layer_tests/instances/common/inverse.cpp @@ -8,22 +8,23 @@ namespace { using ov::test::InverseLayerTest; -const auto shapes = testing::Values(ov::Shape{10, 10}, ov::Shape{5, 7, 7}, ov::Shape{5, 4, 3, 3}, ov::Shape{100, 2, 2}); +const std::vector> input_shapes = { + {{{5, 4, 4}, {{5, 4, 4}}}}, + {{{20, 3, 3}, {{20, 3, 3}}}}, + {{{ov::Dimension{1, 70}, -1}, {{3, 3}, {5, 5}, {4, 4}}}}, + {{{-1, ov::Dimension{1, 70}, -1}, {{3, 3, 3}, {4, 4, 4}, {5, 5, 5}}}}}; + +const auto shapes = testing::ValuesIn(input_shapes); const auto dtypes = testing::Values(ov::element::f32, ov::element::f16, ov::element::bf16); const auto adjoint = testing::Values(false, true); -const auto test_static = testing::Values(true); -const auto test_dynamic = testing::Values(false); - -const auto seed = testing::Values(0u, 1u, 3u); +const auto seed = testing::Values(1, 2, 3); const auto device_cpu = testing::Values(ov::test::utils::DEVICE_CPU); -const auto params_static = ::testing::Combine(shapes, dtypes, adjoint, test_static, seed, device_cpu); -const auto params_dynamic = ::testing::Combine(shapes, dtypes, adjoint, test_dynamic, seed, device_cpu); +const auto params = ::testing::Combine(shapes, dtypes, adjoint, seed, device_cpu); -INSTANTIATE_TEST_SUITE_P(smoke_InverseStatic, InverseLayerTest, params_static, InverseLayerTest::getTestCaseName); -INSTANTIATE_TEST_SUITE_P(smoke_InverseDynamic, InverseLayerTest, params_dynamic, InverseLayerTest::getTestCaseName); +INSTANTIATE_TEST_SUITE_P(smoke_InverseStatic, InverseLayerTest, params, InverseLayerTest::getTestCaseName); } // namespace diff --git a/src/plugins/intel_cpu/tests/functional/shared_tests_instances/skip_tests_config.cpp b/src/plugins/intel_cpu/tests/functional/shared_tests_instances/skip_tests_config.cpp index 800500111d5..ca86d6a69c2 100644 --- a/src/plugins/intel_cpu/tests/functional/shared_tests_instances/skip_tests_config.cpp +++ b/src/plugins/intel_cpu/tests/functional/shared_tests_instances/skip_tests_config.cpp @@ -197,9 +197,7 @@ std::vector disabledTestPatterns() { R"(.*smoke_LPT/InterpolateTransformation.*)", // Issue: 129931 R"(smoke_LPT/ConvolutionTransformation.CompareWithRefImpl/f32_\[.*,3,16,16\]_CPU_f32_rank=4D_fq_on_data=\{level=256_shape=\[1\]_input_low=\{ 0 \}_input_high=\{ 255 \}_output_low=\{ .*18.7 \}_output_high\{ 18.8 \}_precision=\}_fq_on_weights=\{_255_\[6,1,1,1\]_\{ .*1.52806e.*39, .*0.2, .*0.3, .*0.3, .*0.2, .*0.1 \}_\{ 1.52806e.*39, 0.2, 0.3, 0.3, 0.2, 0.1 \}\})", - // Issue: 132494 - R"(.*smoke_Inverse.*bf16.*)", - // Issue: CVS-133173 + // Issue: 133173 R"(.*smoke_ScaledAttn_CPU/ScaledAttnLayerCPUTest.CompareWithRefs/netPRC=bf16.*has_scale=0.*)", R"(.*smoke_LPT_4D/ConvolutionBackpropDataTransformation.CompareWithRefImpl/f32_\[1,8,16,16\]_CPU_f32_\[16,16\]_level=256_shape=\[.*\]_input_low=\{ 0 \}_input_high=\{ 25.5 \}_output_low=\{ 0 \}_output_high\{ 25.5 \}_precision=__255_\[.*\]_\{ -12.7 \}_\{ 12.7 \}_\{\}.*)", R"(.*smoke_LPT_4D/ConvolutionBackpropDataTransformation.CompareWithRefImpl/f32_\[1,8,16,16\]_CPU_f32_\[16,16\]_level=256_shape=\[1,1,1,1\]_input_low=\{ 0 \}_input_high=\{ 255 \}_output_low=\{ -12.7 \}_output_high\{ 12.8 \}_precision=.*)", diff --git a/src/plugins/template/backend/ops/inverse.cpp b/src/plugins/template/backend/ops/inverse.cpp index 38b24c06b32..0490759bc7f 100644 --- a/src/plugins/template/backend/ops/inverse.cpp +++ b/src/plugins/template/backend/ops/inverse.cpp @@ -17,7 +17,7 @@ inline bool evaluate(const std::shared_ptr& op, const auto out_shape = ov::op::v14::shape_infer(op.get(), input_shapes).front().to_shape(); outputs[0].set_shape(out_shape); - ov::reference::inverse::inverse(inputs[0].data(), outputs[0].data(), out_shape, op->get_adjoint()); + ov::reference::inverse(inputs[0].data(), outputs[0].data(), out_shape, op->get_adjoint()); return true; } diff --git a/src/tests/functional/shared_test_classes/include/shared_test_classes/single_op/inverse.hpp b/src/tests/functional/shared_test_classes/include/shared_test_classes/single_op/inverse.hpp index 090ea253075..895875a73c2 100644 --- a/src/tests/functional/shared_test_classes/include/shared_test_classes/single_op/inverse.hpp +++ b/src/tests/functional/shared_test_classes/include/shared_test_classes/single_op/inverse.hpp @@ -13,12 +13,11 @@ namespace ov { namespace test { -using InverseTestParams = typename std::tuple, // input shape + ov::element::Type, // element type + bool, // adjoint + int32_t, // seed + std::string // device_name >; class InverseLayerTest : public testing::WithParamInterface, virtual public SubgraphBaseTest { @@ -31,7 +30,7 @@ protected: void generate_inputs(const std::vector& target_shapes) override; private: - unsigned int m_seed; + int32_t m_seed; }; } // namespace test } // namespace ov diff --git a/src/tests/functional/shared_test_classes/src/single_op/inverse.cpp b/src/tests/functional/shared_test_classes/src/single_op/inverse.cpp index 6c2caacd7f6..050c7205073 100644 --- a/src/tests/functional/shared_test_classes/src/single_op/inverse.cpp +++ b/src/tests/functional/shared_test_classes/src/single_op/inverse.cpp @@ -7,26 +7,26 @@ #include #include +#include "common_test_utils/ov_tensor_utils.hpp" + using namespace ov::test; namespace ov { namespace test { std::string InverseLayerTest::getTestCaseName(const testing::TestParamInfo& obj) { - ov::Shape input_shape; + std::vector input_shape; ov::element::Type element_type; bool adjoint; - bool test_static; - unsigned int seed; + int32_t seed; std::string device_name; - std::tie(input_shape, element_type, adjoint, test_static, seed, device_name) = obj.param; + std::tie(input_shape, element_type, adjoint, seed, device_name) = obj.param; const char separator = '_'; std::ostringstream result; - result << "IS=" << input_shape.to_string() << separator; + result << "IS=" << input_shape[0].first.to_string() << separator; result << "dtype=" << element_type.to_string() << separator; result << "adjoint=" << ov::test::utils::bool2str(adjoint) << separator; - result << "static=" << ov::test::utils::bool2str(test_static) << separator; result << "seed=" << seed << separator; result << "device=" << device_name; @@ -34,68 +34,36 @@ std::string InverseLayerTest::getTestCaseName(const testing::TestParamInfo input_shape; ov::element::Type element_type; bool adjoint; - bool test_static; - std::tie(input_shape, element_type, adjoint, test_static, m_seed, targetDevice) = GetParam(); + std::tie(input_shape, element_type, adjoint, m_seed, targetDevice) = GetParam(); - std::vector in_shapes; - if (!test_static) { - in_shapes.push_back({{}, {input_shape}}); - } else { - in_shapes.push_back({input_shape, {input_shape}}); - } - init_input_shapes(in_shapes); + init_input_shapes(input_shape); - const auto data = std::make_shared(element_type, input_shape); + const auto data = std::make_shared(element_type, input_shape[0].first); data->set_friendly_name("data"); - auto inverse = std::make_shared(data, adjoint); + const auto inverse = std::make_shared(data, adjoint); ov::ResultVector results{std::make_shared(inverse)}; function = std::make_shared(results, ParameterVector{data}, "InverseTestCPU"); } -template -void fill_data(T* dst, const size_t count, const unsigned int seed) { - std::mt19937 gen(seed); - std::uniform_real_distribution dis(0.0f, 10.0f); - - for (size_t i = 0; i < count; i++) { - dst[i] = static_cast(dis(gen)); - } -} - void InverseLayerTest::generate_inputs(const std::vector& targetInputStaticShapes) { inputs.clear(); const auto& func_inputs = function->inputs(); const auto& func_input = func_inputs[0]; const auto& name = func_input.get_node()->get_friendly_name(); const auto& in_prc = func_input.get_element_type(); - auto tensor = ov::Tensor(in_prc, targetInputStaticShapes[0]); - size_t count = std::accumulate(targetInputStaticShapes[0].begin(), - targetInputStaticShapes[0].end(), - 1, - std::multiplies()); - - switch (in_prc) { - case ov::element::f32: - fill_data(tensor.data::value_type>(), count, m_seed); - break; - case ov::element::f16: - fill_data(tensor.data::value_type>(), count, m_seed); - break; - case ov::element::bf16: - fill_data(tensor.data::value_type>(), count, m_seed); - break; - default: - OPENVINO_THROW("Inverse does not support precision ", in_prc, " for the 'data' input."); - } + ov::test::utils::InputGenerateData in_data; + in_data.start_from = 5; + in_data.range = 5; + in_data.resolution = 1; + in_data.seed = m_seed; + auto tensor = ov::test::utils::create_and_fill_tensor(in_prc, targetInputStaticShapes[0], in_data); inputs.insert({func_input.get_node_shared_ptr(), tensor}); } } // namespace test diff --git a/tests/layer_tests/pytorch_tests/test_inverse.py b/tests/layer_tests/pytorch_tests/test_inverse.py index 9afa1b148a0..e8365520655 100644 --- a/tests/layer_tests/pytorch_tests/test_inverse.py +++ b/tests/layer_tests/pytorch_tests/test_inverse.py @@ -77,3 +77,32 @@ class TestInverse(PytorchLayerTest): self._test(aten_inverse_out(), None, "aten::linalg_inv", ie_device, precision, ir_version, trace_model=True, freeze_model=False, kwargs_to_prepare_input={"out": out}) + @pytest.mark.parametrize("shape", [ + (10, 2, 2), + (3, 5, 5), + (10, 10), + (7, 6, 5, 4, 4) + ]) + @pytest.mark.parametrize("dtype", [ + np.float64, + np.float32 + ]) + @pytest.mark.parametrize("seed", [ + 1, 2, 3 + ]) + @pytest.mark.parametrize("out", [ + False, + True + ]) + @pytest.mark.nightly + @pytest.mark.precommit + def test_inverse(self, shape, dtype, seed, out, ie_device, precision, ir_version): + rng = np.random.default_rng(seed) + self.input_tensor = rng.uniform(-10.0, 10.0, shape).astype(dtype) + if not out: + self._test(aten_inverse(), None, "aten::linalg_inv", + ie_device, precision, ir_version, trace_model=True, freeze_model=False) + else: + self._test(aten_inverse_out(), None, "aten::linalg_inv", + ie_device, precision, ir_version, trace_model=True, freeze_model=False, kwargs_to_prepare_input={"out": out}) +