My goal is to minimize a real valued multi-variable function averageFreeEnergyDensity. The function accepts an integer argument n, a set of real parameters j,d,a,h,dx and two real matrices th,ph. I want to find matrices th,ph which minimize the function. FindMinimum minimizes the function well enough but it is too slow at the moment. I am trying to compile the function to accelerate FindMinimum. If I compile the function for arbitrary n it is at least an order of magnitude more slow than if I make a set of functions each of which is defined for a particular n. My current fastest solution is
averageFreeEnergyDensityCompiled["square",n_]:= averageFreeEnergyDensityCompiled["square", n]= Compile[{{j, _Real}, {d, _Real}, {a, _Real}, {h, _Real}, {dx, _Real}, {th, _Real, 2}, {ph, _Real, 2}}, Evaluate[ (*apply boundary conditions*) s[i1_, i2_] := Which[ i1 == n + 1 && i2 == n + 1, s[1, 1], i1 == n + 1, s[1, i2], i2 == n + 1, s[i1, 1], True, {Cos[ph[[i1, i2]]] Sin[th[[i1, i2]]], Sin[ph[[i1, i2]]] Sin[th[[i1, i2]]], Cos[th[[i1, i2]]]} ]; (*define finite difference functions*) dsdx[i1_, i2_] := (s[i1 + 1, i2] - s[i1, i2])/dx; dsdy[i1_, i2_] := (s[i1, i2 + 1] - s[i1, i2])/dx; (*return average free energy density function*) Sum[j (dsdx[i1, i2].dsdx[i1, i2] + dsdy[i1, i2].dsdy[i1, i2]) + a s[i1, i2][[3]]^2 - h s[i1, i2][[3]], {i1, 1, n}, {i2, 1, n}]/ n^2 ]] n. I store the function because I will evaluate it for many values of j,d,a,h,dx. This solution is fast but it throws an error the first time the function is evaluated Block[{n = 100}, AbsoluteTiming[ averageFreeEnergyDensityCompiled["square", n][1, 1, 1, 0, 0.1, Array[0 &, {n, n}], RandomReal[2 Pi, {n, n}]]]] Part::partd: Part specification ph$[[2,1]] is longer than depth of object. >> Part::partd: Part specification th$[[2,1]] is longer than depth of object. >> Part::partd: Part specification ph$[[2,1]] is longer than depth of object. >> General::stop: Further output of Part::partd will be suppressed during this calculation. >> Part function Mathematica gives up and evaluates the expression ignoring Part. The final expression can then be compiled. I want to manually tell Mathematica to ignore Part. I have tried using Inactive,Activate but this slows down the final compiled function by orders of magnitude. I would ideally give an option to Evaluate like Evaluate["insert expression here",Holding->Part] but Evaluate has no options that I know of. Note: The boundary condition I have used here is simple and can be taken out of compile, but I will be using more complicated boundary conditions once I solve this problem. I do not need to use finite difference steps of order greater than 1. I will also be using this code on a triangular grid which will further change the finite difference.
