How Can I Vectorize a Computation in a Python Device Simulation Algorithm?

[maverick280857]

Hi,

I am working on a device simulation algorithm, and am implementing it in Python on my laptop that runs Arch Linux. A particular step requires me to perform a computation of the form

$$\tilde{G} = G \Sigma_{c}^{in} G^{\dagger}$$

where $G$ and $\Sigma_{c}^{in}$ are both $N \times N$ matrices, and $\Sigma_{c}^{in}$ is in fact a diagonal matrix.

N is of the order of 5040

Now, if I do this naively, the command to do it in Python is

Gtilde_n = Gtilde_n + dot( dot(G, sigcIn), conjugate(transpose(G)) )

(The addition is for purposes of sub-band summation, and isn't relevant to the question I'm posing.)

but this takes about 2 hours when it works (sometimes it just crashes my laptop).

Now if I utilize the fact that $\Sigma_{c}^{in}$ is diagonal, then I can do something more naive and write a for loop of the form

              for row in range(0,N):
                      for col in range(0,N):
                              sum = 0
                              for alpha in range(0,N):
                                      sum = sum + G[row,alpha]*sigcIn[alpha,alpha]*conjugate(G[col,alpha])
                              Gtilde_n[row,col] = Gtilde_n[row,col] + sum

But this doesn't look like a nice piece of code: it has 3 nested for loops, and looks a rather ancient way of doing the following calculation

$$\tilde{G}_{ij} = \sum_{\alpha}\sum_{\beta} G_{i\alpha}(\Sigma_{c}^{in})_{\alpha\beta}(G^{\dagger})_{\beta j} = \sum_{\alpha,\beta} G_{i\alpha}(\Sigma_{c}^{in})_{\alpha\alpha}\delta_{\alpha\beta} (G^{\dagger})_{\beta j} =\sum_{\alpha} G_{i\alpha}(\Sigma_{c}^{in})_{\alpha\alpha}G^{*}_{j\alpha}$$

Is there some way to vectorize this?

In general, is there a better way to perform this computation?

 

[Solkar]

$$G_{ij} = \cdots = \sum_{\alpha} G_{i\alpha}(\Sigma_{c}^{in})_{\alpha\alpha}G^{*}_{j\alpha}$$

Is there some way to vectorize this?

Trivially by partioning α's range for each ij, but I would rather try to exploit this
$$G_{i\alpha} G^{*}_{j\alpha} = \left(G_{j\alpha} G^{*}_{i\alpha}\right)^{*}$$ or/and(given Σ is real, what I suspect) premultiplying all $(G_{i})_{\alpha}$ vectors component-wise by a vector made of the Σ diag components. That could well be done in parallel.

In general, is there a better way to perform this computation?

Are there any known algebraic constraints on $G$?

 

[maverick280857]

Thanks for your reply.

Are there any known algebraic constraints on $G$?

No, its just the matrix representation of a Green's function.

 

[Mark44]

What does this notation mean - $\Sigma_{c}^{in}$?

You said that it was a diagonal matrix, but I'm curious about the significance of c and in.

 

[Solkar]

No, its just the matrix representation of a Green's function.

You're sure there's no constraint related to its eigenspaces?

 

[maverick280857]

What does this notation mean - $\Sigma_{c}^{in}$?

You said that it was a diagonal matrix, but I'm curious about the significance of c and in.

So $\Sigma_{c}^{in} = \Sigma_{1}^{in} + \Sigma_{2}^{in}$ is simply the sum of two so called inscattering matrices, computed from the self-energy in a many body problem. There are two terms because there are two contacts. For each contact,

$$\Sigma_{k}^{in} = f(E)\Gamma_k$$

k = 1, 2, and

$$\Gamma_k = i[\Sigma_k - \Sigma_k^\dagger]$$

It turns out that [tex]\Gamma_i[/itex] is a real matrix, i.e. all its elements are real.