def solve(a: ndarray, b: ndarray, out: ndarray | None = None) -> ndarray:
    """
    Solve a linear matrix equation, or system of linear scalar equations.

    Computes the "exact" solution, `x`, of the well-determined, i.e., full
    rank, linear matrix equation `ax = b`.

    Parameters
    ----------
    a : (M, M) array_like
        Coefficient matrix.
    b : {(M,), (M, K)}, array_like
        Ordinate or "dependent variable" values.
    out : {(M,), (M, K)}, array_like, optional
        An optional output array for the solution

    Returns
    -------
    x : {(M,), (M, K)} ndarray
        Solution to the system a x = b.  Returned shape is identical to `b`.

    Raises
    ------
    LinAlgError
        If `a` is singular or not square.

    Notes
    ------
    Multi-GPU usage is only available when compiled with cusolverMP.

    See Also
    --------
    numpy.linalg.solve

    Availability
    --------
    Multiple GPUs, Single CPU
    """
    if a.ndim < 2:
        raise LinAlgError(
            f"{a.ndim}-dimensional array given. "
            "Array must be at least two-dimensional"
        )
    if b.ndim < 1:
        raise LinAlgError(
            f"{b.ndim}-dimensional array given. "
            "Array must be at least one-dimensional"
        )
    if np.dtype("e") in (a.dtype, b.dtype):
        raise TypeError("array type float16 is unsupported in linalg")
    if a.ndim > 2 or b.ndim > 2:
        raise NotImplementedError(
            "cuPyNumeric does not yet support stacked 2d arrays"
        )
    if a.shape[-2] != a.shape[-1]:
        raise LinAlgError("Last 2 dimensions of the array must be square")
    if a.shape[-1] != b.shape[0]:
        if b.ndim == 1:
            raise ValueError(
                "Input operand 1 has a mismatch in its dimension 0, "
                f"with signature (m,m),(m)->(m) (size {b.shape[0]} "
                f"is different from {a.shape[-1]})"
            )
        else:
            raise ValueError(
                "Input operand 1 has a mismatch in its dimension 0, "
                f"with signature (m,m),(m,n)->(m,n) (size {b.shape[0]} "
                f"is different from {a.shape[-1]})"
            )
    if a.size == 0 or b.size == 0:
        return empty_like(b)

    return _thunk_solve(a, b, out)

    Solve a linear matrix equation, or system of linear scalar equations.

    Computes the "exact" solution, `x`, of the well-determined, i.e., full
    rank, linear matrix equation `ax = b`.

    Parameters
    ----------
    a : (M, M) array_like
        Coefficient matrix.
    b : {(M,), (M, K)}, array_like
        Ordinate or "dependent variable" values.
    out : {(M,), (M, K)}, array_like, optional
        An optional output array for the solution

    Returns
    -------
    x : {(M,), (M, K)} ndarray
        Solution to the system a x = b.  Returned shape is identical to `b`.

    Raises
    ------
    LinAlgError
        If `a` is singular or not square.

    Notes
    ------
    Multi-GPU usage is only available when compiled with cusolverMP.

    See Also
    --------
    numpy.linalg.solve

    Availability
    --------
    Multiple GPUs, Single CPU
    
cupynumeric/linalg/linalg.py