# Sparse Cholesky Factorization on FPGA Using Parameterized Model

Authors: Yichun Sun; Hengzhu Liu; Tong Zhou
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Cholesky factorization is a fundamental problem in most engineering and science computation applications. When dealing with a large sparse matrix, numerical decomposition consumes the most time. We present a vector architecture to parallelize numerical decomposition of Cholesky factorization. We construct an integrated analytical parameterized performance model to accurately predict the execution times of typical matrices under varying parameters. Our proposed approach is general for accelerator and limited by neither field-programmable gate arrays (FPGAs) nor application-specific integrated circuit. We implement a simplified module in FPGAs to prove the accuracy of the model. The experiments show that, for most cases, the performance differences between the predicted and measured execution are less than 10%. Based on the performance model, we optimize parameters and obtain a balance of resources and performance after analyzing the performance of varied parameter settings. Comparing with the state-of-the-art implementation in CPU and GPU, we find that the performance of the optimal parameters is 2x that of CPU. Our model offers several advantages, particularly in power consumption. It provides guidance for the design of future acceleration components.


---

## Body

## 1. Introduction

In engineering and science computations, the solution of large sparse linear system of equations is a fundamental problem in most applications. In general, two methods for solving this problem exist: the direct and iterative methods. In the direct method, when focusing on factorization algorithms, we have LU decomposition, QR decomposition, and Cholesky decomposition [1]. It should be noted that Cholesky decomposition is a special form of LU decomposition that deals with symmetric positive definite matrices. Its computational complexity is roughly half of LU algorithm. The Cholesky algorithm is widely used in applications because most matrices involved are symmetric positive definite.The solution of large sparse linear system of equations can be adapted quite easily to parallel architectures. Supercomputers and multiprocessor systems became the main focus of this kind of research [2–4]. However these systems cost highly in making marketable products [5].There are a lot of software- and hardware-based approaches developed to get better solution in Cholesky decomposition. Field-programmable gate arrays (FPGAs) have unique advantages in solving such problems. First, according to the characteristics of the algorithm, a dedicated architecture could be designed [6]. Second, the computational resources, memory, and I/O bandwidth are reconfigurable. Third, the energy consumption and cost are relatively low under the premise of meeting certain performances. Fourth, it provides experimental and verification platform for future coprocessor design. Yang et al., 2010 [7, 8], presented a Cholesky decomposition based on GPUs and FPGAs. The results show that the dedicated FPGA implementation has the highest efficiency. The design in FPGAs has been able to combine parallel floating-point operations with high memory bandwidth in order to achieve higher performance than microprocessors [9].A sparse matrix, in which the ratio of nonzero values is low, is especially suitable for direct methods. In our implementation, we chose the Cholesky factorization method over LU decomposition.The reason why we chose the Cholesky factorization is less operation. The matrix is symmetric and half of factors in matrix need to be calculated and suitable for parallel architectures. In parallel, there is less communication within each interprocessor. There are a lot of work in parallelizing the Cholesky factorization [3, 4]. The memory requirement is significantly lower than other direct methods, because only the lower triangular matrix is used for factorization and the intermediate and final results (the Cholesky factor L) are overwritten in the original matrix.The three algorithms in the implementation of Cholesky factorization are row-oriented, column-oriented, and submatrix algorithms. The difference among these three algorithms is in the order of the three loops [3]. Among these three algorithms, the column-oriented algorithm is the most popular algorithm for sparse matrix. Once each column is regarded as a node, data dependency of each node in Cholesky factorization can be described as an elimination tree [10]. In the tree, the parent nodes depend on their descendants while their descendants must update them before being decomposed. Thus, the process of factorization can be taken as a process of eliminating the tree from bottom to top. Two main operations of the algorithm exist: scaling operation in columns (cdiv()) and updating operation between dependent columns (cmod()). Basing on the updating order, the researchers found that the algorithm can be classified into left-looking and right-looking algorithms.A number of studies in sparse Cholesky factorization optimized in GPU are found. Vuduc and Chandramowlishwaran mapped the submatrix tasks directly to the GPU [11], George scheduled the tasks on the GPU [12], and Lucas et al. used multiple CPU threads to take small-scale computing tasks in parallel and used GPU to speed-up large-scale computing tasks [13]. The small-scale computing tasks are in the bottom of the elimination tree whereas the large-scale computing tasks are in the top of the elimination tree. Lacoste et al. modified the task scheduling module of Cholesky factorization algorithm to provide a unified programming interface on CPU-GPU heterogeneous platforms [14]. Hogg et al. proposed a supernodal algorithm of HLS_MA87 solver for multicore CPU platforms [15]. Chen et al. introduced a GPU version of the CHOLMOD solver for CPU-GPU platforms [16]. Zou et al. improved the generation and scheduling of GPU tasks by queue-based approach and designed a subtree-based parallel method for multi-GPU system [17].To decrease the storage and computation complexity, the sparse matrix would be compressed into a dense format (CSC format, Compressed Sparse Column Format). However, in the process of Cholesky factorization, many original zero entries will turn to be nonzeros. They increase the overhead memory and dynamically change the compressed data structure. To avoid this situation, the matrix should be preprocessed and decomposed in symbol before undergoing numerical factorization (when preprocessed and decomposed in symbol, we could generate nonzero pattern of matrix and its elimination tree. When we reorganize the data, we would get the vector array). Preprocessing exposes parallelism by reordering matrix columns. Symbol decomposition is then executed in advance to identify the sparsity pattern of the result matrix. With the symbol factorization technology being well-developed, numerical decomposition is the most time consuming. This paper discusses and evaluates the performance of numerical decomposition based on the preprocessing and symbol factorization results.After analyzing Cholesky algorithm when dealing with a large sparse matrix, we rewrite the algorithm in a vector version. In order to explore the relationship between parallel scale and performance, to predict the performance on field-programmable gate arrays (FPGAs) platform, we also designed corresponding vector architecture and construct a parameterized model. The performance depends on the underlying architecture and system parameters. The capacity of on-chip memory, the amount of computational resources, and the size of I/O bandwidth are critical parameters. Our model takes these factors into consideration and focused on their relationship. By using the model, we obtain the best performance when choosing the capacity of on-chip memory and the amount of computational resources.The rest of this paper is organized as follows. Section2 introduces the vector version of Cholesky algorithm. Section 3 describes our performance model in ideal, memory-bound, and computational-bound cases. Section 4 proves the consistency of the performance of our model and the implementation; then it analyzes the parameters and tries to find an optimal performance.


