Skip to main content

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"
]];

Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...