The invention generally relates to a technique for transposing nonsymmetric sparse matrices.
Sparse matrix representations are common throughout physical modeling, machine learning and formulations of problems as differential equations. A sparse matrix is a matrix that has a significant number of zeros, as compared to the dimension of the matrix. Due to the large number of zeros, a sparse matrix may be represented in a compressed format for purposes of reducing the, storage and work that is required in processing the matrix.
It is common for computer applications to generate and use the transpose of a matrix. A matrix may be transposed by interchanging its rows and columns. A symmetric matrix is a matrix in which the original matrix and the transposed matrix are identical. Conventional techniques to transpose a non-sparse matrix may be relatively inefficient when applied to transposing a non-symmetric sparse matrix.
In accordance with embodiments of the invention described herein, parallel processing is used to transpose a non-symmetric sparse matrix (also called a “source matrix” herein). Due to the parallel processing, the time to transpose the source matrix is greatly reduced, as compared to serial processing; and additionally, the technique is scalable with future generations of parallel processing architecture.
In the following description, the source matrix is represented in a compressed format. More specifically, in a compressed sparse row (SCR) representation, the source matrix is completely described using three arrays. The first array, called the “Nonzeros array,” sets forth the nonzeros of the source matrix. Another array, called the “ColumnIndices array,” indicates the specific column positions of the nonzeros. A third array, called the “Rows array,” indicates the first nonzero of each of row. Alternatively, the Rows array may indicate the size of each row of the source matrix.
For purposes of clarifying the following discussion, the following exemplary nonsymmetric sparse matrix (called the “exemplary source matrix” herein) to illustrate operations that are used to transpose the matrix:
It is noted that other non-symmetric sparse matrices may be processed according to the techniques described herein, in accordance with embodiments of the invention. For the exemplary source matrix, the transpose is as follows:
Pursuant to the CSR representation format, the nonsymmetric sparse matrix may be described as follows:
For this example, the Nonzeros array sets forth all nine nonzeros of the exemplary source matrix. Due to the nine nonzeros, the Nonzeros array has a length of “9.” The ColumnIndices array, of the same length as the Nonzeros array, sets forth the column position of each nonzero of the Nonzeros array. For example, assuming the array elements in the CSR representation are indexed from 0 to 8, element two of the ColumnIndices array, corresponds to element two, a “3,” of the Nonzeros array. This means that the “3” in the exemplary source matrix is located in column number two of the matrix. Similarly, element eight of the ColumnIndices array, a “2,” indicates that element eight of the Nonzeros array, a “9,” is located in column two of the exemplary source matrix.
Pursuant to the CSR representation, the Rows array denotes the first nonzero element in each row. In this regard, the Rows array includes five elements, each of which corresponds to a row (five total) of the exemplary source matrix. Each element of the Rows array is associated with a particular row and indicates the first nonzero element for that row. For example, element three of the Rows array, a “6,” is associated with row three of the sparse matrix and indicates that element six in the Nonzeros array, a “7,” is the first nonzero element in row three.
The CSR representation is an example of a compressed representation for a sparse matrix in accordance with some embodiments of the invention. However, it is noted that other compressed representations may be used in other embodiments of the invention. For example, a compressed sparse column (CSC) representation may alternatively be used in other embodiments of the invention. The exemplary source matrix has the following CSC representation:
The main difference between the CSR and CSC representations is that the nonzeros are in column order in the CSC representation. The CSC representation contains a Nonzeros array, which indicates the nonzeros of the sparse matrix, and the CSC representations also contains an array called the “RowIndices array,” which indicates the row positions of each corresponding nonzero in the Nonzeros array. Thus, element five of the RowIndices array, a “4,” indicates that element five of the Nonzeros array, a “9,” is located in row four of the exemplary sparse matrix. In the CSC representation, the RowIndices array replaces the ColumnIndices array of the CSR representation. The CSC representation also includes an array called the “Cols array,” which indicates the first nonzero element in each column. For example, element one of the Cols array, a “1,” indicates that a “1,” which is at element one in the Nonzeros array, is the first element in column one of the sparse matrix.
It is noted that the CSR representation of the transpose of the source matrix is identical to the CSC representation of the source matrix (accounting for any change of naming convention).
For purposes of simplifying the discussion herein, the CSR representation is assumed, although other representations (such as the CSC representation) may be used in other embodiments of the invention, as these representations fall within the scope of the appended claims.
As described further below, the computation of the transpose relies on several primitives, each of which may be executed in parallel. These primitives may be viewed as composible patterns of parallel computation, so that the technique may be readily applied in designing a lower level, explicitly threaded and/or vectorized implementation.
In the discussion below, the number of nonzeros in the source matrix is referred to as. “n,” the number of processors that are used to perform the transpose is referred to as “p,” and the number of rows in the source matrix is referred to as “r.” It is also assumed, for purposes of simplifying the discussion herein, that the source matrix is square, although non-square matrices may be processed in accordance with other embodiments of the invention.
More specifically, in accordance with some embodiments of the invention, the processors 40 receive a CSR representation of the source matrix, a representation which includes the Nonzeros array, ColumnIndices array and the Rows array; and the processors 40 execute primitives 50 to generate a CSR representation of the transpose. The CSR representation of the transpose is associated with a newNonzeros array, a newColumnIndices array, and a newRows array, which correspond to the Nonzeros array, ColumnIndices array and Rows array, respectively, of the CSR representation of the source matrix. The computation of the newNonzeros array, newColumnIndices array and newRows array are set forth below for a specific example. In accordance with some embodiments of the invention, the processors 40 execute the following primitives 50, which are described in more detail below: an element wise add primitive 50a, a prefixsum primitive 50b, a multireduce primitive 50c, a multiprefix primitive 50d, a safe forward permute primitive 50e and a safe backward permute primitive 50f.
Referring to
The processors 40 execute the prefixSum primitive 50b to sum all of the elements of the array are summed, while keeping all partial sums. Serially, the prefixSum primitive 50b may be described by the following pseudo code:
In the above pseudo code, the length of each of the C and A matrices is “n.” The prefixSum primitive 50b is executed in parallel by the multiple processors 40, an execution in which p threads are spawned off, and where each thread is responsible for computing a portion of the overall computation. This computation proceeds in three logical phases. The pseudo code, which describes the parallel execution of the prefixSum primitive 50b by the processors 40 is set forth below:
The total time to execute p threads in parallel is (n/p−1)+log2(p)+(n/p)=2*(n/p)+log2(p)−1. If the length of n is much greater than p, then the term n/p dominates and thus, scales linearly with p.
Regarding the multiReduce primitive 50c (see also
In the above pseudo code, the matrix A has a length of “n,” the matrix C has a length of “m,” and the indices in the matrix A fall within the range 0 . . . m. The processors 40 execute the multiReduce primitive 50c in parallel by spawning off p threads where each thread is responsible for computing a portion of the overall computation. The overall computation by the processors proceeds in three logical phases, which is described by the following pseudo code:
The total time to execute the multiReduce primitive 50c with p threads in parallel is (n/p)+m*log 2(p). If the length of n is much greater than m*log 2(p), then the term n/p dominates and this scales linearly with p. If the length of n is not much greater than m*log 2(p), then a filter may be used on the second phase of the computation that only updates those components of partialSums[t] that were last updated. This reduces the complexity to at most m, so the requirement relaxes to n/p being greater than m, which is typical.
The multiprefix primitive 50d is similar to the multiReduce primitive 50c described above, except that the multiprefix primitive 50d keeps track of the partial count at each element of the matrix A. Serially, the multiprefix primitive 50d may be described by the following pseudo code:
In the above pseudo code, the matrix A has a length of “n,” the matrix C has a length of “m,” and the indices in the matrix A fall within the range 0 . . . m. The processors 40 perform the multiPrefix primitive 50d in parallel, in which p threads are spawned off and each is responsible for computing a portion of the overall computation. The overall computation proceeds in three logical phases and may be described by the following pseudo code:
The total time complexity here if p threads are executed in parallel is (n/p)+m*log 2(p)+(n/p=2*(n/p)+m*log 2(p). If the length of n is much greater than m*log 2(p), then the term n/p dominates and this scales linearly with p. If the length of n is not much greater than (n/p)+m*log 2(p), then we can use a filter on the second phase of the computation that only updates those components of partialSums[t] that were last updated. This reduces the complexity to at most m, so the requirement relaxes to n/p being greater than m, which is typical.
Referring to
More specifically, the generation of the transpermute vector may be described by the following pseudo code, which uses the exemplary source matrix for a more specific example:
The eltwiseAdd primitive 50a (see
In the pseudo code that is set forth above, the matrices A, B and C each have a length of “n.” The length of the processors 40 execute the eltwiseAdd primitive 50a in parallel, as p threads are spawned off, and each thread is responsible for computing a portion of the overall computation of size n/p. This means that the algorithm scales linearly with the number of processors 40 since the time to complete across the threads is n/p. The parallel computation may be described by the following pseudo code:
Referring to
The pseudo code related to computing the newColumnIndices, with the specific example given for the exemplary source matrix, is set forth below:
The safeForwardPermute primitive 50e is described as follows. The length of each of the A and Ind matrices is n, the indices in Ind fall within the range of 0 . . . (m−1), and the length of the C matrix is m. Serially, the safeForwardPermute primitive may be described as follows:
The processors 40 execute the safeForwardPermute primitive 50e in parallel by spawning off p threads, where each thread is responsible for computing a portion of the overall computation of size n/p. This means that the algorithm scales linearly with the number of processors since to complete across the threads is n/p. Because Ind has no repetitions, each write to the C matrix is guaranteed not to conflict. The parallel implementation of the safeForwardPermute primitive may be described by the following pseudo code:
The safeBackwardPermute primitive 50f (see
The processors 40 execute the safeBackwardPermute primitive 50f in a parallel fashion by spawning off p threads, where its thread is responsible for computing a portion of the overall computation of size n/p. This means that the algorithm scales linearly with the number of processors since the time to completion across the threads is n/p.
Pseudo code, which describes the parallel computations, may be described as follows:
Referring to
In terms of pseudo code, the processors 40 perform the following operation with a specific example being given for the exemplary source matrix that is set forth herein:
While the invention has been disclosed with respect to a limited number of embodiments, those skilled in the art, having the benefit of this disclosure, will appreciate numerous modifications and variations therefrom. It is intended that the appended claims covet all such modifications and variations as fall within the true spirit and scope of the invention.