## 2. Cholesky Factorization

Cholesky factorization decomposes a symmetric positive definite matrix intoL×LT where L is a lower triangular matrix. The main operations of the algorithm are scale operation in columns (cdiv(k)) and update operation between dependent columns (cmod(k,i)). cdiv(k) divides the subdiagonal elements of column k by the square root of the diagonal element in that column, and cmod(k,i) modifies column k by subtracting a multiple of column i. The algorithm can be classified into left-looking and right-looking algorithm based on updating order. Right-looking algorithm performs all the updates cmod(k,i) using column k immediately after the cdiv(k). By contrast, left-looking algorithm accumulates all necessary updates cmod(k,i) for column i until before the cdiv(i).From Algorithm1 we can note that the parallelism is mainly among the nondata dependency columns and the elements in a column, corresponding to the parallel elimination among different path in the elimination tree(internode parallelism) and the parallel scaling among different subdiagonal elements in a node(intranode parallelism). Besides, because of the columns in dense matrix depending on each other strictly in order, the parallelism among columns only exists in sparse matrix decomposition.Algorithm 1:Right-looking Cholesky algorithm.
Input: A
Output: L
(1)L=A
(2)for  i=1; i≤n; i++  do
(3)l(i,i)=sqrt(l(i,i));               /∗cdiv(k)∗/
(4)for  j=1; j≤n; j++  do
(5)l(j,i)=l(j,i)/l(i,i);
(6)end
(7)for  k=i+1; k≤n; k++  do
(8)if  l(k,i)≠0  then
(9)for  p=k; p≤N; p++  do
(10)l(p,k)=l(p,k)-l(k,i)×l(p,i);    /∗cmod(k,i)∗/
(11)end
(12)end
(13)end
(14)endIn Algorithm1, lines (2), (7), and (9) are a progressive cycle of three layers. The amount of memory access per layer to the array l is O(n) (n is the scale of the matrix). The complexity of memory access is O(n3) whether it is left-looking or right-looking. For the scale of matrix increasing, the complexity is unacceptable when a platform lacks cache space and the on-chip memory is small [18]. Duff and Reid [19] proposed a multifrontal Cholesky algorithm. The method reorganizes the overall factorization of a sparse matrix into a sequence of partial factorizations of dense smaller matrices. The multifrontal algorithm accumulates the update factors into a matrix. In this way, the algorithm reduces the memory access time [20] and the update operation.Algorithm2 shows its content. The decomposition of every column only needs to read its front matrix F and its corresponding update matrix U. The descendants updating effects are accumulating in the matrix U, thus, avoiding the access of its descendants. In Algorithm 2, ⊕ means the extended addition operation between two sparse matrices.Algorithm 2:Multifrontal Cholesky algorithm.
Input: A
Output: L
(1)L=A
(2)for  i=1; i≤N; i++  do
(3)i,i1,…,ir are the abscissas of non-zeros L(:,i);
(4)c1,…,cs are the child nodes of node i;
(5)U=Uc1⊕⋯⊕Ucs; /∗  construct the update matrix of node  i∗/
(6)Fi=ai,iai1,j⋯air,jai1,j⋮air,j⋯0⊕U; /∗   update front matrix: cmod∗/
(7)Fi→li,i0li1,i⋮Ilir,i100Uili,ili1,i⋯lir,i0I   /∗  decompose  Fi: cdiv∗/
(8)endWhile decomposing matrixFi, the update matrix Ui is also generated; specifically, the following equation presents the calculation of matrix Ui.(1)Ui=Fi2:r+1,2:r+1-li1,ili2,i⋮lir,i×li1,ili2,i⋯lir,i.


