Qiskit CircuitStateFn.eval() on more than 16 qubits

Question

The definition of Qiskit CircuitStateFn.eval() is as follows:

def eval(self,
         front: Union[str, dict, np.ndarray,
                      OperatorBase] = None) -> Union[OperatorBase, float, complex]:
    if not self.is_measurement and isinstance(front, OperatorBase):
        raise ValueError(
            'Cannot compute overlap with StateFn or Operator if not Measurement. Try taking '
            'sf.adjoint() first to convert to measurement.')

    # pylint: disable=import-outside-toplevel
    from ..list_ops.list_op import ListOp
    from ..primitive_ops.pauli_op import PauliOp
    from ..primitive_ops.matrix_op import MatrixOp
    from ..primitive_ops.circuit_op import CircuitOp

    if isinstance(front, ListOp) and front.distributive:
        return front.combo_fn([self.eval(front.coeff * front_elem)
                               for front_elem in front.oplist])

    # Composable with circuit
    if isinstance(front, (PauliOp, CircuitOp, MatrixOp, CircuitStateFn)):
        new_front = self.compose(front)
        return new_front.eval()

    return self.to_matrix_op().eval(front)

In the last line to_matrix_op() is called without any arguments. However, in order to work properly for more than 16 qubits, it has to be provided with an additional argument massive = True:

def to_matrix(self, massive: bool = False) -> np.ndarray:
    if self.num_qubits > 16 and not massive:
        raise ValueError(
            'to_vector will return an exponentially large vector, in this case {0} elements.'
            ' Set massive=True if you want to proceed.'.format(2 ** self.num_qubits))

So how do I use CircuitStateFn.eval() for more than 16 qubits??