Implementing the Barabási–Albert model I am writing a code for Barabási–Albert(BA) model with specific node and edges.
The algorithm is almost like [1] as follows:

  1.Add m&lt;N nodes to G.
  2.Connect nodes in G randomly until you get a connected graph.
  3.Create a new node i.
  4.Pick a node j uniformly at random from the graph G. Set p=k(j)/k_tot.
  5.Pick a real number r uniformly at random between 0 and 1.
  6.If r&lt;p then add G(m+step_i,j)=G(j,m+step_i)=1
  7. repeat steps 4-6, k times until each new node has k edges.
  8. repeat steps 3-7 until we have N nodes.


*:The sum of elements of each row or column is the degree(number of edges connected to that node) of each node.

** : k_tot is the total degree of all nodes until step i.
and here is the my implementation in Fortran:

program BA
 implicit none
 integer, parameter :: N = 20
 integer, dimension(n,n) :: A
 double precision :: r,p,r2(2),pj
 integer :: m0,i,j,L0,N1,t,jnod,k,di
 integer :: edge_number,ii,ie,je, degree

 A = 0
 m0 = 5 ! initial nodes
 L0 = 4 ! initial edges
 ii = 0
 ! adding m0 nodes with random connections
 do while (ii &lt; L0)
   do ie = 1, m0
     call random_number(r)
     je = r * (m0 - 1) + 1
     if (A(ie,je) == 0 ) then
       if(ii&gt;=L0) exit
       A(ie,je) = 1
       A(je,ie) = 1
       ii = ii + 1       
     endif
   enddo
 enddo

 edge_number = L0
 N1 = m0
 degree = m0
 ! adding remaining nodes with pereferential attachment
 do t = 1,N-m0
   di = 0
   N1 = m0 + t 
   do while(di &lt; degree)
     call random_number(r)     
     jnod = r * (N1-2)+1
     pj = sum(A(jnod,1:N1-1))/dble(sum(A(1:N1-1,1:N1-1)))       
     call random_number(r)
     if ( r&lt;pj .and. A(N1,jnod)==0 ) then        
        A(N1,jnod) = 1
        A(jnod,N1) = 1
        edge_number = edge_number + 1
        di = sum(A(N1,1:N1))        
     endif
   enddo
 enddo

 do i = 1,N
    write(12,"(*(I5))") (A(i,j),j=1,N)
 enddo
end program


I think the is not efficient at all and as more nodes added, it takes more to get out of loops.Do you have a better idea to make it more efficient?
I'll put some comments inline to the code, and in the end add a complete implementation of the program, with various changes to improve the execution.

program BA
 implicit none
 integer, parameter :: N = 20


For real numbers it is advisable to introduce a real kind with selected_real_kind:

 integer, parameter :: rk = selected_real_kind(15)
 integer :: A(n,n)


Use the real kind introduced above, to declare the reals:

 real(kind=rk) :: r, p, r2(2), pj
 real(kind=rk) :: sumA_q
 integer :: sumA
 integer :: m0, i, j, L0, N1, t, jnod, k, di


