\section{Introduction}\label{sec:introduction}}
\else
\section{Introduction}
\label{sec:introduction}
\fi




\IEEEPARstart{T}{his} demo file is intended to serve as a ``starter file''
for IEEE Computer Society journal papers produced under \LaTeX\ using
IEEEtran.cls version 1.8b and later.
I wish you the best of success.

\hfill mds
 
\hfill August 26, 2015

\subsection{Subsection Heading Here}
Subsection text here.


\subsubsection{Subsubsection Heading Here}
Subsubsection text here.













\section{Conclusion}
The conclusion goes here.








\appendices
\section{Proof of the First Zonklar Equation}
Appendix one text goes here.

\section{}
Appendix two text goes here.


\ifCLASSOPTIONcompsoc
 
  \section*{Acknowledgments}
\else
 
  \section*{Acknowledgment}
\fi


The authors would like to thank...


\ifCLASSOPTIONcaptionsoff
  \newpage
\fi






\section{Introduction}
This demo file is intended to serve as a ``starter file''
for IEEE Computer Society conference papers produced under \LaTeX\ using
IEEEtran.cls version 1.8b and later.
I wish you the best of success.

\hfill mds
 
\hfill August 26, 2015

\subsection{Subsection Heading Here}
Subsection text here.


\subsubsection{Subsubsection Heading Here}
Subsubsection text here.













\section{Conclusion}
The conclusion goes here.







\ifCLASSOPTIONcompsoc
 
  \section*{Acknowledgments}
\else
 
  \section*{Acknowledgment}
\fi


The authors would like to thank...








\section{Introduction}\label{sec:introduction}}



\IEEEPARstart{W}{ith} the rapid developments of digital society, an increasing number of real networks appear in our lives, which are large scale and complex. Real networks can be regarded as complex networks in different fields, such as industri-al networks\cite{jin2010novel}, communication networks\cite{tizghadam2010betweenness}, transporta-tion networks\cite{puzis2013augmented}, social networks\cite{brandes2001faster}, and the Internet of Things\cite{mehmood2016proving}. Research on complex networks has attracted considerable attention from academic fields\cite{2005},\cite{dorogovtsev2010complex}, to deeply study the characteristics of a complex network, network science has emerged\cite{barabasi2013network}.  Network science focuses on network modeling, centrality measures and global characteristics of networks, link prediction and recom-mendation algorithms based on networks, the control and optimization of networks, and network dynamics. \par
 In order to characterize the importance of nodes and edges in the networks, a large volume of research on networks has been devoted to the concept of centrality, such as degree centrality\cite{freeman1978centrality}, closeness centrality\cite{bavelas1950communication}, betweenness centrality\cite{freeman1977set}, Katz centrality\cite{katz1953new} and PageRank centrality\cite{broder2000graph}. Betweenness centrality is a core measure used to estimate the influence of vertices and edges in complex networks, and the calculation of betweenness centrality was incorporated into DARPA HPCS SSCA in 2008, which is a benchmark used exten-sively to evaluate the performance of emerging high-performance computing architectures for graph analy-sis\cite{madduri2009faster},\cite{DARPA}.\par
The scales of complex networks grow rapidly, which makes comparative centrality analysis with more than one million nodes prohibitive. Brandes'  algorithm is the fastest algorithm for computing betweenness centrality\cite{bader2007approximating},\cite{furno2018fast}; however, it cannot fully meet the demand for analysis of complex networks. There is a bottleneck in the existing algorithm, since Brandes' algorithm is based on single source shortest-path traversal algorithms and the accumulation technique\cite{brandes2001faster}. To a certain extent, this bottleneck restricts the acceleration of Brandes' algorithm by parallel computing and thus makes it perform poorly on large-scale networks.\par
In previous studies on the calculation of between-ness centrality through parallel and distributed compu-tation, there are four main technical schemes: solving multiple single-source shortest-path problems at the same time \cite{madduri2009faster},\cite{bader2006parallel},\cite{tan2009parallel}; accelerating Brandes' algorithm with GPUs  \cite{mclaughlin2014scalable},\cite{davidson2014work},\cite{shi2011fast},\cite{jia2012edge}; calculating the shortest paths and then searching for neighbors via parallelized processing after the source node is determined \cite{edmonds2010space}; and dividing and conquering a large network, calculating the subgraphs individually, and synthesizing the calculation results to obtain the betweenness centrality\cite{erdHos2015divide},\cite{wang2013distributed}. These studies have focused on sparse networks and have provided a variety of parallel optimizations of Brandes' algorithm. However, these methods also severely de-generate in performance on dense networks. Exactly computing the betweenness centrality in large-scale net-works faster has become a critical issue in network sci-ence, and there is an urgent need for its application.\par
 In this paper, we propose a more efficient algorithm based on aDjacency mAtrix operations for computing betWeenness ceNtrality, named DAWN. The feasible network scale for the computation of betweenness cen-trality can be extended much further with the rapid de-velopment of hardware. DAWN has the following char-acteristics:\par
