scipy: [Bug] sparse.linalg.LinearOperator.matmat cannot fallback properly to matvec

Hi, LinearOperator cannot properly use matvec to implement matmat. Maybe somehow connected to issue https://github.com/scipy/scipy/issues/2434 because these calls seem to be a bit complicated: matmat -> _matmat -> matvec -> _matvec -> matmat ...

Reproducing code example:

import numpy as np
import scipy as sp
import scipy.sparse.linalg as sla


dim = 8
shp = (dim, dim)

I = np.eye(*shp)

x0 = -sp.pi
xf = sp.pi
x = np.linspace(x0, xf, dim, endpoint=False)
# this is a vector representation of an operator with diagonal matrix
ref = (1-np.exp(-x))**2


def mv(y):
    # multiply this diag matrix by vector
    vv = ref
    return vv*y


def mm(z):
    # some crutches: explicit matmat function
   return np.array([mv(y) for y in z.T])

LO_with_mm = sla.LinearOperator(shp, matvec=mv, matmat=mm)
mm_test = LO_with_mm.matmat(I)
# multiplying by identity should result in the initial diagonal
print('works with explicit matmat:', mm_test - np.diagflat(ref))

# this won't work
LO_withot_mm = sla.LinearOperator(shp, matvec=mv)
mm_test = LO_withot_mm.matmat(I)

Error message:

traceback (most recent call last):

  File "<ipython-input-35-8eb4a3a9fdbf>", line 1, in <module>
    runfile('/home/user/Documents/Code/GniPy/scipy_bug.py', wdir='/home/user/Documents/Code/GniPy')

  File "/usr/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "/usr/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/user/Documents/Code/GniPy/scipy_bug.py", line 43, in <module>
    mm_test = LO_withot_mm.matmat(I)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 326, in matmat
    Y = self._matmat(X)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 468, in _matmat
    return super(_CustomLinearOperator, self)._matmat(X)

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 174, in _matmat
    return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 174, in <listcomp>
    return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])

  File "/usr/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 229, in matvec
    y = y.reshape(M,1)

ValueError: cannot reshape array of size 64 into shape (8,1)

Scipy/Numpy/Python version information:

1.0.0 1.14.0 sys.version_info(major=3, minor=6, micro=4, releaselevel='final', serial=0)