Performance Tuning: Construction of Matrix with Summations

Question

I am trying to solve an eigenvalue problem of a large matrix. The issue is that it takes too much time to construct this large matrix. This is the code that I built:

msa = 1.711*10^6;
msb = 0;
mu0 = 1.256*10^-6;
h0 = 0.1/mu0;
aa = 8.3*10^-12;
ab = 0;
alpa = 0.0019;
alpb = 0;
a = 10*10^-9;
rcyl = Sqrt[a^2/(2*Pi)];

nmax = 10;
n2 = 4;
ncap = (2*n2 + 1)^2;

gx[n_] := 2*n*Pi/a; 
gy[m_] := 2*m*Pi/a;
gxy = Table[{gx[n], gy[m]}, {n, -n2, n2}, {m, -n2, n2}];
g = ArrayReshape[gxy, {(2*n2 + 1)^2, 2}];

msg[{x_, y_}] := If[x == 0 && y == 0, (msa + msb)/2, (msa - msb)* BesselJ[1, (Sqrt[x^2 + y^2]*rcyl)]/(Sqrt[x^2 + y^2]*rcyl)];

q[{x_, y_}] := If[x == 0 && y == 0, ((2*aa/(mu0*h0*msa)))/2, ((2*aa/(mu0*h0*msa)))* BesselJ[1, (Sqrt[x^2 + y^2]*rcyl)]/(Sqrt[x^2 + y^2]*rcyl)];

sum1[kx_, ky_, i_, j_] := Sum[msg[g[[i]] - g[[l]]]* q[g[[l]] - g[[j]]]*(({kx, ky} + g[[j]]).({kx, ky} + g[[l]]) - (g[[i]] - g[[j]]).(g[[i]] - g[[l]])), {l, (2*n2 + 1)^2}];

bxy[kx_, ky_] := Table[If[(g[[j]] + {kx, ky}) == {0, 0}, KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]/(2*h0) + sum1[kx, ky, i, j], KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]*((g[[j, 2]] + ky)*(g[[j, 2]] + ky))/(h0*((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky}))) + sum1[kx, ky, i, j]], {i, 1, (2*n2 + 1)^2}, {j, 1, (2*n2 + 1)^2}];

Next, I used arbitrary kx and ky to construct this bxy matrix. However, because of sum1, it takes a lot of time to construct it. When I used $n2=4$ and $kx=ky=\dfrac{2\pi}{3a}$,

k1=Pi*2/(3*a);
k2=Pi*2/(3*a);

Absolutetiming[bxy[k1,k2]]
(*{501.605,{{9.55609*10^8,4.84554*10^8...}}*)

It took about 502 seconds. It may not look bad, but later, I want to use $n2$ more than 10. Considering that the time would increase exponentially, I cannot even imagine how long it will take with $n2=10$.

After some investigation, I realized that it is caused by sum1 function.

bxy[kx_, ky_] := Table[If[(g[[j]] + {kx, ky}) == {0, 0}, KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]/(2*h0), KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]*((g[[j, 2]] + ky)*(g[[j, 2]] + ky))/(h0*((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky})))], {i, 1, (2*n2 + 1)^2}, {j, 1, (2*n2 + 1)^2}];

k1 = Pi*2/((3)*a);
k2 = Pi*2/((3)*a);

AbsoluteTiming[bxy[k1, k2]]
(*{2.29623,{{6.37254,1.62134...}}}*)

So without sum1 function, it only takes about 2 seconds. I wonder if there is any way to speed up the summation. Thank you in advance!

P.S. I am not really familiar with this forum, so if I did not follow any rules, please let me know.

Answer 1

At the moment, a lot of the computations are performed in exact arithmetic although you want finite precision results in the end. Making this tiny change

 a = N[10 10^-9];

along with defining gby

 g = Tuples[{gx /@ Range[-n2, n2], gy /@ Range[-n2, n2]}];

should speed up computations a lot.

Moreover, you compute many quantities involving the Bessel function multiple times. The Bessel function is rather expensive, so one should its evaluation to a minimum. In the code below, I do that by setting the matrices msgmatrix and qmatrix.

Moreover, you can use replace the Table in bxy by ParallelTable.

Here is the modified code; it executes in about 1.32 seconds on my Quad Core machine:

msa = 1.711 10^6;
msb = 0.;
mu0 = 1.256 10^-6;
h0 = 0.1/mu0;
aa = 8.3 10^-12;
alpa = 0.0019;
a = N[10 10^-9];
rcyl = Sqrt[a^2/(2 Pi)];

nmax = 10;
n2 = 4;
ncap = (2 n2 + 1)^2;

gx[n_] := 2  Pi/a n;
gy[m_] := 2 Pi/a m;
g = Tuples[{gx /@ Range[-n2, n2], gy /@ Range[-n2, n2]}];

msg[{x_, y_}] := If[Sqrt[x^2 + y^2] < 1. 10^-8,
   (msa + msb)/2,
   (msa - msb) BesselJ[
      1, (Sqrt[x^2 + y^2] rcyl)]/(Sqrt[x^2 + y^2] rcyl)
   ];

q[{x_, y_}] := If[Sqrt[x^2 + y^2] < 1. 10^-8,
   (2 aa/(mu0 h0 msa))/2,
   (2 aa/(mu0 h0 msa)) BesselJ[
      1, (Sqrt[x^2 + y^2] rcyl)]/(Sqrt[x^2 + y^2] rcyl)
   ];
msgmatrix = Table[msg[gi - gj], {gi, g}, {gj, g}];
qmatrix = Table[q[gi - gj], {gi, g}, {gj, g}];
sum1[kx_, ky_, i_, j_] := Sum[
   msgmatrix[[i, l]] msgmatrix[[l, j]] (({kx, ky} + g[[j]]).({kx, ky} + g[[l]]) - (g[[i]] - g[[j]]).(g[[i]] - g[[l]])),
   {l, ncap}];

bxy[kx_, ky_] := ParallelTable[
   If[Norm[(g[[j]] + {kx, ky})] < 10^-8,
    KroneckerDelta[i, j] + msgmatrix[[i, j]]/(2 h0) + sum1[kx, ky, i, j], 
    KroneckerDelta[i, j] + msgmatrix[[i, j]] (g[[j, 2]] + ky)^2/(h0 ((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky}))) + sum1[kx, ky, i, j]],
   {i, 1, ncap}, {j, 1, ncap}];

k1 = N[Pi 2/(3*a)];
k2 = N[Pi 2/(3*a)];
A = bxy[k1, k2];

Of course, a lot more can be done. For example, sum1 or bxy could be compiled with Compile...