(1) DAWN is suitable for sparse networks, dense networks, unweighted networks, weighted networks, undirected networks and directed networks.\par
(2) The time complexity of DAWN depends only on the number of nodes and is insensitive to the density of the networks. Hence, its performance advantage is more prominent in large-scale networks.\par
(3) The operations of DAWN are mainly matrix op-erations which can be accelerated with GPUs, and im-provements to hardware will bring greater computation efficiency.\par
In Section 2, we introduce the theoretical foundation of betweenness centrality and its typical algorithms on complex networks. In Section 3, we describe the novel algorithm based on the adjacency matrix. In Section 4, we demonstrate the experimental results and a comparative performance analysis with existing algorithms. Finally, we conclude the work of this paper and propose future problems in Section 5.
\section{Theoretical Foundation}
In this section, we describe the theoretical foundation of this paper, in three aspects: the concept and properties of betweenness centrality, typical algorithms for computing betweenness centrality and an introduction to GPUs.
\subsection{Betweenness Centrality}
Bavelas et el. proposed betweenness centrality for social network analysis and research\cite{bavelas1950communication}. Freeman created a rigorous mathematical definition of betweenness centrality and proposed a method of evaluating the influence of nodes based on betweenness centrality\cite{freeman1978centrality}. This method is based on the fact that a critical node with a high frequency on the shortest paths of a network has a greater impact in the network than a subordinate node with a low frequency on the shortest paths of the network, and and it has complexities of $\Theta(n^3)$ in time and $\Theta(n^2)$ in space, where n is the number of nodes in the network.\par
Brandes proposed a fast betweenness centrality algorithm that requires $\Theta(m+n)$ space and runs in $\Theta(mn)$ and $\Theta(mn+n^2log_2n)$ time on unweighted and weighted networks, respectively, where $m$ is the number of links \cite{brandes2001faster}. This algorithm has been integrated into a notable open-source project in network analysis "NetworkX" \cite{hagberg2008exploring}.\par
Construct a connected network $G=(V,E)$ according to the physical structure of the network, where $V$ represents the set of nodes in the network and $E$ represents the set of edges. We use $n$ and  $m$ to denote the numbers of nodes and edges in the networks, respectively. Betweenness centrality is usually divided into node betweenness centrality and edge betweenness centrality. The definition of node betweenness is the ratio of the number of shortest paths passing through $v$ between $v_i$ and $v_j$ to the total number of shortest paths between two nodes in the network. The betweenness of node $v$ in the network is defined as:
\begin{equation}
	Bet(v)=\sum_{v_i,v_j \in V} \frac {\sigma(v_i,v_j \mid v)}{\sigma(v_i,v_j)}, \label{1}
\end{equation}
where $\sigma(v_i,v_j )$ represents the number of shortest paths passing between $v_i$ and $v_j$ and $\sigma(v_i,v_j \mid v)$ represents the number of shortest paths passing through $v$ between $v_i$ and $v_j$.\par
The definition of edge betweenness is similar to that of node betweenness, except that the object is changed from nodes to edges. The betweenness of edge e is defined as:\par
\begin{equation}
	Bet(e)=\sum_{v_i,v_j \in V} \frac {\sigma(v_i,v_j \mid e)}{\sigma(v_i,v_j)}. \label{2}
\end{equation}

To facilitate comparative experiments, we adopt the normalized form of betweenness centrality in this paper, which is defined as follows:\par

\begin{subequations}
\begin{align}
		bet(v)&=\frac{2}{(n-1)(n-2)} \times Bet(v),\\
		bet(e)&=\frac{2}{n(n-1)} \times Bet(e).
		\label{3}
\end{align}
\end{subequations}

	Wang et al. proved that the sums of the normalized form of betweenness centrality satisfy the following identities in a connected network\cite{wang2020betweenness},\cite{WWZMC21}:
	
	\begin{subequations}
		\begin{align}
			\sum_{v \in V}bet(v)&=l-1,\\
			\sum_{e \in E}bet(e)&=l.
			\label{4}
		\end{align}
	\end{subequations}
