def METHOD_NAME(out, *fargs):
    # Compute M_C.
    if diff_var is None:
        sh = out.shape
        _out = nm.zeros((sh[0], sh[1], sh[2], sh[2]), dtype=out.dtype)
    else:
        _out = out
    status = _fun(_out, *fargs)
    # Consistent element mass matrices M_C.
    mc = _out if beta in (0.0, 1.0) else _out.copy()
    # Diagonals of lumped element mass matrices M_L.
    if lumping == 'row_sum':
        mld = mc.sum(-1)
    elif lumping == 'hrz':
        n_bf = mc.shape[-1] // self.region.dim
        det = fargs[2][0][0]
        jrho = det * rho[..., 0, 0]
        diag = nm.einsum('...ii->...i', mc)
        intrho = jrho.sum(axis=1) # Element masses to preserve.
        intdiag = diag[:, 0, :n_bf].sum(axis=1)
        alpha = intrho / intdiag
        mld = alpha[:, None, None] * diag
    else:
        raise ValueError('unknown lumping mode! (%s)' % lumping)
    if beta == 1.0:
        # M_A = M_L.
        eye = nm.eye(mld.shape[-1])
        _out[:] = nm.einsum('can,nm->canm', mld, eye)
    elif beta > 0.0:
        # M_A = (1 - beta) * M_C + beta * M_L.
        _out[:] = (1.0 - beta) * mc
        outd = nm.einsum('...ii->...i', _out) # View to diagonal.
        outd += beta * mld
    else:
        # beta == 0.0, so M_A = M_C in RMM term_mode -> do nothing.
        pass
    if term_mode == 'RMM':
        # out contains M_A, mld contains A_e diagonal.
        ime = nm.linalg.inv(_out)
        # ce = nm.einsum('cik,ckl,clj->cij', ae, ime, ae)
        _out[:] =  mld[..., None] * ime * mld[:, None, :]
    if diff_var is None:
        earg = self.einfos[None].eargs[2]
        step_cache = earg.arg.evaluate_cache.setdefault('dofs', {})
        cache = step_cache.setdefault(0, {})
        dofs = earg.get_dofs(cache, step_cache, state.name + '.dofs')
        dofs = dofs.reshape((sh[0], sh[2])) # Does a copy!
        # Could be faster by exploiting that mld is a vector.
        out[..., 0] = nm.einsum('caij,cj->cai', _out, dofs)
    return status