## 3. Vector Version and Architecture

This section chiefly introduces our vector version of the Cholesky algorithm and its implementation on FPGA.


### 3.1. Vector Version of Cholesky Algorithm

In Algorithm3, we rewrite multifrontal algorithm in a more detailed fashion. U(j) is a storage structure, the information of a series jth node in all the update matrices U. M(j) counts the number of its update matrices, which means the amount of child nodes of the jth node. p is an array that records the corresponding parent node. For instance, p(j) is the parent node of the jth node. In line (6) of Algorithm 3, the jth column of L is calculated by the first column of Fj. Lines (7) to (12) work on the generation of update matrix to its parent node and store it in U(p(j)). Notably, the update matrix is symmetric; thus, only a lower part of the matrix is generated and stored.Algorithm 3:Multifrontal Cholesky algorithm in a vector version.
Input: A
Output: L
(1)for  j=1; j≤n; j++  do
(2)Fj=ai,i⋯air,j⋮air,j0
(3)for  i=1; i≤M(j); i++  do
(4)Fj=Fj⊕Uj⋅Ui;
(5)end
(6)L(∗,j)=Fj(∗,1)/sqrt(Fj(1,1));
(7)M(p(j))++;
(8)for  k=j+1; k≤n; k++  do
(9)for t=k; t≤n; t++  do
(10)U(p(j))⋅UM(p(j))(t,k)=Fj(t,k)-L(t,j)×L(k,j);
(11)end
(12)end
(13)endThe execution of lines (2) to (3) in Algorithm 3 is illustrated in a vector fashion. In Algorithm 4 we exploited the strip-mined technology. z is the size of Fj while m is the amount of lanes in our vector architecture. k is the number of child nodes of this processed node. Furthermore, U(i;∗,t) is the tth column of update matrix from the ith child. Each updating of a column is assigned to a lane at one time. The calculation of L and the updating of matrix U are in a similar vector processing method.Algorithm 4:Updating algorithm of front matrix Fj.
Input: Fj, U
Output: Fj
(1)for  i=1; i≤kn; i++  do
(2)low=0;
(3)VL=mod(z,m);
(4)for  p=1; p≤z/m; p++  do
(5)for  t=low; t≤low+VL-1; t++  do
(6)Fj(∗,t)=Fj(∗,t)+U(i;∗,t);
(7)end
(8)      low = low +VL;
(9)VL = m;
(10)end
(11)end


### 3.2. Architecture

The primary architecture of vector machine containsk-independent processing elements (PEs). We divide cdiv and cmod modules for each PE. Figure 1 gives the details of the implementation. It is mainly composed of 3 computation modules and a series of storage components. The computation modules are Fet_num, Sqrt_Div, and Mul_Sub. The storage components include 5 FIFOs (First Input First Output, implemented by RAM) and RAMs for increasing buffer when data is transforming. Fet_num gets data from F_FIFO  (it stores frontal matrix which is Fj in Algorithm 3) and Uc_FIFO (it stores update matrix which is the output F_j in Algorithm 4) and combines with the data and the index that is located in data to finish extended addition operation. Then the operation in Sqrt_Div module starts up. The operation in Sqrt_Div is to finish the decomposition of column, which includes division and extraction of a root. The implementation of division and extraction of a root are based on SRT algorithm (it is usually called Digit Recurrence algorithm) and they share the implementation of SRT algorithm. Lt_FIFO stores the results of L(t,j). The output data from Sqrt_Div is divided into two parts: L(∗,j) is the final results and it is stored in L_FIFO. L(k,j) is stored in LRAM for updating. Data in LRAM is for Mul_Sub to generate update matrix. Mul_Sub includes multiplication and subtraction operation. Algorithm 4 is implemented in this module. The initial data in Cholesky factorization is obtained from DDR memory and the results after factorization are also stored in DDR.Figure 1
The implementation of vector architecture.The primary architecture of vector machine containsk-independent processing elements (PEs). We divide cdiv and cmod modules for each PE. Figures 2(a) and 2(b) show the two architectures. Each module has m lanes. From the two figures, we could conclude the data flow in one PE as described in Algorithm 3. For the extended addition between F and U, each lane processes one column at one time. More concretely, Figure 2(a) shows when mod(z,m)=2, the first two columns of F and corresponding columns in U are processed first. After that, all lanes are fully parallel most of the time. In Figure 2(b), only the first column of F needs to be processed when calculating the corresponding column of L. The rest of the columns in F need to be used when generating the updating matrix for its parent node.Figure 2
The architecture ofcmod and cdiv.
(a)
c
m
o
d
(b)
c
d
i
v


## 3.1. Vector Version of Cholesky Algorithm