where $l$ is the average length of the shortest paths in the network.\par
\subsection{Typical Algorithms of Betweenness Centrality}
At present, the most known algorithm that exactly com-putes the betweenness centrality is Brandes' algorithm and Matthias et al. proposed its optimization\cite{bentert2018adaptive}. It is computationally too expensive for us to calculate be-tweenness centrality based on Brandes' algorithm in large-scale networks with tens of thousands of nodes. To solve the computation problem, Matteo Riondato et al. proposed an approximate Brandes algorithm (the ran-domized approximate Brandes algorithm, named the RA-Brandes algorithm)\cite{riondato2016fast}.\par
The motivation to complete operations faster makes approximation algorithms a viable choice in scenarios with low precision requirements. Therefore, approximate calculation methods have become an important research field. In this subsection, we introduce Brandes' algorithm, the K-shell approximation algorithm and the RA-Brandes algorithm.\par
\subsubsection{Brandes' Algorithm}
Given network $G=(V,E)$, the betweenness centrality of node $v$ is defined as:
\begin{equation}
	Bet(v)=\sum_{s\not=t\not=v\in V} \frac{\sigma_st(v)}{\sigma_st}=\sum_{s\not=t\not=v\in V}Bet_st(v).\label{5}
\end{equation}
Brandes proposed the following formula to calculate betweenness centrality:
\begin{equation}
	Bet_{s\cdot}(v)=\sum_{w:v\in P_s(w)} \frac{\sigma_{sv}}{\sigma_{sw}}[1+Bet_{s\cdot}(w)],\label{6}
\end{equation}
$Bet_{s\cdot}(v)$represents the betweenness centrality of node $v$, which has the shortest path to other nodes, with node $s$ as the source node in the network. As the network becomes increasingly dense, the performance of Brandes' algorithm tends to degenerate, and it runs in $\Theta(n^{2k})$ time on both unweighted and weighted networks, where $k$ is the length of the diameter. \par
\subsubsection{K-Shell Approximation Algorithm}
In 1983, Seidman introduced the concept of K-core decomposition\cite{seidman1983network}. In 2010, Maksim Kitsak et al. proposed the K-shell decomposition algorithm, which can be used for fast approximate calculation of betweenness centrality\cite{kitsak2010identification}.\par
The K-shell decomposition algorithm layers the network according to degree, forms subnets with degree $k$, and extracts some nodes to form a new subnet within each subnet. The new subnets have some of the characteristics of the original networks, and betweenness centrality indices can be calculated quickly, however they are approximate. The open-source project NetworkX integrates this algorithm to perform the approximate computation of betweenness centrality.\par
\subsubsection{RA-Brandes Algorithm}
Matteo Riondato et al. proposed the RA-Brandes algorithm\cite{riondato2016fast}, which appended a sampling strategy to Brandes' algorithm by selecting a certain subset of the entire network as a sampling set to calculate the betweenness centrality of node $v$.\par
In the open-source project Neo4j, the RA-Brandes algorithm uses two sampling strategies\cite{Neo4j}. The first sampling strategy randomly selects nodes from the network to form a sample set. The sampling probability can be adjusted, and the default value is $\frac{{log}_{10}{N}}{e^2}$. If the probability is 1, the algorithm runs in the same way as the exact betweenness centrality algorithm. The second sampling strategy is to select nodes by degree. The mean degree of the nodes is calculated, and only the nodes whose degrees are higher than the mean are selected to form the sample set.\par
The RA-Brandes algorithm pays more attention to nodes with higher degrees, and the K-shell algorithm can retain nodes with diverse degrees.\par
\subsection{GPU Acceleration}
A graphics processing unit (GPU) is a visual processor focused on graphics processing. Researchers have found that GPUs are more suitable for large-scale computing than central processing units (CPUs). This led to the birth of the general-purpose graphics processing unit (GPGPU).\par
NVIDIA has developed a dedicated GPU program-ming language, CUDA, and released high-performance computing libraries such as cuBLAS\cite{CUDAtoolkit}. CPU, GPU, FPGA, etc. processors each have fields in which they perform best\cite{OverviewofHeterogeneous}, and GPUs are more suitable for linear algebra.
In experiments, the matrix multiplication function in cuBLAS takes only a ten thousandth of the time of a function written based on Openmp. As the scale of the matrix grows, this gap continues to increase. In view of the powerful performance of GPUs in dealing with matrix operations, we utilize a GPU to accelerate the computa-tion in this paper.
\section{Designs of the Betweenness Centrality Algorithm}
DAWN adopts a GPU for computational acceleration and uses Openmp and CUDA/C++ to compile. This requires $\Theta (n^2)$ space and $\Theta (n^{2.387})$ time (approximate values in the limit state), where $n$ is the number of nodes in the network. DAWN addresses the issue of calculating the betweenness centrality index in large-scale networks. \par
Classic complex network models mainly include the Erdős–Rényi random network model\cite{erdos1960evolution}, Barabási–Albert scale-free network model\cite{barabasi1999emergence},\cite{albert1999diameter},\cite{faloutsos1999power}, and the Watts-Strogatz small world network model\cite{watts1998collective},\cite{amaral2000classes}. We choose the ER and BA models as the reference network models for algorithm testing. In this section, we introduce an algorithm for unweighted networks, an algorithm for weighted networks, the theoretical design of auxiliary functions and the parallel computing model of the algorithms.


\subsection{Unweighted Networks}
The unweighted network $G=(V,E)$ is conveniently described as an adjacency matrix:
\begin{equation}
	A=
	\begin{bmatrix}
		a_{11}&\cdots&a_{1n}\\
		\vdots&\ddots&\vdots\\
		a_{n1}&\cdots&a_{nn}
	\end{bmatrix}.
	\label{7}
\end{equation}
The entry $a_{ij}=1$ in adjacency matrix A means that there is an edge with a distance of 1 between node $i$ and node $j$. We strip node $v$ and the edges connecting it to other nodes from the network, which means setting all the values in row $v_i$ and column $v_j$ to 0 in adjacency matrix A to obtain a matrix B with node $v$ stripped out. We define the iterative formulas of matrices $A^{(k)}$,$ B^{(k)}$ and $C^{(k)}$ as:
\begin{equation}
		\begin{cases}
		A^{(k+1)}=A^{(k)}\times A^{(1)},\\
		B^{(k+1)}=B^{(k)}\times B^{(1)},  \quad\text {$1\leq k < k+1 \leq n-1$} \\
		C^{(k+1)}=A^{(k+1)} - B^{(k+1)}.
	\end{cases}
\label{8}
\end{equation}
Although the definitions of the three matrices are similar, the slight gaps between them are the focus of the DAWN algorithm. We define them as follows:\par
(1) The element $a_{ij}^{(k)}$ in matrix $A^{(k)}=\left(a_{ij}^{(k)}\right)$ represents the number of paths with length $k$ between nodes $i$ and $j$ in the network.\par
(2) The element $b_{ij}^{(k)}$ in matrix $B^{(k)}=\left(b_{ij}^{(k)}\right)$ represents the number of paths with length $k$ that do not pass through node $v$ in going between nodes $i$ and $j$ in the network.\par
(3) The element $c_{ij}^{(k)}$ in matrix $C^{(k)}=\left(c_{ij}^{(k)}\right)$ represents the number of paths with length $k$ that pass through node $v$ in going between nodes $i$ and $j$ in the network.\par
Given the above definitions of the matrices, we propose a new betweenness centrality calculation formula:\par
\begin{equation}
	{Bet}_{AM}(v)=\sum_{1 \leq k \leq n-1}^{1\leq i \cap j\leq n \cap i\not = j} F_{ij}(k). 
	\label{9}
\end{equation}\par
We define the function $F_{ij}\left(k\right)$ as:
\begin{equation}
	F_{ij}(k)=
	\begin{cases}
		\begin{aligned}
				&\frac{c_{ij}^{(k)}}{a_{ij}^{(k)}},&\quad\text{$a_{ij}^{(k)}\not =0 \cap i\not=j \cap \sum_{1 \leq p \leq k-1}a_{ij}^{(p)}=0$}  \\
			&0.&\quad\text{Otherwise}
		\end{aligned}
	
	\end{cases}
	\label{10}
\end{equation}\par
To better analyze betweenness centrality, we suggest using the normalized formulas (11) and (12) for calculation. We define the formula used for undirected networks as: 
\begin{equation}
	bet_{AM}(v)=\frac{2}{(n-1)\times(n-2)}\sum_{1 \leq k \leq n-1}^{1\leq i \cap j\leq n \cap i\not = j}F_{ij}(k).
	\label{11} 
\end{equation}
Similiary, the formula used for directed networks is:
\begin{equation}
	bet_{AM}(v)=\frac{1}{(n-1)\times(n-2)}\sum_{1 \leq k \leq n-1}^{1\leq i \cap j\leq n \cap i\not = j}F_{ij}(k).
	\label{12} 
\end{equation}\par
We use pseudocode to describe the algorithm above. 
\begin{algorithm}
		\caption{ Unweighted Networks} 
		\label{ABCLFRS}
		\textbf{Input}: $ A^{(n)},B^{(n)}$\\
		\textbf{Output}: $bet_{AM}(\nu)$
		\begin{algorithmic}[1]
			\Function {MAIN}{$A^{(n)},B^{(n)},bet_{AM}(\nu)$}
			\State $bet_{AM}(\nu) \gets 0 $
			\For{$i = 1 \to n-2 $}
			\State $A^{(i+1)} \gets A^{(i)} \times A$
			\State $B^{(i+1)} \gets B^{(i)} \times B$
			\State $C^{(i+1)} \gets A^{(i+1)} - B^{(i+1)}$
			\EndFor
			\For{$k = 1 \to n $}
			\State $sum \gets sum + F(k)$
			\EndFor
			\State $bet_{AM}(\nu) \gets \cfrac{2}{(n-1)\times(n-2)} \times sum$
			\State \Return{$ bet_{AM}(\nu)$}
			\EndFunction
		\end{algorithmic} 
\end{algorithm}  
\subsection{Weighted Networks}
The weighted network $G=(V,E,W)$ is conveniently described as an adjacency matrix:
\begin{equation}
	A_w=
	\begin{bmatrix}
		w_{11}&\cdots&w_{1n}\\
		\vdots&\ddots&\vdots\\
		w_{n1}&\cdots&w_{nn}
	\end{bmatrix}.
	\label{13}
\end{equation}
Let w be the weight indices of the edges. We assume that $w\in N$, where $N$ is the set of natural numbers. In matrix $A_w$, $w_{ij}$ represents the weight of the edge between nodes $i$ and $j$. First, we define matrices $A_w^{(k)}$, $B_w^{(k)}$, $C_w^{(k)}$ and $H_{A_w}^{(k)}$:\par
(1) Matrix $A_w^{(k)}$ represents the number of paths with weight k between nodes $i$ and $j$ in the network. Matrix $A_w^{(1)}$ represents the adjacency matrix corresponding to the subnetwork formed by the edges with weight 1 in the adjacency matrix $A_w$. \par
(2) Matrix $B_w^{(k)}$ represents the number of paths with weight $k$ that do not pass through node $v$ between nodes$i$ and $j$ in the network.  Matrix $B_w^{(1)}$  represents the adjacency matrix corresponding to the subnetwork formed by the edges with weight 1 that do not pass through node $v$ in the adjacency matrix $A_w$.\par
(3) Matrix $C_w^{(k)}$  represents the number of paths with weight $k$ that pass through node $v$ between nodes $i$ and $j$ in the network. Matrix $C_w^{(k)}$ represents the adjacency matrix corresponding to the subnetwork formed by the edges with weight 1 that pass through node $v$ in the adjacency matrix $A_w$.\par
(4) Matrix $H_{A_w}^{(k)}$ represents the adjacency matrix corresponding to the subnetwork formed by the edges with weight $k$ in the adjacency matrix $A_w$. If there is no edge with weight $k$ in the network, then $H_{A_w}^{(k)}$ is a null matrix, and we also define $H_{A_w}^{(1)}$ as a null matrix.
Second, we define iterative formulas for matrix $A_w^{(k)}$ as follows:\par
\begin{equation}
	\begin{cases}
		\begin{aligned}
		A_w^{(k+1)}=&A_w^{(k)}\times A_w^{(1)} + 	H_{A_w}^{(k+1)}, &\quad\text{$1\leq k < k+1\leq w_0$} \\
		A_w^{(k)}=&A_w^{(r)}\times A_w^{(\lambda w_0)}.  &\quad\text{$w_0< k < n$}
	\end{aligned}
	\end{cases}
	\label{14}
\end{equation}\par
Third, similar to matrix $A_w^{(k)}$, we define iterative formulas for matrix $B_w^{(k)}$ as follows:
\begin{equation}
		\begin{cases}
		\begin{aligned}
			B_w^{(k+1)}=&B_w^{(k)}\times B_w^{(1)} + 	H_{B_w}^{(k+1)}, &\quad\text{$1\leq k < k+1\leq w_v$} \\
			B_w^{(k)}=&B_w^{(r)}\times B_w^{(\lambda w_v)}.  &\quad\text{$w_v< k < n$}
		\end{aligned}
	\end{cases}
		\label{15}	
\end{equation}\par
$H_{Bw}^{(k)}$ represents the adjacency matrix corresponding to the subnetwork formed by the edges with weight $k$ in the adjacency matrix $B_w^{(k)}$, and $H_{Bw}^{(1)}$ is a null matrix. $w_0$ represents the upper limit of the weight of the weighted network represented by matrix $A_w^{(k)}$, and $w_v$ represents the upper limit of the weight of the weighted network represented by matrix $B_w^{(k)}$.
Finally, according to \eqref{8}, \eqref{14} and \eqref{15}, we define the matrix $C_w^{(k)}$ as:
\begin{equation}
	C_w^{(k)}=A_w^{(k)}-B_w^{(k)}.
	\label{16}	
\end{equation}\par
We define the variable parameters $r$ and $\lambda$ in the above formulas as:
\begin{equation}
	\begin{cases}
		r=k \text{ mod } w_0, \lambda=\lfloor\frac{k-r}{w_0}\rfloor, &\quad\text{In formulas \eqref{14}}\\
		r=k \text{ mod } w_v,
		\lambda=\lfloor\frac{k-r}{w_v}\rfloor. &\quad\text{In formulas \eqref{15}}\\
	\end{cases}
	\label{17}	
\end{equation}\par
The calculation processes of the unweighted and weighted networks are the same, except for the different matrix definitions and iterative formulas. Finally, we calculate the betweenness centrality of unweighted and weighted networks with \eqref{9}$\sim$\eqref{12}.\par
\subsection{Auxiliary Functions}
For a network with $n$ nodes, DAWN needs at most $k_0$ loops to complete the calculation, where $k_0$  is the maximum length of the shortest paths. We build auxiliary functions to reduce the number of loops, which are the search function based on the matrix and the threshold function.
\subsubsection{Search Function for the Matrix}
We define the matrix $\Gamma^{(n)}$ to register the value of the function $F(k)$. According to the definition of $F(k)$, the new operations do not affect the betweenness centrality indices after the matrix is filled, so we suggest ending the loop at this time and calculating ${bet}_{AM}(v)$.\par
For example, we define the matrix $\Gamma^{(8)}$ as:
\begin{equation}
		\Gamma^{(8)}=
	\begin{bmatrix}
		\#&\#&\#&\#&\#&\#&\#&\#\\
		\#&\#&\ast&\ast&\ast&\ast&\ast&\ast\\
	    \#&\ast&\#&\ast&\ast&\ast&\ast&\ast\\
		\#&\ast&\ast&\#&\ast&\ast&\ast&\ast\\
		\#&\ast&\ast&\ast&\#&\ast&\ast&\ast\\
		\#&\ast&\ast&\ast&\ast&\#&\ast&\ast\\
		\#&\ast&\ast&\ast&\ast&\ast&\#&\ast\\
		\#&\ast&\ast&\ast&\ast&\ast&\ast&\#\\
	\end{bmatrix}.
	\label{18}
\end{equation}
"$\#$" represents that the value of this position is not included in the calculation of the betweenness centrality indices. "$\ast$" represents that the value is to be calculated, and if the value does not appear, it will finally be recorded as 0. $\Gamma^{(8)}$ above corresponds to the operations of betweenness centrality for node 0 in the networks. The values in row $i_0$ and column $i_0$ are not included in the calculations. We do not consider the situation in which the path starts and ends overlap in calculating the shortest path, so we exclude the case of $i=j$. \par
\subsubsection{Threshold Function}
The threshold function is used to estimate the number of loops based on the average path length of the network. We define the threshold of the Erdős–Rényi network as:
\begin{equation}
	K_0=\lceil \mu \times \frac{ln(n)}{ln\langle deg\rangle}\rceil.
		\label{19}
\end{equation}
$\left\lceil\bullet\right\rceil$ is the ceiling function, $\mu$ is defined as the tolerance coefficient $(\mu\geq1)$, deg represents the average degree of the network, and $\frac{ln(n)}{ln\langle deg\rangle}$ represents the average shortest path length in the network. $\langle deg\rangle$ can be expressed as:
\begin{equation}
	\langle deg\rangle=p \times (n-1),
	\label{20}
\end{equation}
where $p$ represents the connection probability of any pair of nodes in the network.\par
The degree distribution of the BA network obeys $p\left(k\right)\sim k^{-\gamma}$, where $\gamma$ is the power exponent. We define the threshold of the Barabasi-Albert network as:
\begin{equation}
	\begin{cases}
		K_0=\lceil \mu \times ln(ln(n)) \rceil, &\quad\text{$2<\gamma<3$}\\
		K_0=\lceil \mu \times \frac{ln(n)}{ln(ln(n))} \rceil. &\quad\text{$2<\gamma<3$}\\
	\end{cases}
	\label{21}
\end{equation}\par
After adding the auxiliary functions, the maximum number of loops is reduced from $n-1$ to $K_0$, and the time complexity of the algorithm is also reduced. When the number of loops far exceeds $K_0$, special circumstances such as isolated nodes or disconnected subnets in the networks need to be considered. In these circumstances, there is no connected path between some nodes. Even if the operation is performed $n-1$ times, the matrix $\Gamma^{(n)}$ will not be filled, and the calculation can be ended earlier.
\par
\subsection{Parallel Computing Model of the Algorithm}
We propose a parallel computing model named EVENT parallelism for betweenness centrality calculation, which aggregates the advantages of data parallelism and task parallelism. Event parallelism means dividing and conquering the target event and generating subevents, assigning subevents to each thread, selecting the running processor according to the characteristics of each program in the subevents, supporting data sharing and for load balancing across threads.\par
The betweenness centrality calculations of the nodes in the networks are independent of each other. We define the betweenness centrality calculation of the entire network as the target event and the betweenness centrality calculation of each node as a subevent. The calculation time is inversely proportional to the number of threads enabled, and thread conflicts are effectively avoided.\par
\subsection{Optimized Matrix Algorithm}
We achieve acceleration of the matrix algorithms through data precomputation and matrix block multiplication, thereby significantly reducing the time complexity of the algorithm.\par
\subsubsection{Precomputation}
The matrix $A^{(k)}$ is constant for nodes in the network, and the betweenness centrality calculation of each node requires the data of matrices $A^{(1)}\sim A^{(k)}$. We calculate this part of the data before dividing and conquering the event and then set it as shared data. These data do not need to be kept in memory; instead, they can be read from the disk when needed. \par
This reduces the number of matrix multiplication operations from $2k_0n$ to $k_0(n+1)$, when calculating the betweenness centrality of all nodes; here, $k_0$ represents the number of times the algorithm runs until the matrix $\Gamma^{(k)}$ is filled. This optimization does not fundamentally change the spatial complexity of the algorithm, but it significantly reduces its temporal complexity. For large-scale networks, we recommend storing this part of the data on disk to save the relatively limited available memory.\par
\subsubsection{Matrix Block Multiplication}
We adopt the Coppersmith-Winograd matrix multiplication algorithm to propose a novel implementation by GPUs\cite{strassen1969gaussian},\cite{coppersmith1987matrix},\cite{williams2012multiplying}. This novel implementation uses 7 multiplications and 3 nested loops ($\Theta(n^2)$) which are calculated by GPUs and CPUs, respectively, and it can be executed iteratively. The following is the pseudocode for the proposed novel implementation.\par
\begin{algorithm}
	\caption{ Matrix Block Multiplication} 
	\textbf{Input}: $ A^{(n)},B^{(n)}$\\
	\textbf{Output}: $C^{(n)}$
	\begin{algorithmic}[1]
		\Function {MAIN}{$A^{(n)},B^{(n)},C^{(n)}$}
		\State parallel section 1:
		\State $S1\gets A21+A22$, $T1\gets B12-B11 $
		\State $S3\gets A11-A21$, $T3\gets B22-B12$ 
		\State $S2\gets S1-A11 $, $T2\gets B22-T1 $
		\State parallel section 2:
		\State $M5\gets S1 \times T1$
		\State $M7\gets S3 \times T3$
		\State $M6\gets S2 \times T2$
		\State $M1\gets A11 \times B11$
		\State $M2\gets A12 \times B21$
		\State parallel section 3:
		\State $C11\gets M1+M2$, $U2\gets M1+M6 $
		\State $U3\gets U2+M7$, $C22\gets U3+M5$ 
		\State $U4\gets U2+M5 $, $S4\gets A12-S2 $
		\State $T4\gets T2-B21 $
		\State parallel section 4:
		\State $M4\gets S4 \times B22$
		\State $M3\gets A22 \times T4$
		\State parallel section 5:
		\State $C21\gets U3-M4$, $C12\gets U4+M3 $
		\State \Return{$ C^{(n)}$}
		\EndFunction
	\end{algorithmic} 
\end{algorithm}
We believe that this order of calculations achieves the best balance between temporal and spatial complexity. It requires $\Theta(3n^2+7\left(\frac{n}{2}\right)^2)$ space and $\Theta(n^\omega)$ time, where $ 2.387<\omega<2.81$. Changing the order of the calculations may lead to calculation errors. The details of the implementation code are given in the appendix.\par
When the order of the matrix is odd, we add a new row and column to the matrix in which all values are set to 0. This operation, which requires at most $\log_2{n}$ time, does not affect the result of matrix block multiplication and extends the applicability of the algorithm. The temporal complexity of the algorithm slightly increases to $\Theta({(n+\log_2{n)}}^\omega)$.\par
\section{Experimental Results}
We analyze the experimental results in this section, which is divided into two subsections: the introduction of the experimental environment and parameters, the experimental results and analysis. The calculation of betweenness centrality for a network with millions of nodes and connection probabilities at the level of tens of percent can be achieved on a single-GPU system. The experimental results show that a multi-GPU system achieves nearly linear speedup.\par
\subsection{Experimental Environment Introduction}
We compiled the algorithmic solution and compared it with the code based on NetworkX. NetworkX is a Python package for complex networks analysis; it has excellent performance and a wide range of application scenarios\cite{hagberg2008exploring}. We noticed that the maximum of matrices order in multiplications that TASLE V100 can support is 45000 when DAWN is running, and this upper limit is slightly floating. When the scale of the network is above this upper limit, matrix block multiplication becomes an inevitable choice. Table 1 shows some of the parameters of the test machine.\par
\begin{table}[h]
	\centering
	\caption{Partial Parameter Table of the Test Machine}
	\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|c|l|}
	\hline 
	Hardware &	 Parameters\\
	\hline  
	CPU	&  Intel Xeon Gold 6151(2)\\
	\hline  
	RAM	& 1TB\\
	\hline  
	GPU	& NVIDIA TESLA V100(4)\\
	\hline  
	IDE	& Visual Studio 2019\\
	\hline  
	OS	& Windows Server 2019\\
	\hline  
	Toolkit	& CUDA 11.4\\
	\hline  
