By its anatomy the DGEMM ($C = \alpha AB + \beta C$) is one of the most optimizeable routines in computer science. For historic reason this routine is implemented in FORTRAN or the implementation provides at least an interface which is compatible the old FORTRAN one. If you want to call it from C you have to take care of the function arguments because all of them need to be a pointer (even the scalar values $\alpha$ and $\beta$ or the sizes of the matrices).
In order to avoid this converting-to-FORTRAN problem, one can also use the so called CBLAS interface which is provided by all main BLAS implementations (ATLAS, OpenBLAS, Intel MKL) as well. This interface allows you to pass all function arguments in a C style way. But it has one big trap door. If you use the classical row-oriented C-style way for storing two dimensional data, the CBLAS interface will copy all data in the FORTRAN colum-major format which costs time and memory. So if you use CBLAS try to use the colum-major storage as in Fortran.
The question for the most optimized way is easy to answer: If you have an license of the Intel Compiler Suite or at least of the Intel MKL use this as BLAS library. If not than use OpenBLAS. Both reach more than 90% of the theoretical peak performance of the computer but OpenBLAS has some problems with its performance on multi-socket systems.