Complete the Magic Square
MATL, 36 bytes
`nZ@[]etGg)GXz-yt!hs&amp;ytXdwPXdhsh&amp;-ha

The input is an \$ n \times n\$ matrix, with \$0\$ for the unknown numbers.
The code keeps generating random \$ n \times n\$ matrices formed by the numbers \$1, \dots, n^2\$ until one such matrix meets the required conditions. This procedure is guaranteed to finish with probability one.
This is a terrible approach, in that:

running time is random, and unbounded;
average running time increases as \$(n^2)!\$ (that's more than exponentially);
and so it very likely times out in the online interpreter.

... but hey, it's the shortest answer so far!
(Don't) try it online.
See below a sped-up animated GIF of an example run which happened to take about 2 minutes, here compressed to a few seconds.

Explanation
`         % Do...while
  n       %   Number of elements. This implictly takes the input in the first
          %   iteration, or uses the candidate solution from the previous iteration.
          %   Let this number be denoted as N
  Z@      %   Random permutation of integers 1, 2, ..., N
  []e     %   Reshape as a square matrix. This yields a candidate solution
  t       %   Duplicate
  Gg)     %   Push input, convert to logical, index: this produces a column vector
          %   of the entries of the candidate solution that correspond to nonzero
          %   entries in the input matrix
  GXz     %   Push input, take its nonzero elements. Gives a column vector
  -       %   Element-wise difference (*). This will be all zeros for a valid
          %   solution
  y       %   Duplicate second-top object from the stack, that is, the candidate
          %   solution
  t!      %   Duplicate, transpose
  h       %   Concatenate horizontally
  s       %   Sum of columns. This also gives the sum of rows, thanks to the
          %   concatenated, transposed copy. The result is a two-element row
          %   vector (**)
  &amp;y      %   Duplicate third-top object from the stack: the candidate solution
  tXd     %   Duplicate, extract diagonal as a column vector
  wPXd    %   Swap, flip vertically, extract diagonal. This gives the anti-diagonal
          %   as a column vector
  h       %   Concatenate horizontally
  s       %   Sum of each column. This gives the sum of the diagonal and that  
          %   of the anti-diagonal
  h       %   Concatenate horizontally with (**)
  &amp;-      %   Matrix of all element-wise differences. This will be a matrix of
          %   zeros for a valid solution (***)
  h       %   Concatenate (*) and (***) horizontally. Since sizes do not match,
          %   both (*) and (***) are first linearized to row vectors, and the
          %   result is a row vector
  a       %   Any. This gives true if any element is non-zero
          % End (implicit). A new iteration is run if the top of the stack is true
          % Display (implicit). The candidate solution from the last iteration is
          % the valid solution