matrix - Can any one help me make my program work faster?

In order to construct the matrix h1+h2, I have used three Table commands. I am using fxn[n1a_, n1b_, n2a_, n2b_] as a function in a loop. I need my program to be a lot faster than what it is now, since this is only a part of my main program. Could any one tell me how I can make it work faster? In other words, what are the points I should pay attention to in Mathematica, in case of timing?

Description:

In the original program I am using an iteration method in which I need to construct the matrix h1+h2, calculate the new eigenvectors in each step, reconstruct the matrix h1+h2, until I observe convergence in eigenenergies.

avec, is a guessed vector, and I am using it as an initial vector to construct the matrix h1+h2.

Later on, I need to calculate the eigenvectors of this matrix, choose the vectors related to the minimum eigenvalues and put them instead of the previous avec.

h1 is the kinetic term of my matrix and h2 is the exchange potential (the formula which I have used came from the math I did for Hartree-Fock approximation).

fxn is a function I constructed based on the pre-calculated integrals in my code (so that I do not need to calculate the integrals in each iteration).

The problem now is the calculation for h2, which consumes most of the time. I assume the problem is the way I am using the Table commands. But I couldn't think of a more efficient way.

  avec = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
  nvec = {1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6};
  svec = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
  kdfxn[i_, j_] := If[i == j, 1, 0]

  ne = 3;
  n\[Mu] = 12;

  \[Delta] = -50;
  \[Beta] = 1;

   fxn[n1a_, n1b_, n2a_, n2b_] := N[ Which[
   n1a == n1b && n2a == n2b,
   1/6 (1 - 3/(n1a^2 \[Pi]^2) - 3/(n2a^2 \[Pi]^2)),
   n1a == n1b && n2a != n2b,
   (4 (1 + (-1)^(n2a + n2b)) n2a n2b)/((n2a^2 - n2b^2)^2 \[Pi]^2),
   n1a != n1b && n2a == n2b,
   (4 (1 + (-1)^(n1a + n1b)) n1a n1b)/((n1a^2 - n1b^2)^2 \[Pi]^2),

   True, -((
    32 (-1 + (-1)^(n1a + n1b)) (-1 + (-1)^(
   n2a + n2b)) n1a n1b n2a n2b)/((n1a^2 - n1b^2)^2 (n2a^2 - 
   n2b^2)^2 \[Pi]^4))]]

hmat = Table[
n0 = nvec[[n\[Mu]0]];
n1 = nvec[[n\[Mu]1]];
h1 = (n0^2 \[Pi]^2 + \[Beta] svec[[n\[Mu]0]]) kdfxn[n\[Mu]0, 
  n\[Mu]1] ;

h2 = Total[Table[
  n2 = nvec[[n\[Mu]2]];
  n3 = nvec[[n\[Mu]3]];
  Table[
   temp1 = 
    fxn[n1, n0, n3, n2] avec[[j, n\[Mu]3]] avec[[j, n\[Mu]2]] ;
   temp2 = 
    fxn[n1, n3, n0, n2] avec[[j, n\[Mu]3]] avec[[j, n\[Mu]2]];
   temp = \[Delta] (temp1 - temp2), {j, 1, ne}],
  {n\[Mu]2, 1, n\[Mu]},

  {n\[Mu]3, 1, n\[Mu]}], Infinity];
  h1 + h2
 , {n\[Mu]0, 1, n\[Mu]}, {n\[Mu]1, 1, n\[Mu]}] // Timing

---

I want to use the code @Henrik Schumacher gave me to find my new avecs (which are coming from the eigenvectors of hcmat) to recalculate hcmat and iterate. The problem is I can do it manually by substituting new avec instead of the previous one, but as I start using table it doesn't work. I don't know why it's acting like that.