\end{tabular}
\end{table}
\subsection{Results and Analysis}
Unweighted undirected networks (UWD) and weighted directed networks (WD) are the models with the highest and lowest operating efficiency for Brandes' algorithm, respectively. We show the running time of NetworkX and DAWN on Erdős–Rényi networks and Barabasi-Albert networks with 10000 nodes. The unit for the running time is seconds, with 6 significant digits reserved. The time shown in the table is the average time of 5 tests.
\begin{table}[h]
	\centering
	\caption{Calculation Time on Multiple Network Models}
		\renewcommand{\arraystretch}{1.3}
	\begin{tabular}{|c|c|c|c|}
		\hline 
		MODEL & Connection Probability & Brandes' &	DAWN\\
		\hline  
		\multirow{4}{*} {UWD-ER} &	0.05&	7196.63 & 15467.1\\
				\cline{2-4}  
		&0.1 & 16934.5 & 15694.3\\
				\cline{2-4}  
		&0.2 &	38541.3	& 15413.7\\
				\cline{2-4}  
		&0.5 &	97671.5	& 15573.4\\
				\hline  
		\multirow{4}{*} {WD-ER} &	0.05 &	41164.8 & 15854.2\\
				\cline{2-4}  
	& 0.1 &	84756.4 & 15891.8\\
				\cline{2-4}  
	& 0.2 &	181734 & 15803.0\\
				\cline{2-4}  
	& 0.5 &	458414 & 15872.5\\
					\hline  
		\multirow{4}{*} {UWD-BA} & 0.05 &	17614.5 & 20874.3\\
				\cline{2-4}  
		& 0.1 & 41169.7 & 17502.8\\
				\cline{2-4}  
		& 0.2 &	83451.6 & 17279.6\\
				\cline{2-4}  
		& 0.5 &	217589 & 15327.4\\
							\hline 
	\end{tabular}