In Algorithm3, we rewrite multifrontal algorithm in a more detailed fashion. U(j) is a storage structure, the information of a series jth node in all the update matrices U. M(j) counts the number of its update matrices, which means the amount of child nodes of the jth node. p is an array that records the corresponding parent node. For instance, p(j) is the parent node of the jth node. In line (6) of Algorithm 3, the jth column of L is calculated by the first column of Fj. Lines (7) to (12) work on the generation of update matrix to its parent node and store it in U(p(j)). Notably, the update matrix is symmetric; thus, only a lower part of the matrix is generated and stored.Algorithm 3:Multifrontal Cholesky algorithm in a vector version.
Input: A
Output: L
(1)for  j=1; j≤n; j++  do
(2)Fj=ai,i⋯air,j⋮air,j0
(3)for  i=1; i≤M(j); i++  do
(4)Fj=Fj⊕Uj⋅Ui;
(5)end
(6)L(∗,j)=Fj(∗,1)/sqrt(Fj(1,1));
(7)M(p(j))++;
(8)for  k=j+1; k≤n; k++  do
(9)for t=k; t≤n; t++  do
(10)U(p(j))⋅UM(p(j))(t,k)=Fj(t,k)-L(t,j)×L(k,j);
(11)end
(12)end
(13)endThe execution of lines (2) to (3) in Algorithm 3 is illustrated in a vector fashion. In Algorithm 4 we exploited the strip-mined technology. z is the size of Fj while m is the amount of lanes in our vector architecture. k is the number of child nodes of this processed node. Furthermore, U(i;∗,t) is the tth column of update matrix from the ith child. Each updating of a column is assigned to a lane at one time. The calculation of L and the updating of matrix U are in a similar vector processing method.Algorithm 4:Updating algorithm of front matrix Fj.
Input: Fj, U
Output: Fj
(1)for  i=1; i≤kn; i++  do
(2)low=0;
(3)VL=mod(z,m);
(4)for  p=1; p≤z/m; p++  do
(5)for  t=low; t≤low+VL-1; t++  do
(6)Fj(∗,t)=Fj(∗,t)+U(i;∗,t);
(7)end
(8)      low = low +VL;
(9)VL = m;
(10)end
(11)end


## 3.2. Architecture

The primary architecture of vector machine containsk-independent processing elements (PEs). We divide cdiv and cmod modules for each PE. Figure 1 gives the details of the implementation. It is mainly composed of 3 computation modules and a series of storage components. The computation modules are Fet_num, Sqrt_Div, and Mul_Sub. The storage components include 5 FIFOs (First Input First Output, implemented by RAM) and RAMs for increasing buffer when data is transforming. Fet_num gets data from F_FIFO  (it stores frontal matrix which is Fj in Algorithm 3) and Uc_FIFO (it stores update matrix which is the output F_j in Algorithm 4) and combines with the data and the index that is located in data to finish extended addition operation. Then the operation in Sqrt_Div module starts up. The operation in Sqrt_Div is to finish the decomposition of column, which includes division and extraction of a root. The implementation of division and extraction of a root are based on SRT algorithm (it is usually called Digit Recurrence algorithm) and they share the implementation of SRT algorithm. Lt_FIFO stores the results of L(t,j). The output data from Sqrt_Div is divided into two parts: L(∗,j) is the final results and it is stored in L_FIFO. L(k,j) is stored in LRAM for updating. Data in LRAM is for Mul_Sub to generate update matrix. Mul_Sub includes multiplication and subtraction operation. Algorithm 4 is implemented in this module. The initial data in Cholesky factorization is obtained from DDR memory and the results after factorization are also stored in DDR.Figure 1
The implementation of vector architecture.The primary architecture of vector machine containsk-independent processing elements (PEs). We divide cdiv and cmod modules for each PE. Figures 2(a) and 2(b) show the two architectures. Each module has m lanes. From the two figures, we could conclude the data flow in one PE as described in Algorithm 3. For the extended addition between F and U, each lane processes one column at one time. More concretely, Figure 2(a) shows when mod(z,m)=2, the first two columns of F and corresponding columns in U are processed first. After that, all lanes are fully parallel most of the time. In Figure 2(b), only the first column of F needs to be processed when calculating the corresponding column of L. The rest of the columns in F need to be used when generating the updating matrix for its parent node.Figure 2
The architecture ofcmod and cdiv.
(a)
c
m
o
d
(b)
c
d
i
v


## 4. Performance Model

