Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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 5: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 5:31
add comment

2 Answers 2

up vote 7 down vote accepted

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.

share|improve this answer
1  
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: mathematica.stackexchange.com/a/1101/219 –  faleichik May 6 '12 at 19:01
add comment

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 \[Mu]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];
\[Mu]M = \[Mu] IdentityMatrix[NN];

or, even

{kxm, \[Mu]M} = {A1 k #, \[Mu] #}& @ 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[{{\[Mu], _Real}, {NNN, _Integer}}, 
Module[{}, 

  mm = DiagonalMatrix[MM - B2 k k - B1 (Range[NN] Pi/L)^2];
  {kxm, \[Mu]M} = {A1 k #, \[Mu] #}& @ 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[{
         {\[Mu]M + mm, 0, kz, kxM}, 
         {0, \[Mu]M + mm, kxM, -kz}, 
         {kz, kxM, \[Mu]M - mm, 0}, 
         {kxM, -kz, 0, \[Mu]M - mm}
       }];

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

 ] (* End Module  *)
]  (* End Compile *)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.