Default.Dot.cs

using System;
using System.Diagnostics;
using System.Diagnostics.CodeAnalysis;
using System.Linq;
using System.Runtime.CompilerServices;
using System.Threading.Tasks;
using NumSharp.Utilities;

namespace NumSharp.Backends
{
    public partial class DefaultEngine
    {

        /// <remarks>https://numpy.org/doc/stable/reference/generated/numpy.dot.html</remarks>
        public override NDArray Dot(NDArray left, NDArray right)
        {
            //Dot product of two arrays.Specifically,
            //If both a and b are 1 - D arrays, it is inner product of vectors(without complex conjugation).
            //If both a and b are 2 - D arrays, it is matrix multiplication, but using matmul or a @ b is preferred.
            //V If either a or b is 0 - D(scalar), it is equivalent to multiply and using numpy.multiply(a, b) or a *b is preferred.
            //If a is an N - D array and b is a 1 - D array, it is a sum product over the last axis of a and b.
            //If a is an N - D array and b is an M - D array(where M >= 2), it is a sum product over the last axis of a and the second-to - last axis of b:
            //  dot(a, b)[i, j, k, m] = sum(a[i, j,:] * b[k,:, m])
            // An external BLAS, when one is installed (see TensorEngine.Blas). np.dot is NOT
            // matmul's dispatcher — the two disagree on e.g. a strided matrix times a vector — so
            // a backend gets its own entry point here, ahead of every shape branch below.
            // Read the property ONCE: it is settable at any time from any thread, and testing it
            // then calling it would be a null-dereference every time a concurrent OpenBlasEngine.Disable()
            // landed between the two reads (measured: ~2% of products under a hostile flipper).
            var blas = Blas;
            NDArray blasResult;
            if (blas != null && blas.TryDot(left, right, out blasResult))
                return blasResult;

            var leftshape = left.Shape;
            var rightshape = right.Shape;
            var isLeftScalar = leftshape.IsScalar;
            var isRightScalar = rightshape.IsScalar;

            if (isLeftScalar && isRightScalar)
            {
                return Multiply(left, right);
            }

            //If either a or b is 0-D (scalar), it is equivalent to multiply and using numpy.multiply(a, b) or a * b is preferred.
            if (isLeftScalar || isRightScalar)
            {
                return Multiply(left, right);
            }

            //If both a and b are 2-D arrays, it is matrix multiplication, but using matmul or a @ b is preferred.
            if (leftshape.NDim == 2 && rightshape.NDim == 2)
            {
                return MultiplyMatrix(left, right);
            }

            //If both a and b are 1-D arrays, it is inner product of vectors (without complex conjugation).
            // Fused single-pass kernel (no temp product array) — see Default.Dot.Fused.cs.
            if (leftshape.NDim == 1 && rightshape.NDim == 1)
            {
                return DotInner1D(left, right);
            }

            //If a is an N-D array and b is a 1-D array, it is a sum product over the last axis of a and b.
            // This is exactly matmul's contraction (a's last axis with the promoted column vector), so
            // route through Matmul — it preserves the input dtype (NumPy: dot(int32, int32) -> int32)
            // and contracts the correct (last) axis, where np.sum used the wider accumulator (int64)
            // and a hard-coded axis=1 (wrong for N-D).
            if (leftshape.NDim >= 2 && rightshape.NDim == 1)
            {
                return Matmul(left, right);
            }
            // If a is 1-D and b is 2-D, treat a as row vector and do matrix multiply, then squeeze
            // (n,) dot (n, m) -> (m,)
            if (leftshape.NDim == 1 && rightshape.NDim == 2)
            {
                if (leftshape[0] != rightshape[0])
                    throw new ArgumentException($"shapes ({leftshape[0]},) and ({rightshape[0]},{rightshape[1]}) not aligned: {leftshape[0]} (dim 0) != {rightshape[0]} (dim 0)");

                // Reshape 1D to row vector (1, n), multiply, squeeze back to 1D
                var leftReshaped = left.reshape(1, leftshape[0]);
                var result = MultiplyMatrix(leftReshaped, right);
                var squeezed = result.reshape(rightshape[1]);
                leftReshaped.Dispose(); // row-vector alias of the input — the input keeps its own ref
                result.Dispose();       // the returned view holds its own counted ref on the product buffer
                return squeezed;
            }

            //left cant be 0 or 1 by this point
            //If a is an N-D array and b is an M-D array (where M>=2), it is a sum product over the last axis of a and the second-to-last axis of b:
            //dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])
            if (rightshape.NDim >= 2)
            {
                return DotNDMD(left, right);
            }

            throw new NotSupportedException();
        }
    }
}