Basing on the analysis of the Cholesky algorithm, we construct a performance model on FPGA platform to predict the performance of our designed vector architecture. The performance is related to the underlying architecture and its parameters. When ignoring the implementation architecture, we can predict the peak performance with on-chip memory, computational resources, andI/O bandwidth.Figure3 depicts the nonzeros pattern of L and its elimination tree. The elimination tree of the matrix L is defined to be the structure with size(L,1) nodes and node p is the parent of j if and only if the following is established.(2)p=min⁡i>j∣lij≠0.Figure 3
Nonzeros pattern ofL and its elimination tree.
(a)
(b)In the elimination tree, each column in the matrix is presented as a node. In Figure3 the numbers on the right side of the elimination tree represent the levels of the corresponding nodes. The path with a dashed line is the critical path. Given the existing dependency between parent node and its child node, the amount of nodes in the critical path, which is determined by the algorithm and the nonzeros pattern of the matrix, has the least number of nodes that should be processed sequentially.To construct the model on FPGA, we define these variables and their descriptions as listed in Table1.Table 1
Variables that we define and their descriptions.
AOs
The amount of operations for decomposing matrixA
RW
The data volume transferred between DDR3 and FPGA on-chip memory
BW
The available bandwidth between DDR3 and FPGA
OPV
Byte/Flop ratio, that isAOs/RW
T
c
d
i
v
The total computing time forcdiv of a PE
T
c
m
o
d
The total computing time forcmod of a PE
T
n
o
d
e
Computing time of decomposing one node
n
The amount of nonzeros in matrixL
Si
The amount of nodes in theith level of the elimination tree, i=1,2,…,n
nz
The amount of operations for decomposing matrixA
η
(
L
∗
j
)
The amount of the nonzero entries injth column of L
η
(
L
)
The amount of the nonzero entries ofLAssuming that data are initially stored in DDR3, the result will finally be written back to DDR3. Therefore, RW contains the input data, output data, and intermediate values. The intermediate values are forced to store and load in DDR3 when needed because of the small on-chip memory size. ThusRW=Input+output+2×temp. For a given problem, input and output data volumes are fixed, only temp is determined by the problem itself and FPGA on-chip memory.The amount of operations (AO) is determined by the amount of nonzero entries, specifically, forN×N matrix. It takes one square operation and η(L∗j)-1 division operation to scale one column (cdiv). Given that the update and front matrices are symmetric, only the lower part of the matrix should be taken into consideration. The amount of multiplication operations is 1/2η(L∗j)×η(L∗j)+1, and the AOs for addition in cmod are the same. As presented in the following equation, AO is the sum of all these operations.(3)AOs=∑j-1N1+ηL∗j-1+2×12×ηL∗j×ηL∗j+1=2ηL+∑j-1NηL∗j2.The chip resources contain on-chip memory and computational resources. The on-chip memory is the size of the memory that could be loaded immediately. It is a critical factor in RW. Moreover, computational resources determine the extent of parallel execution. In the next two subsections, we construct two different performance models for both memory- and computational-bounded cases.


### 4.1. Memory-Bound

Assuming FPGA owns enough computational resources and the system is totally memory-bound, the performance is determined by the memory access time; thus, the time of decomposition can be expressed as(4)TM=RWBW.For DDR3, when the data storage format and the access pattern are determined, BW will be determined. RW is a function of on-chip memory and the matrices. When memory bounded, the peak performance is described as(5)PM=AOsTM=AOsRW×BW.


### 4.2. Computational-Bound

Figure4(a) demonstrates the decomposition process of the matrix. The matrix is shown in Figure 3 with two PEs. We note that the bottleneck is shown in the first PE and we conclude that module cdiv and cmod can be executed in parallel; the total time is T=max{Tcdiv,Tcmod}.Figure 4
The parallel and pipelining execution.
(a)
The flow under the data dependency
(b)
The execution ofcdiv with m lanes in PEFrom Figure4(b), (clksqrt+clkdiv+⌈(η(L∗j)-2)/m⌉) is needed to finish cdiv when the operations are pipeline executed. Similarly, cmod will take (clkmult+clkadd+⌈(η(L∗j)+1)×η(L∗j)/2/m⌉) clocks. clksqrt,clkdiv,clkmult, and clkadd are the clocks that square root, division, multiplication, and addition operation take. In our implementation, our operation modules of square root, division, multiplication, and addition need 16, 27, 1, and 8 clocks, respectively. Generally, when η(L∗j)>8m, clkcmod>clkdiv, then we have (6)T=max⁡Tcdiv,Tcmod=Tcmod.Assuming that resources are enough,Tnode will be independent from the scale of nonzeros. They are executed in parallel and the time it takes to finish a square root, division, multiplication, and addition operation sequentially is the same. Besides that, only the time of nodes in the critical path will be counted into the amount of time. In this ideal case, the total time is TC1=∑j=1nclkcmod(ij) and the corresponding performance reaches the roofline value [21].To be more realistic, we assume one PE containsm lanes. A total of m nonzeros exist, executed in parallel after the diagonal element being processed in one node. In this case, we assume that one PE requires b units of resources. When the computational resources are bound and FPGA available resources are v units, one FPGA can contain k=⌊v/b⌋ PEs. It means that it supports k=⌊v/b⌋ nodes processed mostly in parallel.Assuming that the PE, which reaches bottleneck, is assigned to processi1,i2,…,ir nodes, we indicate that the amount of clocks is clks=∑j=1rclkcmod(ij) clocks. If the FPGA working frequency is f, the total process time is T=clks/f.Given that the largest number of nodes processed in parallel isk, decomposing the ith level needs time for processing ⌈Si/k⌉ nodes. For the whole elimination tree, a total of ∑i=1n⌈Si/k⌉(Nnode) nodes are processed sequentially at most. Consequently, when computational resources are bounded, the computation time is shown as (7)TC2=∑j=1Nnodeclkcmodijf.N
n
o
d
e
≥
n, which means that TC2≥TC1. So when computational resources and data dependency are taken into consideration, the process time is described as (8)TC=max⁡TC1,TC2=TC2.Furthermore, the performance is described as(9)PC=AOsTC.


