1
$\begingroup$

Im trying to compute a triple cross product of vectors a,b, and c in real space and integrate over the entire space. The result is a term in the hamiltonian for an electronic system so there are indices for the band number too (you can think of a and b as gradient of wavefunctions in real space and of c as a vector field): $$ h_{nm} = \int \boldsymbol{a}_n(x)\times\boldsymbol{c}(x)\times\boldsymbol{b}_m(x)\,d x %![equation](https://i.stack.imgur.com/VeGe3.gif) $$ where $n$ and $m$ are indices from i.e. 0 to 50. I can write the triple product as $$ \boldsymbol{A}\times\boldsymbol{B}\times\boldsymbol{C} = \boldsymbol{B}(\boldsymbol{A}\cdot\boldsymbol{C})-\boldsymbol{C}(\boldsymbol{A}\cdot\boldsymbol{B}) %![eq](https://i.stack.imgur.com/aiEle.gif) $$ which should be faster.

So for the first term in python I'm doing

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((N,v,g), dtype = 'complex')
b = np.ones((N,v,g), dtype = 'complex')
c = np.ones((v,g), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):
         H_vnn += np.array([c[v,i]*(np.dot(a[:,:,i],b[:,:,i].T)) for v in range(3)])

which, when doesnt run into memory error, takes forever. Im pretty sure this is not the most efficient way of doing it, but the only alternative I can think of is looping over the g index which seems to be even more expansive. Thank to everybody and particular to those who can help me

|cite|improve this question
$\endgroup$
  • $\begingroup$ where does the integration come in? It looks you're just summing over $n$ and $m$. So should your first equation not just be a double sum? $\endgroup$ – Walter Jun 2 '17 at 12:39
  • $\begingroup$ Ops you are right; i have written the wrong code. Instead of the double for loop, is one single loop: for i in range(g): $\endgroup$ – al_j Jun 2 '17 at 12:42
  • 1
    $\begingroup$ It may be worth taking a look at numpy.einsum. (I do not have time for more than this hint or checking whether it actually applies here.) $\endgroup$ – Wrzlprmft Jun 2 '17 at 13:09
  • $\begingroup$ In your code, the sub-expression (np.dot(a[:,:,i],b[:,:,i].T) does not depend on v, so could be pre-computed. And, of course (if you use C-ordering of indices), the index i is the fastest running index in a[:,:,i] or b[:,:,i], so your dot product is over the slower indices, implying extremely inefficient memory access (cache misses all the time). Thus suggests to re-order the indices for a and b as in a = np.ones((g,N,v), dtype = 'complex') (and change the code accordingly). $\endgroup$ – Walter Jun 5 '17 at 8:14
0
$\begingroup$

Try this

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((g,N,v), dtype = 'complex')
b = np.ones((g,N,v), dtype = 'complex')
c = np.ones((g,v), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):
    adotb  = np.dot(a[i,:,:],b[i,:,:].T)
    H_vnn += np.array([c[i,v]*adotb) for v in range(3)])
|cite|improve this answer
$\endgroup$
  • $\begingroup$ It improves speed indeed, not as much as I expected but thank you very much! $\endgroup$ – al_j Jun 6 '17 at 8:57
  • $\begingroup$ How much then? (btw, I appreciate an up-vote or acceptance tick) Also, are you sure the .T does actually what you want? $\endgroup$ – Walter Jun 6 '17 at 17:34
  • $\begingroup$ At least I don't get a memory error anymore so thank you a lot! I am now focusing on optimizing the problem I am working on so to not have all this massive workload, if i find some good hints I will write them here. Anyway, yes, this is what I was looking for $\endgroup$ – al_j Jun 8 '17 at 16:13

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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