Can SetDelayed (:=) be compiled?

Question

Is it possible to put an If statement as below within a compile (see below)? I received a warning about SetDelayed.

  eef = Compile[{{μ, _Real}, {NNN, _Integer}}, 
   Module[{NN = NNN}, 
     kd[m_, n_] := If[m == n, 1, 0];
     cosfun[m_, n_] := 
       If[m == n, 
          0, 
          (1 - Cos[(m + n) π])/( (m + n) π) 
          + ((Cos[(m - n) π] - 1)/( (m - n) π))
       ];

      mm = Table[kd[m, n] (MM - B2 k k - B1 (m Pi/L)^2),  
              {m, 1, NN}, {n, 1, NN}];
      kz = Table[cosfun[m, n] (-I A2 m Pi/L), {m, 1, NN}, {n, 1, NN}];
      kxM = Table[kd[m, n] A1 k, {m, 1, NN}, {n, 1, NN}];
      μM = Table[kd[m, n] μ, {m, 1, NN}, {n, 1, NN}];

      HH = ArrayFlatten[{
             {μM + mm, 0 mm, kz, kxM}, 
             {0 mm, μM + mm, kxM, -kz}, 
             {kz, kxM, μM - mm, 0 mm}, 
             {kxM, -kz, 0 mm, μM - mm}
           }];

      ees = Table[Eigenvalues[HH], {k, -.1, .1, .01}]

     ] (* End Module  *)
    ]  (* End Compile *)

where A1, A2, B1, B2, and MM are global variables.

It would help to have some test values for your global variables. Also, $L$ is also a global variable, I think. — Verbeia, Mar 12 '12 at 05:30
You can define cosfun and kd (why redefine KroneckerDelta?) before with compile and use them in your compiled function without problems. — FJRA, Mar 12 '12 at 05:31

score 10 · Accepted Answer · edited Aug 13 '13 at 11:05

Compile can't handle SetDelayed.

In your specific case, you might be able to avoid the need for SetDelayed altogether through the use of Boole.

kd[m_, n_] :=If[m == n, 1, 0] can be replaced by Boole[n==m]. (or as FJRA pointed out in comments, by KroneckerDelta.)

cosfun[m_, n_] := 
 If[m == n, 0, (1 - Cos[(m + n) π])/((m + n) π) + 
  ((Cos[(m - n) π] - 1)/((m - n) π))]

can be replaced by

Boole[m == n] (1 - Cos[(m + n) π])/((m +  n) π) + 
 ((Cos[(m - n) π] - 1)/((m - n) π)) (-I A2 m Pi/L)

Giving

eenew = Compile[{{μ, _Real}, {NNN, _Integer}}, 
  Module[{NN = NNN}, (* do you even need this placeholder for the input? *)
   mm = Table[Boole[m == n] (MM - B2 k k - B1 (m Pi/L)^2), {m, 1, NN}, {n, 1, NN}]; 
   kz =  Table[Boole[m != n] (1 -  Cos[(m + n) π])/((m + n) π) + 
     ((Cos[(m - n) π] - 1)/((m - n) π)) (-I A2 m Pi/L), 
     {m, 1, NN}, {n, 1, NN}]; 
   kxM =   Table[Boole[m == n] A1 k, {m, 1, NN}, {n, 1, NN}]; 
  μM = Table[Boole[m == n] μ, {m, 1, NN}, {n, 1, NN}]; 
   HH =    ArrayFlatten[{{μM + mm, 0 mm, kz, kxM}, 
      {0 mm, μM + mm, kxM, -kz}, {kz, kxM, μM - mm,0 mm}, 
      {kxM, -kz, 0 mm, μM - mm}}]; 
  ees = Table[Eigenvalues[HH], {k, -.1, .1, .01}]]]

What is k for in this last line? There is no way for the iterator to matter for HH.

This doesn't give the error involving SetDelayed anymore.

The Indeterminate expression you were getting had to do with letting $m \to n$ in kz. I corrected your code, so this won't happen. — rcollyer, Mar 12 '12 at 14:31
Note that SetDelayed is included in this list of compilable functions: http://mathematica.stackexchange.com/a/1101/219 — faleichik, May 06 '12 at 19:01

score 8 · Answer 2 · edited Apr 13 '17 at 12:55

Since your matrices are all of the same size, I would construct them in a different manner entirely, and not use SetDelayed at all.

Starting with mm, kxm, and μM, I would note that they are all diagonal, and except for mm, they are all multiples of the identity matrix. So, I would use this instead,

kxm = A1 k IdentityMatrix[NN];
μM = μ IdentityMatrix[NN];

or, even

{kxm, μM} = {A1 k #, μ #}& @ IdentityMatrix[3];

and for mm which depends on the position in the diagonal

mm = DiagonalMatrix[MM - B2 k k - B1 (Range[NN] Pi/L)^2];

If you have not encountered it before, the construction of mm may seem a little odd, but it relies on the fact that some mathematical operations are vectorizable in Mathematica, e.g. q {1, 2, 3, 4, 5}^2 returns

{q, 4 q, 9 q, 25 q}

Also, DiagonalMatrix is a built-in function for constructing just this type of matrix, and it is compilable (as is Range, too).

For kxm, I would embed you If statement directly in Table, as follows

kxM =  
 Table[
  If[m == n, 0, 
   (-2 n + (m + n) Cos[(m - n)Pi] + (-m + n) Cos[(m + n) Pi])/((m^2 - n^2)Pi^2),
  {m, 1, NN}, {n, 1, NN}
 ];

Lastly, ArrayFlatten has the property that you can use scalars in place of constant matrices, provided that there is at least one matrix in both the same row and column that Mathematica can use to determine the size from. For example,

 ArrayFlatten[{{IdentityMatrix[2], 5},{b^2, IdentityMatrix[2]}}]

{{1, 0, 5, 5}, {0, 1, 5, 5}, {b^2, b^2, 1, 0}, {b^2, b^2, 0, 1}}

but this the size of the c block in this matrix can not be determined to be anything but $1 \times 1$,

ArrayFlatten[{{IdentityMatrix[3], 5}, {0, c}}]

{{1, 0, 0, 5}, {0, 1, 0, 5}, {0, 0, 1, 5}, {0, 0, 0, c}}

Using this, I would rewrite your code as

eef = Compile[{{μ, _Real}, {NNN, _Integer}}, 
Module[{}, 

  mm = DiagonalMatrix[MM - B2 k k - B1 (Range[NN] Pi/L)^2];
  {kxm, μM} = {A1 k #, μ #}& @ IdentityMatrix[3];

  kxM =  
   Table[
    If[m == n, 0, 
     (-2 n + (m + n) Cos[(m - n)Pi] + (-m + n) Cos[(m + n) Pi])/((m^2 - n^2)Pi^2),
    {m, NNN}, {n, NNN}
   ];


  HH = ArrayFlatten[{
         {μM + mm, 0, kz, kxM}, 
         {0, μM + mm, kxM, -kz}, 
         {kz, kxM, μM - mm, 0}, 
         {kxM, -kz, 0, μM - mm}
       }];

  ees = Table[Eigenvalues[HH], {k, -.1, .1, .01}]

 ] (* End Module  *)
]  (* End Compile *)

Can SetDelayed (:=) be compiled?

2 Answers2

Linked