### 4.3. Summary

When considering I/O bandwidth and computational resources, we can conclude that the following is the performance equation according to (4), (5), (6), (7), (8), and (9).(10)P=min⁡PC,PM=min⁡AOsTC,AOsTM=AOs×min⁡1TC,1TM.The analysis and summary are demonstrated in Table2. Once the parameters of FPGA and matrix are ascertained, we can obtain the peak performance and bottleneck.Table 2
The analysis of sparse Cholesky factorization.
Condition
Bottleneck
Performance
Conclusion
TC<TM
Chip memory, I/O bandwidth
OPV×BW
Performance is proportional to I/O bandwidth and Flop/Byte ratio
TC>TM
Computational resources
AOs×f∑j=1Nnodeclkcmodij
Performance is negative correlation toNnode and clkcmodSpecifically, when the system is computational-bound, the performance is determined byNnode and clkcmod. Table 3 discusses the influence of these two factors. From the table, these two factors are subject to the number of PE(k) and the number of lanes(m) in PE. Greater values of k and m do not guarantee a better performance. When k=max⁡(Si) and m=max⁡{(η(L∗j)+1)×η(L∗j)/2}, factor Nnode and clkcmod will reach their ultimate values.Table 3
The analysis of factors in computational-bound cases.
Factor
Calculation
Best value
Ultimate value
Nnode
N
n
o
d
e
=
∑
i
=
1
n
S
i
k
k=max⁡(Si)
Nnode=n
clkcmod
c
l
k
m
u
l
t
+
c
l
k
a
d
d
+
(
η
(
L
∗
j
)
+
1
)
×
η
L
∗
j
/
2
m
m
=
max
⁡
η
L
∗
j
+
1
×
η
L
∗
j
2
Tnode=clkmult+clkadd+1


## 4.1. Memory-Bound

Assuming FPGA owns enough computational resources and the system is totally memory-bound, the performance is determined by the memory access time; thus, the time of decomposition can be expressed as(4)TM=RWBW.For DDR3, when the data storage format and the access pattern are determined, BW will be determined. RW is a function of on-chip memory and the matrices. When memory bounded, the peak performance is described as(5)PM=AOsTM=AOsRW×BW.


## 4.2. Computational-Bound

Figure4(a) demonstrates the decomposition process of the matrix. The matrix is shown in Figure 3 with two PEs. We note that the bottleneck is shown in the first PE and we conclude that module cdiv and cmod can be executed in parallel; the total time is T=max{Tcdiv,Tcmod}.Figure 4
The parallel and pipelining execution.
(a)
The flow under the data dependency
(b)
The execution ofcdiv with m lanes in PEFrom Figure4(b), (clksqrt+clkdiv+⌈(η(L∗j)-2)/m⌉) is needed to finish cdiv when the operations are pipeline executed. Similarly, cmod will take (clkmult+clkadd+⌈(η(L∗j)+1)×η(L∗j)/2/m⌉) clocks. clksqrt,clkdiv,clkmult, and clkadd are the clocks that square root, division, multiplication, and addition operation take. In our implementation, our operation modules of square root, division, multiplication, and addition need 16, 27, 1, and 8 clocks, respectively. Generally, when η(L∗j)>8m, clkcmod>clkdiv, then we have (6)T=max⁡Tcdiv,Tcmod=Tcmod.Assuming that resources are enough,Tnode will be independent from the scale of nonzeros. They are executed in parallel and the time it takes to finish a square root, division, multiplication, and addition operation sequentially is the same. Besides that, only the time of nodes in the critical path will be counted into the amount of time. In this ideal case, the total time is TC1=∑j=1nclkcmod(ij) and the corresponding performance reaches the roofline value [21].To be more realistic, we assume one PE containsm lanes. A total of m nonzeros exist, executed in parallel after the diagonal element being processed in one node. In this case, we assume that one PE requires b units of resources. When the computational resources are bound and FPGA available resources are v units, one FPGA can contain k=⌊v/b⌋ PEs. It means that it supports k=⌊v/b⌋ nodes processed mostly in parallel.Assuming that the PE, which reaches bottleneck, is assigned to processi1,i2,…,ir nodes, we indicate that the amount of clocks is clks=∑j=1rclkcmod(ij) clocks. If the FPGA working frequency is f, the total process time is T=clks/f.Given that the largest number of nodes processed in parallel isk, decomposing the ith level needs time for processing ⌈Si/k⌉ nodes. For the whole elimination tree, a total of ∑i=1n⌈Si/k⌉(Nnode) nodes are processed sequentially at most. Consequently, when computational resources are bounded, the computation time is shown as (7)TC2=∑j=1Nnodeclkcmodijf.N
n
o
d
e
≥
n, which means that TC2≥TC1. So when computational resources and data dependency are taken into consideration, the process time is described as (8)TC=max⁡TC1,TC2=TC2.Furthermore, the performance is described as(9)PC=AOsTC.


