'Default matrix-matrix multiplication handler.
Falls back on the user-defined _matvec method, so defining that will
define matrix multiplication (though in a very suboptimal way).'

def _matmat(self, X):

    return np.hstack([self.matvec(col.reshape((-1), 1)) for col in X.T])
