  We present a general framework to compute upper and lower bounds for
linear-functional outputs of the exact solutions of the Poisson equation based
on reconstructions of the field variable and flux for both the primal and
adjoint problems. The method is devised from a generalization of the
complementary energy principle and the duality theory. Using duality theory,
the computation of bounds is reduced to finding independent potential and
equilibrated flux reconstructions. A generalization of this result is also
introduced, allowing to derive alternative guaranteed bounds from
nearly-arbitrary H(div;{\Omega}) flux reconstructions (only zero-order
equilibration is required). This approach is applicable to any numerical method
used to compute the solution. In this work, the proposed approach is applied to
derive bounds for the hybridizable discontinuous Galerkin (HDG) method. An
attractive feature of the proposed approach is that superconvergence on the
bound gap is achieved, yielding accurate bounds even for very coarse meshes.
Numerical experiments are presented to illustrate the performance and
convergence of the bounds for the HDG method in both uniform and adaptive mesh
refinements.