\end{table}
The scale, density and complexity of the networks will have a significant impact on the efficiency of Brandes' algorithm. As the connection probability increases, the calculation time of DAWN on the networks decreases. Katzav et al. proposed the phenomenon of network diameter shrinkage\cite{katzav2015analytical}. As the density increases, the diameter of the network shrinks and the calculation time decreases.\par
A. McLaughlin et al. studied how to accelerate Brands' algorithm with GPUs and used the edge traversal speed to evaluate the algorithm performance\cite{mclaughlin2014scalable}. Most of the data used in the research of A. McLaughlin et al. were represented by relatively sparse UWD networks. The public dataset released by Stanford University includes dense networks, such as: Social Circles(Google+)\cite{leskovec2012learning}; Email-Eu-core Temporal Network\cite{paranjape2017motifs}; Dynamic Face-to-Face Interaction Networks\cite{kumar2021deception}; Social Network(MOOC User Action Dataset)\cite{kumar2021deception}; Disease-gene Association Network\cite{davis2019comparative}; Enhanced Tissue-specific Gene Interaction Networks\cite{wang2018network}.\par
\begin{figure}[ht]
	\subfloat[Single-GPU System]{\includegraphics[width=3.5in]{1.pdf}
		\label{a}}
	\\
	\subfloat[Multi-GPU System]{\includegraphics[width=3.5in]{2.pdf}
		\label{b}}
	\caption{Calculation Time of betweenness centrality on GPUs}
	\label{fig_sim}
