Speed up my RGF algorithm

Hi everyone,

I have this function called solve_transport_negf which allows me to solve a 1D transport model with Green's formalism.

I am trying to speed up this routine that takes approximately 15 seconds on my device, to less than 1 second. Do you think it's possible?

I am attaching the whole code below with a profiler as well.

The most critical steps are the vector operations inside the two for loops for calculating g and GR. A smart move would be to use parfoor but each iteration is based on the value of a previous one. I tried to vectorize as much as possible but the final result still takes too long.

It is doable, but I need some help to see how. Thank you in advance

function [E, TEL, TER, gEL, gER] = solve_transport_negf(varargin)
    tollim = 1e-6;
    e = 1.6022e-19;

    mat = varargin{1}{1}{1};
    Nx = mat.Nx;
    dE = mat.dE;
    V = mat.V;
    x = mat.x;
    dx = x(2) - x(1);

    H = calc_discrete_hamiltonian(Nx, mat, V);
    Emin = min(V);
    Emax = max(V) + 0.3;
    E = Emin - 5 * dE:dE:Emax;
    NE = length(E);
    
    TEL = zeros(1, NE);
    TER = zeros(1, NE);
    gEL = zeros(NE, Nx);
    gER = zeros(NE, Nx);

    TL = -H(1, 2);
    TR = -H(Nx - 1, Nx);
    
    for IE = 1:NE
        DL = E(IE) - H(1, 1);
        kL = 1 / dx * acos(-DL / (2 * TL));
        Sigma11R = -TL * exp(1i * kL * dx);
        Gamma11 = 1i * (Sigma11R - conj(Sigma11R));

        DR = E(IE) - H(Nx, Nx);
        kR = 1 / dx * acos(-DR / (2 * TR));
        SigmaNNR = -TR * exp(1i * kR * dx);
        GammaNN = 1i * (SigmaNNR - conj(SigmaNNR));

        SigmaR = sparse(Nx, Nx);
        SigmaR(1, 1) = Sigma11R;
        SigmaR(Nx, Nx) = SigmaNNR;

        g = zeros(1, Nx);
        GR = zeros(Nx, Nx);
        g(Nx) = 1 ./ (E(IE) - H(Nx, Nx) - SigmaR(Nx, Nx));
        
        for I = Nx - 1:-1:2
            g(I) = 1 ./ (E(IE) - H(I, I) - H(I, I + 1) * g(I + 1) * H(I + 1, I));
        end
        
        GR(1, 1) = 1 / (E(IE) - H(1, 1) - SigmaR(1, 1) - H(1, 2) * g(2) * H(2, 1));
        
        for I = 2:Nx
            GR(I, 1) = g(I) * H(I, I - 1) * GR(I - 1, 1);
            GR(I, I) = g(I) + g(I) * (-H(I, I - 1)) * GR(I - 1, I - 1) * (-H(I - 1, I)) * g(I);
        end

        gEL(IE, :) = real(1e9 / (2 * pi * dx) * diag(GR(:, 1) .* Gamma11 * ctranspose(GR(:, 1))));
        gER(IE, :) = real(-1e9 / (2 * pi * dx) * 2 * imag(diag(GR)));
        
        TEL(IE) = GR(Nx, 1) * Gamma11 * conj(GR(Nx, 1)) * GammaNN;
        TER(IE) = GR(1, Nx) * GammaNN * conj(GR(1, Nx)) * Gamma11;
    end
    gER = gER - gEL;
end

[Voss]

Please provide:

• example inputs to the function, saved in a mat-file and uploaded here (use the paperclip button to upload a file)
• an example usage of the function, i.e., show how solve_transport_negf is called in your code
• the definition of any other functions (e.g., calc_discrete_hamiltonian) necessary to run solve_transport_negf

[jo]

Thank you for the answer. Actually, most of the code is quite irrelevant to speed up because it does not take much time and is also part of the WF function which we compare to. The key section is below, which differes from WF and which takes most of the time.

E is a (1,Nx) string, H is a full matrix and SigmaR is a sparse matrix with only (1,1) and (Nx,Nx) elements different from 0

        g = zeros(1, Nx);
        GR = zeros(Nx, Nx);
        g(Nx) = 1 ./ (E(IE) - H(Nx, Nx) - SigmaR(Nx, Nx));
        
        for I = Nx - 1:-1:2
            g(I) = 1 ./ (E(IE) - H(I, I) - H(I, I + 1) * g(I + 1) * H(I + 1, I));
        end
        
        GR(1, 1) = 1 / (E(IE) - H(1, 1) - SigmaR(1, 1) - H(1, 2) * g(2) * H(2, 1));
        
        for I = 2:Nx
            GR(I, 1) = g(I) * H(I, I - 1) * GR(I - 1, 1);
            GR(I, I) = g(I) + g(I) * (-H(I, I - 1)) * GR(I - 1, I - 1) * (-H(I - 1, I)) * g(I);
        end

[Mike Croucher]

@Giovanni you are unlikley to get any help unless you provide everything required for us to run your code on our machines. This is because it is extremely difficult to make code go faster if we cannot run it.

[jo]

I see, please find attached the whole code needed to run the program. In main.m you can switch solver from WF to NEGF. The function that has to be modified is solve_transport_negf