Contracting tensors

Given $F^{\mu\nu}$ defined as in QED: $$F^{\mu\nu} = \begin{bmatrix} 0&E_x&E_y&E_z\\ -E_x&0&-B_z&B_y\\ -E_y&B_z&0&-B_x\\ -E_z&-B_y&B_x&0 \end{bmatrix}$$ when contracting it in the lagrangian as $F^{\mu\nu}F_{\mu\nu}$ how do I compute it? Simply as a product of matrices?

Your first equation involves a common, but misleading, abuse of notation. The left-hand side shows a generic contravariant component $F^{\mu\nu}$ of the field strength tensor and the right hand side shows all components. It would be better to write something like $[F^{\mu\nu}]$ on the left-hand side (simply writing $F$ would be ambiguous as to whether the components are co- or contravariant).
Now, it should be absolutely clear, how to compute $F^{\mu\nu}F_{\mu\nu}$. The Einstein sum convention tells you to sum over all indices that occur once as upper and once as lower index. In other words one has to compute $\sum_\mu \sum_\nu F^{\mu\nu}F_{\mu\nu}$ (here without implicit summation!), where $F^{\mu\nu}$ and $F_{\mu\nu}$ are simply numbers. Don't forget, that $F^{\mu\nu} \ne F_{\mu\nu}$, but you have to contract with the metric tensor first to lower the indices (again using the sum convention): $$F_{\mu\nu} = F^{\alpha\beta}g_{\mu\alpha}g_{\nu\beta}.$$