\end{figure}
These dense networks have appeared in datasets in recent years, which, in some sense, indicates that dense networks have gradually become a research focus in the field of complex network analysis. Therefore, we believe that it is also necessary to evaluate the performance of our algorithm on dense networks.
The results include experiments running on the UWD-ER model with a connection probability of 0.1. We show that the betweenness centrality calculation time DAWN operations can be completed on a single-GPU system, and multi-GPU acceleration has no significant advantage when the number of nodes in the network is fewer than $2^{15}$. We present separate calculation times for a network with a number of nodes in the range of $2^{10}$ to $2^{15}$.\par

Figure 1(b) shows that DAWN achieves approximately linear speedup on the multi-GPU systems. A. McLaughlin et al. gave the calculation time of sampling and edge-parallel methods running on delaunay-n20 via 192 GPUs \cite{mclaughlin2014scalable},\cite{bader2013graph}, which is a random triangulation network with over a million nodes and an average degree of about 3. DAWN has performance advantages within the scope of the experiment and running on an Erdős–Rényi network with connection probability of $10\%$.\par

The distribution of the nodes and edges in a real network is uneven. We performed validation experiments using the real datasets mentioned in this section, and the results show that the computation times for some nodes exceed the values shown in Figure 1, while the times for other nodes are lower than these values. If the node to be calculated is a marginal node or an isolated node in the network,the calculation time will drop significantly. The maximum value of the calculation time for a node is related to the network diameter, but this rarely exists in the network.\par
In the experiments, DAWN shows three performance advantages:\par
(1) On networks with the same number of nodes, as the density increases, DAWN gradually shows greater advantages in computational efficiency;\par
(2) On networks with the same number of nodes, DAWN has a more obvious advantage over Brandes' algorithm in terms of computational efficiency in the directed or weighted network.\par
(3) For networks with an increasing number of nodes, DAWN can continue to process their adjacency matrices even on a single-GPU system and achieves approximately linear speedup on multi-GPU systems.\par

\par
















\section{Conclusion}
We propose a novel algorithm for the betweenness centrality that uses a parallel computing method based on the adjacency matrix to improve the calculation efficiency. It addresses the insufficient calculation efficiency of Brandes' algorithm on dense networks with thousands of nodes and solves the urgent problem in the field of complex network analysis: how to exactly compute the betweenness centrality index on large networks faster. DAWN can be widely used in large-scale network application fields such as social network analysis, information analysis and recommendation on the internet, biological network structure analysis, and communication network vulnerability analysis.








\appendices
\section{}
 Please refer to the GITHUB library for the verification code of this paper.\\
 https://github.com/lxrzlyr/Betweenness-centrality.git
\ifCLASSOPTIONcompsoc
 
  \section*{Acknowledgments}
\else
 
  \section*{Acknowledgment}
\fi


\ifCLASSOPTIONcaptionsoff
  \newpage
\fi










