3]_ and [4]_ have more reference
implementations in Python.

.. versionadded:: 1.6.0

References
----------
.. [1] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian
   Q Weinberger, Andrew Gordon Wilson, "GPyTorch: Blackbox Matrix-Matrix
   Gaussian Process Inference with GPU Acceleration" with contributions
   from Max Balandat and Ruihan Wu. Available online:
   https://github.com/cornellius-gp/gpytorch

.. [2] J. Demmel, P. Koev, and X. Li, "A Brief Survey of Direct Linear
   Solvers". In Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der
   Vorst, editors. Templates for the Solution of Algebraic Eigenvalue
   Problems: A Practical Guide. SIAM, Philadelphia, 2000. Available at:
   http://www.netlib.org/utk/people/JackDongarra/etemplates/node384.html

.. [3] R. Scheibler, E. Bezzam, I. Dokmanic, Pyroomacoustics: A Python
   package for audio room simulations and array processing algorithms,
   Proc. IEEE ICASSP, Calgary, CA, 2018.
   https://github.com/LCAV/pyroomacoustics/blob/pypi-release/
   pyroomacoustics/adaptive/util.py

.. [4] Marano S, Edwards B, Ferrari G and Fah D (2017), "Fitting
   Earthquake Spectra: Colored Noise and Incomplete Data", Bulletin of
   the Seismological Society of America., January, 2017. Vol. 107(1),
   pp. 276-291.

Examples
--------
Multiply the Toeplitz matrix T with matrix x::

        [ 1 -1 -2 -3]       [1 10]
    T = [ 3  1 -1 -2]   x = [2 11]
        [ 6  3  1 -1]       [2 11]
        [10  6  3  1]       [5 19]

To specify the Toeplitz matrix, only the first column and the first
row are needed.

>>> import numpy as np
>>> c = np.array([1, 3, 6, 10])    # First column of T
>>> r = np.array([1, -1, -2, -3])  # First row of T
>>> x = np.array([[1, 10], [2, 11], [2, 11], [5, 19]])

>>> from scipy.linalg import toeplitz, matmul_toeplitz
>>> matmul_toeplitz((c, r), x)
array([[-20., -80.],
       [ -7.,  -8.],
       [  9.,  85.],
       [ 33., 218.]])

Check the result by creating the full Toeplitz matrix and
multiplying it by ``x``.

>>> toeplitz(c, r) @ x
array([[-20, -80],
       [ -7,  -8],
       [  9,  85],
       [ 33, 218]])

The full matrix is never formed explicitly, so this routine
is suitable for very large Toeplitz matrices.

>>> n = 1000000
>>> matmul_toeplitz([1] + [0]*(n-1), np.ones(n))
array([1., 1., 1., ..., 1., 1., 1.], shape=(1000000,))

rG