### PARAMETER ###
param n;

### VARIABLE ###
var y {1..n, 1..n, 1..n*n} binary;
var constant integer;

### CONSTRAINT ###
subject to OneEach {i in 1..n, j in 1..n}: sum {k in 1..n*n} y[i,j,k] = 1;
subject to AllDifferent {k in 1..n*n}: sum {i in 1..n, j in 1..n} y[i,j,k] <= 1;

s.t. Row1Sum {k in 1..n * n}: sum {j in 1..n} k * y [1, j, k] = constant;
s.t. Row2Sum {k in 1..n * n}: sum {j in 1..n} k * y [2, j, k] = constant;
s.t. Row3Sum {k in 1..n * n}: sum {j in 1..n} k * y [3, j, k] = constant;
s.t. Col1Sum {k in 1..n * n}: sum {i in 1..n} k * y [i, 1, k] = constant;
s.t. Col2Sum {k in 1..n * n}: sum {i in 1..n} k * y [i, 2, k] = constant;
s.t. Col3Sum {k in 1..n * n}: sum {i in 1..n} k * y [i, 3, k] = constant;

s.t. q1: ((k*k)/2)*(k*k+1) = constant*n;

### OBJECTIVE ###
minimize Magic_Constant: constant;

data;
param n:= 3;
option solver cplex;
solve;

display Magic_Constant;

for{i in 1..n}
{
   for {j in 1..n}
   {
   	printf "%3d ", sum {k in 1..n*n} k * y[j,i,k];
   }
    	printf "\n";
}