This question has two parts: first, is the problem in NP, and second, is it NP-hard?
For the first part, I have a positive answer with a non-obvious proof. (Thanks to Suresh for pointing out an earlier error.)
Consider the following way to formalize the question as a decision problem:
UNRESTRICTED MAGIC SQUARE COMPLETION
Input: positive integer $n$ given in unary, list of integers with their positions in an $n$ by $n$ grid
Question: do there exist integers for the remaining positions in the grid, so that the arrangement forms a magic square?
If we add the restriction that each of the integers $1,2,\dots,n^2$ must occur precisely once in the magic square, then the resulting MAGIC SQUARE COMPLETION decision problem is obviously in NP. The definition of a magic square in the 1911 Encyclopædia Britannica, following Euler, has this restriction; in contrast, the Wikipedia article currently uses the terminology "normal magic square" and reserves "magic square" for the unrestricted version.
With an $n$ by $n$ grid, at least $n$ numbers must be given, otherwise the answer is trivially "YES" for the unrestricted version. The size of the input can therefore be assumed to require more than $n$ bits in this case. For the normal version, it is possible that there are inputs that require few bits but which have no solution; to avoid such complications I have specified that $n$ is given in unary.
The argument uses a bound on the possible size of integers that appear in solutions. In the normal case this bound is obviously $n^2$, but in the general case it is not a priori obvious that such a bound exists. It turns out that an exponential bound exists.
Theorem (Tyszka, Theorem 12): Any system of Diophantine equations involving equations of the form $x_i = 1$ and $x_i = x_j + x_k$, for $i,j,k \in \{1,2,\dots,n\}$, either has no integer solution, or has a solution in which every $x_i$ is an integer and at most $\sqrt{5}^{n-1}$ in absolute value.
This also appeared as Theorem 4.7 in:
Cipu recently announced an asymptotically better bound of $2^n$. (Note that the smallest possible bound is $2^{n-1}$.) The argument builds on a bound on the determinant of a matrix, due to Waldi.
Theorem (Cipu, 2011): Any system of Diophantine equations involving equations of the form $x_i = 1$ and $x_i = x_j + x_k$, for $i,j,k \in \{1,2,\dots,n\}$, either has no integer solution, or has a solution in which every $x_i$ is an integer and at most $2^n$ in absolute value.
Even more recently, Freitas, Friedland and Porta showed that the $2^{n-1}$ bound holds, as a corollary of their more general Theorem 1.1.
This yields the following:
Corollary: If an instance of UNRESTRICTED MAGIC SQUARE COMPLETION of size $N$ has a solution, then it has a solution using only numbers up to $2^{O(N^2)}$.
This means that one can use $O(N^4)$ space to guess a solution up to the bound, and then check in $O(N^8)$ time whether it is a solution. The precise polynomial depends on whether one uses Tyszka's or Cipu's result for the bound, and how one expresses the Diophantine system of equations representing the magic square. In either case, the number of variables required in a Diophantine system of the form posed by Tyszka is no more than $n^2 + 2(n+1)(n-2)+1 = 3n^2-2n-3$. This is achieved with $n-2$ variables for partial sums for each row, column, and diagonal and a grand total variable linking these together. Beyond the magic square itself, a further polynomial number of variables is required for the numbers in the square: a number requiring $m$ bits can be represented by using $O(m^2)$ intermediate variables.
For the second part of the question, as far as I can tell either version of MAGIC SQUARE COMPLETION should be NP-hard, but I do not have reductions. It is worth noting that there are procedures to construct normal magic squares of arbitrarily large size; moreover, the number of normal magic squares seems to grow superpolynomially with $n$ (see OEIS A006052) so the underlying language does not seem to be sparse.
Using Papadimitriou's bound on the solutions of an instance of INTEGER LINEAR PROGRAMMING, one can also show that the version where the numbers must all be non-negative is also in NP.
Theorem (Papadimitriou): Let $A$ be an $r \times s$ matrix and $b$ an $r$-vector, both with entries from $\{-a, -a+1, \dots, a-1, a\}$. Then if $Ax = b$ has a solution in non-negative integers, it also has one where every component is in $\{0,1,\dots,s(ra)^{2r+1}\}$.
The constraints forming the magic square can be expressed as a linear program using $a=1$, with $s=n^2+1$ variables and $r=2n+2$ inequalities. This yields a somewhat larger bound than the bound above where negative numbers are allowed, but the certificate is still of polynomial size. (Thanks to Tony Tan for pointing out the result of Papadimitriou.)
- Christos H. Papadimitriou, On the complexity of integer programming, JACM 28 765–768, 1981. (link)