Table[chmat = 
    With[{ccfxn = cfxn}, 
    Compile[{{nm, _Integer}, {ne, _Integer}, {d, _Real}, {b, \
   _Real}, {avec, _Real, 2}, {nvec, _Real, 1}, {svec, _Real, 1}}, 

        Block[{n0, n1, n2, n3, h1, h2, tmp, tmp2}, 
        Table[n0 = Compile`GetElement[nvec, nm0];
        n1 = Compile`GetElement[nvec, nm1];
        tmp = 0.;
        Do[n2 = Compile`GetElement[nvec, nm2];
        n3 = Compile`GetElement[nvec, nm3];
        tmp2 = (ccfxn[n1, n0, n3, n2] - ccfxn[n1, n3, n0, n2]);

        Do[tmp += 
        tmp2 Compile`GetElement[avec, j, nm3] Compile`GetElement[

           avec, j, nm2], {j, 1, ne}], {nm2, 1, nm}, {nm3, 1, 
        nm}];
        d tmp + 
        If[nm0 == 
          nm1, (n0^2 Pi^2 + b Compile`GetElement[svec, nm0]), 
          0.], {nm0, 1, nm}, {nm1, 1, nm}]], 
       CompilationTarget -> "C", 
       CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
       RuntimeOptions -> "Speed"]];


       hmat2 = chmat[n\[Mu], ne, \[Delta], \[Beta], avec, nvec, svec]; //
       AbsoluteTiming // First
            hmat = hmat2

       eout = Eigensystem[hmat];
       evals = eout[[1]];
       evecs = eout[[2]];
       sortedevals = Sort[evals];
       bvec = Chop[Table[bpos = Position[evals, sortedevals[[i]]][[1, 1]];
       evecs[[bpos]], {i, 1, ne}]];

       Table[If[Sign[bvec[[i, 1]]] == Sign[avec[[i, 1]]], avec = ( bvec) , 
       avec = (Sign[avec[[i, 1]]] bvec) ], {i, 1, ne}];
       normalizedavec = Table[avec[[i]]/norm[avec, i], {i, 1, ne}];
       avec = normalizedavec;
       en = Total[Take[sortedevals, ne]], {j, 1, 5}] 

In the rest of the code, I am finding the first ne(3) minimum eigenenergies and finding the the related eigenvectors, to construct bvec. by doing some modification on bvec and normalizing it ,I could get the new avec I want to use instead of the previous one. Now my problem is doing the iteration. I would appreciate it if you could help me with this.

I want to calculate this integrals

 sol = Assuming[
 A \[Element] Reals && 0 <= s1 <= 1 && 1 <= n1a <= 10 && 

 1 <= n2a <= 10 && 1 <= n1b <= 10 && 1 <= n2b <= 10, 
 Expand@FullSimplify[
 Integrate[
  4 A Exp[-B (s2 - s1)] Sin[n1a \[Pi] s1] Sin[n1b \[Pi] s1] Sin[
    n2a \[Pi] s2] Sin[n2b \[Pi] s2], {s2, s1, 1}]]];

 sol1 = Assuming[
  A \[Element] Reals && 0 <= s2 <= 1 && 1 <= n1a <= 10 && 
  1 <= n2a <= 10 && 1 <= n1b <= 10 && 1 <= n2b <= 10, 
  Expand@FullSimplify[

  Integrate[
   4 A Exp[-B (s1 - s2)] Sin[n1a \[Pi] s1] Sin[n1b \[Pi] s1] Sin[
    n2a \[Pi] s2] Sin[n2b \[Pi] s2], {s1, s2, 1}]]];

(int[n1a_, n2a_, n1b_, n2b_] = 
 Integrate[sol, {s1, 0, 1}, 
  Assumptions -> 
  B \[Element] Reals && A \[Element] Reals && 1 <= n1a <= 10 && 
  1 <= n2a <= 10 && 1 <= n1b <= 10 && 1 <= n2b <= 10] + 
  Integrate[sol1, {s2, 0, 1}, 

  Assumptions -> 
  B \[Element] Reals && A \[Element] Reals && 1 <= n1a <= 10 && 
  1 <= n2a <= 10 && 1 <= n1b <= 10 && 1 <= n2b <= 10])

and replace n1a in the int function with n1a+0.0001 . Then I need to use the result of int after the replacement, as the material for code :code=result in:

        "cfxn = Block[{n1a, n1b, n2a, n2b}, 
         With[{code = "result"}, 
         Compile[{{n1a, _Integer}, {n1b, _Integer}, {n2a, _Integer}, {n2b, \
          _Integer}}, code, CompilationTarget -> "C"]]];" 

but I am getting these errors:

CompiledFunction::cfse: Compiled expression CompiledFunction[{10,11.,5464},<<7>>,LibraryFunction[/Users/delaramnematollahi/Library/Mathematica/ApplicationData/CCompilerDriver/BuildFolder/Delarams-MacBook-Pro-3755/compiledFunction3.dylib,compiledFunction3,{<<1>>},{Real,2}]] should be a machine-size real number.

which I do not understand why.

Could you please help me with this?

Answer

Compile can be your friend. Fortunately, your code is almost perfectly prodedural and thus can be compiled with only few changes.

avec = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
nvec = {1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6};

svec = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
ne = 3;
nμ = 12;
δ = -50;
β = 1;

cfxn = Block[{n1a, n1b, n2a, n2b},
   With[{code = Which[
       n1a == n1b && n2a == n2b,
       Evaluate[N[1/6 (1 - 3/(n1a^2 π^2) - 3/(n2a^2 π^2))]],


       n1a == n1b && n2a != n2b,
       Evaluate[N[(4 (1 + (-1)^(n2a + n2b)) n2a n2b)/((n2a^2 - n2b^2)^2 π^2)]],

       n1a != n1b && n2a == n2b,
       Evaluate[N[(4 (1 + (-1)^(n1a + n1b)) n1a n1b)/((n1a^2 - n1b^2)^2 π^2)]],       

       True, 
       Evaluate[N[-((32 (-1 + (-1)^(n1a + n1b)) (-1 + (-1)^(n2a + n2b)) n1a n1b n2a n2b)/((n1a^2 - n1b^2)^2 (n2a^2 - n2b^2)^2 π^4))]]
       ]},

    Compile[{{n1a, _Integer}, {n1b, _Integer}, {n2a, _Integer}, {n2b, _Integer}},
     code,
     CompilationTarget -> "C"
     ]
    ]];

Here I tried to enforce that cfxn uses floating point numbers internally (therefore the Evaluate) such that less Integers have to be type cast to Reals.

chmat = With[{ccfxn = cfxn},
   Compile[{{nm, _Integer}, {ne, _Integer}, {d, _Real}, {b, _Real}, {avec, _Real, 2}, {nvec, _Real, 1}, {svec, _Real, 1}},
    Block[{n0, n1, n2, n3, h1, h2, tmp, tmp2},

     Table[
      n0 = nvec[[nm0]];
      n1 = nvec[[nm1]];
      tmp = 0.;
      Do[
       n2 = nvec[[nm2]];
       n3 = nvec[[nm3]];
       tmp2 = (ccfxn[n1, n0, n3, n2] - ccfxn[n1, n3, n0, n2]);
       Do[tmp += tmp2 avec[[j, nm3]] avec[[j, nm2]], {j, 1, ne}],
        {nm2, 1, nm}, {nm3, 1, nm}];

      d tmp + If[nm0 == nm1, N[(n0^2 Pi^2 + b svec[[nm0]])], 0.],
      {nm0, 1, nm}, {nm1, 1, nm}]
     ],
    CompilationTarget -> "C",
    CompilationOptions -> {"InlineCompiledFunctions" -> True}
    ]];

Edit: My former implementation of chmat did not produce the correct result. I removed it and also optimized away the Total: Summing up the results in a fixed Real register has the advantage that no auxiliary arrays have to be allocated and that no indexing into such arrays takes place within Total.

You can compute hmat with

hmat2= chmat[nμ, ne, δ, β, avec, nvec, svec];//AbsoluteTiming//First

hmat==hmat2

(* 0.003411 *)
(* True *)

On my machine, this results in a 600-fold speed-up.

Addendum:

This need not to be the end of the story. By deactivating some security checks and replacing Part with CompileGetElement` I even get a speed-up factor of over 2500 (compared to the original code). This is how it looks:

chmat = With[{ccfxn = cfxn}, 
   Compile[{{nm, _Integer}, {ne, _Integer}, {d, _Real}, {b, _Real},  {avec, _Real, 2}, {nvec, _Real, 1}, {svec, _Real, 1}},

    Block[{n0, n1, n2, n3, h1, h2, tmp, tmp2},
     Table[
      n0 = Compile`GetElement[nvec, nm0];
      n1 = Compile`GetElement[nvec, nm1];
      tmp = 0.;
      Do[
       n2 = Compile`GetElement[nvec, nm2];
       n3 = Compile`GetElement[nvec, nm3];
       tmp2 = (ccfxn[n1, n0, n3, n2] - ccfxn[n1, n3, n0, n2]);
       Do[

         tmp += tmp2 Compile`GetElement[avec, j, nm3] Compile`GetElement[avec, j, nm2],
        {j, 1, ne}],
       {nm2, 1, nm}, {nm3, 1, nm}];
      d tmp + If[nm0 == nm1, (n0^2 Pi^2 + b Compile`GetElement[svec, nm0]), 0.],
      {nm0, 1, nm}, {nm1, 1, nm}]
     ],
    CompilationTarget -> "C",
    CompilationOptions -> {"InlineCompiledFunctions" -> True},
    RuntimeOptions -> "Speed"
    ]];