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 \Gamma_i[/itex] is a real matrix, i.e. all its elements are real.