## 4.3. Summary

When considering I/O bandwidth and computational resources, we can conclude that the following is the performance equation according to (4), (5), (6), (7), (8), and (9).(10)P=min⁡PC,PM=min⁡AOsTC,AOsTM=AOs×min⁡1TC,1TM.The analysis and summary are demonstrated in Table2. Once the parameters of FPGA and matrix are ascertained, we can obtain the peak performance and bottleneck.Table 2
The analysis of sparse Cholesky factorization.
Condition
Bottleneck
Performance
Conclusion
TC<TM
Chip memory, I/O bandwidth
OPV×BW
Performance is proportional to I/O bandwidth and Flop/Byte ratio
TC>TM
Computational resources
AOs×f∑j=1Nnodeclkcmodij
Performance is negative correlation toNnode and clkcmodSpecifically, when the system is computational-bound, the performance is determined byNnode and clkcmod. Table 3 discusses the influence of these two factors. From the table, these two factors are subject to the number of PE(k) and the number of lanes(m) in PE. Greater values of k and m do not guarantee a better performance. When k=max⁡(Si) and m=max⁡{(η(L∗j)+1)×η(L∗j)/2}, factor Nnode and clkcmod will reach their ultimate values.Table 3
The analysis of factors in computational-bound cases.
Factor
Calculation
Best value
Ultimate value
Nnode
N
n
o
d
e
=
∑
i
=
1
n
S
i
k
k=max⁡(Si)
Nnode=n
clkcmod
c
l
k
m
u
l
t
+
c
l
k
a
d
d
+
(
η
(
L
∗
j
)
+
1
)
×
η
L
∗
j
/
2
m
m
=
max
⁡
η
L
∗
j
+
1
×
η
L
∗
j
2
Tnode=clkmult+clkadd+1


## 5. Experimental Evaluation

The matrices used in this paper are part of the University of Florida matrix collection, which is involved in the field of electronics-magnetics, structure design, and so on. The data format is double floating-point (64 bits), needing 8 bytes per entry when using Compressed Sparse Column format and ignoring the index.


### 5.1. Accuracy of Performance Modeling

This paper implements a double floating-point version of sparse Cholesky algorithm on Xilinx Virtex-7 FPGA VC709 evaluation board (xc7vx485t-2ffg1761), which tightly integrates a high-performance ARM cortex with programmable logic devices. We implement this version in 4 different scales (1 PE, 1 module in 1 PE. 1 PE, 2 modules in 1 PE. 2 PE, 1 module in 1 PE. 2 PE, and 2 modules in 1 PE) to verify our performance module. The memory chip in our platform is 6.7 MB and the frequency is 200 MHz.We adopt the following steps in our experiments: (1) predicting the performance of our model using commonly used matrices, (2) measuring the actual performance with the matrices in our implementation, and (3) assessing the consistency of predicted and measured values. The results are shown in Figure 5 for accuracy evaluation and Figure 6 for comparison of performance different rates. In Figure 5, m is the number of PEs and k is the number of modules in one PE. The percentage of relative difference between a predicted value and measured value is calculated by (E-V)/V×100, where E is the predicted value and T is the measured value. There are many optimizations in hardware implementation (like DDR storage access mode, hardware configuration parameters). After optimizing the implementation and selecting some general class of matrices so they can get a better agreement, we have the following observations from our experimental data.Figure 5
Accuracy evaluation on different scales.
(a)
m
=
1; k=1
(b)
m
=
1; k=2
(c)
m
=
2; k=1
(d)
m
=
2; k=2Figure 6
Comparison of performance difference rates.The accuracy of model is critical since we need to use this model to infer large-scale parallel performance. We focus on performance and the difference rate of the performance. Figure5 shows the comparisons between the measured performance and predicted performance by our model of four classes matrix sets (shipsec1, nd24k, Trefethen_2000b, and nd3k) at different parallel scales (1 PE 1 module, 1 PE 2 module, 2 PE 1 module, and 2 PE 2 module). The measured performances are relatively stable and accurate. From [22–24] this established literature on performance model, we can conclude that the predicted and measured values have good consistency. Figure 6 presents the performance difference rates. As the scale of parallelism increased, the difference rates dropped form 10 percent to less than 5 percent.As shown in Figure5, the estimated performance is higher than the tested performance. We simplify data processing in the computational-bound model. Furthermore in the memory-bound model, a deviation in the bandwidth is found between DDR3 and FPGA.


### 5.2. Optimization