Often it is better to use more than single characters for variables (for example something like iNode or iEdge, this makes it visible what the iterator is actually used for.

 integer :: edge_number, ii, ie, je, degree

 A = 0
 m0 = 5 ! initial nodes
 L0 = 4 ! initial edges
 ii = 0
 ! adding m0 nodes with random connections
 rand_cnx: do ! while is actually deprecated, you anyway use 'if exit' in the
              ! loop, so I'd stick to that one. Using a block label to allow
              ! the identification of the loop in the exit statement.
   do ie = 1, m0
     call random_number(r)


You could create a bunch of m0 random numbers all at once here, and use them then subsequently. That might be a little faster, but probably has no really big impact.

     je = r * (m0 - 1) + 1
     if (A(ie,je) == 0 ) then
       A(ie,je) = 1
       A(je,ie) = 1
       ii = ii + 1
       if (ii&gt;=L0) EXIT rand_cnx


No need to compute another random number here, if the limit was reached, leave the loop. Use a block label to directly leave the outer loop.

     end if
   end do
 end do rand_cnx

 edge_number = L0
 N1 = m0
 degree = m0

 ! Precompute sum of submatrix up to m0-1
 sumA = sum(A(1:m0-1, 1:m0-1))

 ! adding remaining nodes with pereferential attachment
 do t = 1,N-m0
   di = 0
   N1 = m0 + t


The summation over the subarray does not depend on the random numbers (only entries outside the subarray at N1 are modified, if I am not mistaken) and can be pre-computed. As N1 is linearly increasing, we can update the sum of the submatrix up to N1, by just adding the newly filled column and row from the previous iteration.

   ! Update sum of submatrix with the previously computed column and row.
   sumA = sumA + sum(A(1:N1-1,N1-1)) + sum(A(N1-1,1:N1-2))
   sumA_q = 1.0_rk / real(sumA, kind=rk)

   do
     call random_number(r)     
     jnod = r * (N1-2)+1
     pj = sum(A(jnod,1:N1-1))*sumA_q
     call random_number(r)
     if ( r&lt;pj .and. A(N1,jnod)==0 ) then        
        A(N1,jnod) = 1
        A(jnod,N1) = 1
        edge_number = edge_number + 1
        di = sum(A(N1,1:N1))
        if (di &gt;= degree) EXIT
     end if
   end do
 end do

 do i = 1,N
    write(12,"(*(I5))") (A(i,j),j=1,N)
 end do
end program BA


With a little more memory effort, you could also precompute all sum(A(jnod,1:N1-1)), use them for the summation of the complete submatrix and in the computation of pj, instead of recalculating the sum everytime the random number picks a specific jnod. This would result in turning the di into an array of length N. Implementing this led me to find that the first additional node has to be connected to all previous nodes, as you set degree=m0. In this case no randomness seems to be involved? I catch this special case in the implementation below:

program BA
  implicit none
  integer, parameter :: N = 20

  integer, parameter :: rk = selected_real_kind(15)
  integer :: A(n,n)

  real(kind=rk) :: r, pj
  real(kind=rk) :: sumA_q
  integer :: sumA
  integer :: m0, i, j, L0, N1, t, jnod
  integer :: di(N)

  integer :: edge_number, ii, ie, je, degree

  A = 0
  di = 0

  m0 = 5 ! initial nodes
  L0 = 4 ! initial edges

  ! adding m0 nodes with random connections
  ii = 0
  rand_cnx: do

    do ie = 1, m0
      call random_number(r)
      je = nint(r * (m0 - 1)) + 1
      if (A(ie,je) == 0 ) then
        A(ie,je) = 1
        A(je,ie) = 1
        di(ie) = di(ie) + 1
        di(je) = di(je) + 1
        ii = ii + 1
        if (ii&gt;=L0) EXIT rand_cnx
      end if
    end do

  end do rand_cnx

  edge_number = L0
  N1 = m0
  degree = m0

  ! Precompute sum of submatrix up to m0-1
  sumA = sum(di(:m0-1)) - sum(A(:m0-1,m0))

  ! adding remaining nodes with pereferential attachment
  do t = 1,N-m0
    N1 = m0 + t

    ! Update sum of submatrix with the previously computed column and row.
    sumA = sumA + sum(A(:N1-2,N1-1)) + di(N1-1)
    sumA_q = 1.0_rk / real(sumA, kind=rk)

    need_all: if (N1-1 == degree) then

      ! Need all nodes to fulfill degree, no randomness required?
      A(N1,:N1-1) = 1
      A(:N1-1,N1) = 1
      edge_number = edge_number + N1-1
      di(:N1-1) = di(:N1-1) + 1
      di(N1) = di(N1) + N1-1

    else need_all

      do
        call random_number(r)
        jnod = nint(r * (N1-2))+1
        ! Only consider not yet connected nodes.
        if (A(jnod,N1) == 0) then
          pj = di(jnod)*sumA_q
          call random_number(r)
          if ( r&lt;pj ) then
             A(N1,jnod) = 1
             A(jnod,N1) = 1
             edge_number = edge_number + 1
             di(N1) = di(N1)+1
             di(jnod) = di(jnod) + 1
             if (di(N1) &gt;= degree) EXIT
          end if
        end if
      end do

    end if need_all

  end do

  do i = 1,N
     write(12,"(*(I5))") (A(i,j),j=1,N)
  end do

end program BA


No guarantees, that this is in any way correct. I guess, this could be further increased by just picking the random nodes out of the not yet connected ones. This would imply that you need to keep an index list of all 0 values for the current node, and then use jnod to pick out of those. This would be more involved book-keeping but would reduce the number of guesses, you need to do.