We present a domain decomposition finite element technique for efficiently generating lower and upper bounds to outputs which are linear functionals of the solutions to symmetric or nonsymmetric second-- order elliptic linear partial differential equations in two space dimensions. The method is based upon the construction of an augmented Lagrangian, in which the objective is a quadratic "energy" reformulation of the desired output, and the constraints are the finite element equilibrium equations and intersubdomain continuity requirements. The bounds on the output for a suitably fine "truth--mesh" discretization are then derived by appealing to a dual maxmin relaxation evaluated for optimally chosen adjoint and hybrid--flux candidate Lagrange multipliers generated by a K--element coarser "working--mesh" approximation. Independent of the form of the original partial differential equation, the computation on the truth mesh is reduced to K decoupled subdomain--local, symmetric Neumann pro...