To obtain a balance in resources in FPGA and the performance, we set the memory chip to 64 MB, withm PEs and each PE having k modules.We compare the performance of our optimum implementation with four state-of-the-art supernode algorithms. One is the HSL_MA87 solver for multicore CPU platforms [15] while the other is the GPU version of CHOLMOD solver for CPU-GPU platforms [16]. The rest are subtree-based parallel methods for CPU-GPU platforms [17]. And the factorization time obtained is shown in Table 4. We have the following conclusions in the performance comparison.Table 4
Performance comparison with HSL_MA87, GPU version of CHOLMOD, and supernode.
Matrix
HSL_MA87 times (s)
CHOLMOD times (s)
One-GPU times (s)
Two-GPU times (s)
Ours
nd3k
2.02
2.92
1.27
1.12
1.96 (m=2, k=256)
nd24k
28.56
22.17
14.63
10.12
10.08 (m=8, k=32)
Trefethen_20000b
12.63
8.49
5.47
3.94
3.58 (m=4, k=64)The frequency of HSL_MA87, CHOLMOD, and GPU is 2.5 GHz, 3.2 GHz, and 2.93 GHz (CPU is 2.93 GHz, GPU is 1.15 GHz in [16]). Our frequency is 200 MHz. After testing matrix sets nd3k, nd24k, and Trefethen_20000b, our performance is on the same level as GPU performance. When our work is implemented as a coprocessor (i.e., when the frequency increases, our model can guide the design of coprocessor), our implementation has a greater advantage in performance.


## 5.1. Accuracy of Performance Modeling

This paper implements a double floating-point version of sparse Cholesky algorithm on Xilinx Virtex-7 FPGA VC709 evaluation board (xc7vx485t-2ffg1761), which tightly integrates a high-performance ARM cortex with programmable logic devices. We implement this version in 4 different scales (1 PE, 1 module in 1 PE. 1 PE, 2 modules in 1 PE. 2 PE, 1 module in 1 PE. 2 PE, and 2 modules in 1 PE) to verify our performance module. The memory chip in our platform is 6.7 MB and the frequency is 200 MHz.We adopt the following steps in our experiments: (1) predicting the performance of our model using commonly used matrices, (2) measuring the actual performance with the matrices in our implementation, and (3) assessing the consistency of predicted and measured values. The results are shown in Figure 5 for accuracy evaluation and Figure 6 for comparison of performance different rates. In Figure 5, m is the number of PEs and k is the number of modules in one PE. The percentage of relative difference between a predicted value and measured value is calculated by (E-V)/V×100, where E is the predicted value and T is the measured value. There are many optimizations in hardware implementation (like DDR storage access mode, hardware configuration parameters). After optimizing the implementation and selecting some general class of matrices so they can get a better agreement, we have the following observations from our experimental data.Figure 5
Accuracy evaluation on different scales.
(a)
m
=
1; k=1
(b)
m
=
1; k=2
(c)
m
=
2; k=1
(d)
m
=
2; k=2Figure 6
Comparison of performance difference rates.The accuracy of model is critical since we need to use this model to infer large-scale parallel performance. We focus on performance and the difference rate of the performance. Figure5 shows the comparisons between the measured performance and predicted performance by our model of four classes matrix sets (shipsec1, nd24k, Trefethen_2000b, and nd3k) at different parallel scales (1 PE 1 module, 1 PE 2 module, 2 PE 1 module, and 2 PE 2 module). The measured performances are relatively stable and accurate. From [22–24] this established literature on performance model, we can conclude that the predicted and measured values have good consistency. Figure 6 presents the performance difference rates. As the scale of parallelism increased, the difference rates dropped form 10 percent to less than 5 percent.As shown in Figure5, the estimated performance is higher than the tested performance. We simplify data processing in the computational-bound model. Furthermore in the memory-bound model, a deviation in the bandwidth is found between DDR3 and FPGA.


## 5.2. Optimization

To obtain a balance in resources in FPGA and the performance, we set the memory chip to 64 MB, withm PEs and each PE having k modules.We compare the performance of our optimum implementation with four state-of-the-art supernode algorithms. One is the HSL_MA87 solver for multicore CPU platforms [15] while the other is the GPU version of CHOLMOD solver for CPU-GPU platforms [16]. The rest are subtree-based parallel methods for CPU-GPU platforms [17]. And the factorization time obtained is shown in Table 4. We have the following conclusions in the performance comparison.Table 4
Performance comparison with HSL_MA87, GPU version of CHOLMOD, and supernode.
Matrix
HSL_MA87 times (s)
CHOLMOD times (s)
One-GPU times (s)
Two-GPU times (s)
Ours
nd3k
2.02
2.92
1.27
1.12
1.96 (m=2, k=256)
nd24k
28.56
22.17
14.63
10.12
10.08 (m=8, k=32)
Trefethen_20000b
12.63
8.49
5.47
3.94
3.58 (m=4, k=64)The frequency of HSL_MA87, CHOLMOD, and GPU is 2.5 GHz, 3.2 GHz, and 2.93 GHz (CPU is 2.93 GHz, GPU is 1.15 GHz in [16]). Our frequency is 200 MHz. After testing matrix sets nd3k, nd24k, and Trefethen_20000b, our performance is on the same level as GPU performance. When our work is implemented as a coprocessor (i.e., when the frequency increases, our model can guide the design of coprocessor), our implementation has a greater advantage in performance.



---

*Source: 3021591-2017-10-17.xml*