Wolfram Language (Mathematica), 652 bytes
f[m_, n_] := Total@LinearProgramming[Table[-1, m*n], SparseArray[ Flatten@MapIndexed[{#2[[1]], #} -> 1 &, Flatten[{Table[ m*i + j + k + 1, {i, 0, n - 1}, {j, 0, m - 3}, {k, 0, 2}], Table[ m*(i + k) + j + 1, {i, 0, n - 3}, {j, 0, m - 1}, {k, 0, 2}], Table[ m*(i + k) + j + k + 1, {i, 0, n - 3}, {j, 0, m - 3}, {k, 0, 2}], Table[ m*(i + k) + j - k + 1, {i, 0, n - 3}, {j, 2, m - 1}, {k, 0, 2}]}, 2], {2}], {(m - 2) n + m (n - 2) + 2 (m - 2) (n - 2), m*n}], Table[{2, -1}, (m - 2) n + m (n - 2) + 2 (m - 2) (n - 2)], Table[{0, 1}, m*n], Integers] Using the built-in LinearProgramming. Gives the results from 1x1 to 9x9 in 40 seconds on TIO.