Code that creates a scale-free (preferential attachment) graph edge list - Is there a better or more efficient way to code this?

Hello, I am trying to create a scale-free graph (network) edge list. By edge list, I mean M-by-2 matrix with first column representing originating node and the second column representing destination node. I am working with undirected graph, so there is no distinction between origination or destination nodes. For example, [1 2;5 11], this 2X2 matrix represents two edges (equivalent to four connections), node 1 is connected to node 2, and node 5 is connected to node 11. Each of nodes 1, 2, 5, and 11 all has single connection. Since the graph is undirected, [1 2] and [2 1] are identical and mean the same thing, an edge between node 1 and node 2. Scale-free network by the idea of preferential attachment is explained in detail in this link: http://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model

Below is my code to create this network given the total population, N, and average number of connections, d, each node/individual has. I am looking for a more efficient way to implement this because I need to create a network of couple million nodes and using for-loop with this many nodes will take too much time.

Any suggestions would be greatly appreciated! Thanks in advance.

%%%%%%%%%%%%%%CODE%%%%%%%%%%%%%%%%%%%

N=624; %number of nodes
d=4; %average # of connections 

%We first start with fully connected network with d+1 nodes.  Which ensures the average degree of d.  For example, if d=2, then we generate [1 1; 1 2; 1 3; 2 1; 2 2; 2 3; 3 1; 3 2; 3 3], which is [s t], then we get rid of self-loops such as [1 1; 2 2; 3 3] and repetitions such as [1 2; 2 1].  After these implementation, the resulting edge list, [s t] will look like [1 2;1 3;2 3], and each node has d=2 connections as we intended.

s=floor((d+1:d+d*(d+1))/(d+1))';
t=(mod(d+1:d+d*(d+1),d+1)+1)';

LoopIndex=find(s==t);
s(LoopIndex)=[];
t(LoopIndex)=[];

[C a]=unique(sort([s t],2),'rows');
RepeatIndex=setdiff(1:length(s),a);
s(RepeatIndex)=[];
t(RepeatIndex)=[];

%Following two if statements and declaration of the 'degree' variable computes the number of connections each node has.  We need this, because new nodes that will be added in to the small fully connected network we created above proportionally according to their degree (number of connections)

f=accumarray(s,ones(length(s),1));
if length(f)<d+1
    for mm=1:(d+1-length(f))
        f(length(f)+1)=0;
    end
end
g=accumarray(t,ones(length(t),1));
if length(g)<d+1
    for mm=1:(d+1-length(g))
        g(length(g)+1)=0;
    end
end

degree=f+g;

%Following for-loop will add the remaining nodes we need (N-(d+1) nodes) to the small fully connected network one by one.  As a new node joins the network, it will connect to the existing nodes proportional to their degree (which is the idea of preferential attachment, new incoming individuals are likely to connect to the popular nodes, nodes with higher degree, than less popular ones).  New incoming node will make d/2 connections, because each added edge introduces two connections (one for originating node and another for destination node), in order to ensure that the network has the average degree of d that we specified. (when d is odd, average connections of the network will be approximately d, not exactly d, since we choose to connect with floor(d/2) nodes).

for i=d+2:N

prob=degree/sum(degree);
nodeProb=[(1:length(prob))' prob];

Connect=datasample(nodeProb(:,1),floor(d/2),'Replace',false,'Weights',prob)
s=[s;ones(floor(d/2),1)*i];
t=[t;Connect];

f=accumarray(s,ones(length(s),1));
if length(f)<i
    for mm=1:(i-length(f))
        f(length(f)+1)=0;
    end
end
g=accumarray(t,ones(length(t),1));
if length(g)<i
    for mm=1:(i-length(g))
        g(length(g)+1)=0;
    end
end
degree=f+g;
end

edges = [s t];

[Mark Spong]

I ran your code and it works great. I can't comment about whether or not there is a simpler way.