### Symmetric Sparse Matrix

Sparse matrices are the matrices populated primarily with zeros. Iterative linear solvers, such as the Lanczos-Arnoldi procedure, are among the most important computational tools for solving large symmetric (Hermitian) sparse linear problems. The fundamental operation in the iterative methods is the sparse matrix vector multiplication (SpMV).

Sparse matrix can be stored in several formats. The most flexible one is the coordinate storage scheme (COO). In this scheme, one use three arrays (two integral and one scalar):

###

###

In this scheme, the $n_p$ processors are aligned as a $n_d \times n_d$ 2D grid. For symmetric matrices, only half of the grid is needed. Then, for $n_p$ processors, $ n_d (n_d + 1)/2 \le n_p \implies n_d \le \text{floor}\left( \frac{\sqrt{1+n_p}-1}{2}\right)$. It is useful to assign each processor a 2D coordinate $(I, J)$, where $I,J = 1,2,3,\cdots, n_d$. We will discuss the processor assignment later.

The vector $x$ is stored on the diagonal and the matrix $A$ is distributed across the processor grid. More specifically, vector element $x_i$ is stored in processor $(I, I)$, where $I = \mod(i, n_d) + 1$. Matrix element $A_{ij}$ is distributed to processor $(I, J)$, $I = \mod(i, n_d)+1, J = \mod(j,n_d)+1$, forming a local block matrix $A_{IJ}$.

The local blocks need not be square matrix, except the diagonal blocks. Let $n_I$ be the dimensionality of the $I$-th diagonal blocks. It is easy to see, the row and column dimensionality of the block $A_{IJ}$ is $n_I$ and $n_J$, respectively. \[ n_I = \text{floor}\left(\frac{n-1}{n_d}\right) + \vartheta\left(n-\text{floor}\left(\frac{n-1}{n_d}\right)\cdot n_d - I \ge 0\right) \]

Matrix block $A_{IJ}$ is on processor $ p = F(I-J-1)+J-1$, where $F(x) = (x+1) (2 n_d - x)/2$. In practice, one knows the local process id $p$ and wants to determine which matrix block is stored on it: \[

\begin{split}

I &= 1+p+1+G(p+1)-F(G(p+1)-1); \\

J &= p+1-F(G(p+1)-1).

\end{split}

\] where $G(x) = \text{floor}\left(n_d - \sqrt{\left(n_d+\frac{1}{2}\right)^2-2x} - \frac{1}{2}\right)$.

To store the node information, we introduce a FORTRAN module (a struct is a similar structure in C):

The problem of the DM scheme is that the number of processors in each row & column communicator is not the same. Indeed, it is not well balanced: some rows & columns have $n_d$ (ndiag) processors, some of them have only 1. This problem can be solved by using the so-called Balanced Column Major scheme.

In the BCM scheme, each row or column has $\frac{n_d+1}{2}$ blocks.

The corresponding matrix blocks distribution on processor $p$ is,

if $ p\le \text{ceiling}\left(\frac{n_d+1}{2}\right)\left[ n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]-1$ then \[ \begin{split}

J &= \text{floor}\left(\frac{p}{\text{ceiling}\left(\frac{n_d+1}{2}\right)}\right)+1; \\ I &= \mod\left( p, \text{ceiling}\left(\frac{n_d+1}{2}\right)\right)+J;

\end{split} \]

else,

\[ \begin{split} J &= n_d+2-\text{ceiling}\left(\frac{n_d+1}{2}\right)+\text{floor}\left( \frac{q}{\text{floor}\left(\frac{n_d+1}{2}\right)}\right); \\ I &= \mod\left( \mod\left(q, \text{floor}\left(\frac{n_d+1}{2}\right)\right) +J-1, n_d\right)+1; \end{split} \] where $q = p-\text{ceiling}\left(\frac{n_d+1}{2}\right)\left[n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]$.

The nodeinfo.mod can be modified accordingly.

\[ A x + y \to y \Rightarrow T x + (x^t L)^t + y \to y \]

Sparse matrix can be stored in several formats. The most flexible one is the coordinate storage scheme (COO). In this scheme, one use three arrays (two integral and one scalar):

row(nnz), col(nnz), mat(nnz)to store the row indices, column indices and the values of the nonzero matrix elements. It's easy to convert a coordinate storage scheme into other schemes such as the compressed column/row scheme (CCS/CRS). For SpMV, coordinate scheme is as efficient as CCS/CRS. It demands $2 n_{nz} + n_{nz}$ memory comparing with CCS/CRS $(n_{nz} + n) + n_{nz}$, where $n$ is the dimensionality of the matrix and $n_{nz}$ is the number of nonzeros, and $n^2 \gg n_{nz} \gg n$ for most useful cases. Only half of the matrix is stored, because the matrix is symmetric (Hermitian).

### Serial Algorithm

subroutine SpMV(nnz, row, col, mat, n, x, y, ierr) implicit none ! A*x + y -> y ! n is the dimension of the vector x, y ! nnz is the number of nonzeros integer, intent(in) :: nnz, n, row(nnz), col(nnz) double precision, intent(in) :: mat(nnz), x(n) double precision, intent(inout) :: y(n) integer, intent(out) :: ierr integer :: i do i = 1, nnz y(row(i)) = y(row(i)) + mat(i)*x(col(i)) if(row(i) .ne. col(i)) then y(col(i)) = y(col(i)) + x(row(i))*conjg(mat(i)) end if end do ierr = 0 end subroutine

### Parallel Algorithm

Here we only talk about MPI parallelism. In reality, it is most efficient to use both MPI and OpenMP. The OpenMP parallelism can be used in the series part.

#### 2D cyclic distribution

A parallel SpMV depends on the the distribution of the processors across the processors. 1D distribution has a serious load balancing problem. An simple and popular scheme is the 2D cyclic distribution.Left: A naive 2D distribution. Center: A cyclic 2D distribution. Right: A cyclic 2D distribution with non-square block matrices. |

The vector $x$ is stored on the diagonal and the matrix $A$ is distributed across the processor grid. More specifically, vector element $x_i$ is stored in processor $(I, I)$, where $I = \mod(i, n_d) + 1$. Matrix element $A_{ij}$ is distributed to processor $(I, J)$, $I = \mod(i, n_d)+1, J = \mod(j,n_d)+1$, forming a local block matrix $A_{IJ}$.

The local blocks need not be square matrix, except the diagonal blocks. Let $n_I$ be the dimensionality of the $I$-th diagonal blocks. It is easy to see, the row and column dimensionality of the block $A_{IJ}$ is $n_I$ and $n_J$, respectively. \[ n_I = \text{floor}\left(\frac{n-1}{n_d}\right) + \vartheta\left(n-\text{floor}\left(\frac{n-1}{n_d}\right)\cdot n_d - I \ge 0\right) \]

#### Processor Assignment: Diagonal Major (DM) Scheme

In the DM scheme, the processors are assigned to the matrix blocks according to their diagonal coordinate:\begin{split}

I &= 1+p+1+G(p+1)-F(G(p+1)-1); \\

J &= p+1-F(G(p+1)-1).

\end{split}

\] where $G(x) = \text{floor}\left(n_d - \sqrt{\left(n_d+\frac{1}{2}\right)^2-2x} - \frac{1}{2}\right)$.

To store the node information, we introduce a FORTRAN module (a struct is a similar structure in C):

module nodeinfo implicit none integer :: nproc, myrank, myrow, mycol, ndiag, idiag; integer :: proc_comm, row_comm, col_comm; contains subroutine initnodeinfo use mpi implicit none integer :: ierr ! determine ndiag, nproc and rank call mpi_comm_size(mpi_comm_world, nproc, ierr) call mpi_comm_rank(mpi_comm_world, myrank, ierr) ndiag = floor(0.5*sqrt(1.+8.*nproc)-0.5) color = 0 if(ndiag*(ndiag+1)/2>rank) then color = 1 endif call mpi_comm_split(mpi_comm_world, color, myrank, proc_comm, ierr) call mpi_comm_size(proc_comm, nproc, ierr) call mpi_comm_rank(proc_comm, myrank, ierr) ! determine myrow and mycol i.e. (I,J) the block ID mycol = 1+rank-F(G(myrank,ndiag),ndiag) myrow = G(myrank,ndiag)+mycol+1 idiag = myrow - mycol ! create row_comm and col_comm, diag_comm call mpi_comm_split(proc_comm, myrow, abs(idiag), row_comm, ierr) call mpi_comm_split(proc_comm, mycol, abs(idiag), col_comm, ierr) call mpi_comm_split(proc_comm, idiag, myrank, diag_comm, ierr) end subroutine integer function F(x, nd) implicit none integer, intent(in) :: x, nd F = (x+1)*(2*nd-x)/2 end function integer function G(x, nd) implicit none integer, intent(in) :: x, nd G = floor(nd-0.5-sqrt((nd+0.5)**2-2*x)-0.5) end function integer function nloc(n, nd, i) implicit none integer, intent(in) :: n, nd, i nloc = (n-1)/nd if(n.ge.((n-1)/nd)*nd+i) then nloc = nloc + 1 end if end function end module

The problem of the DM scheme is that the number of processors in each row & column communicator is not the same. Indeed, it is not well balanced: some rows & columns have $n_d$ (ndiag) processors, some of them have only 1. This problem can be solved by using the so-called Balanced Column Major scheme.

#### Processor Assignment: Balanced Column Major (BCM) Scheme

In the BCM scheme, each row or column has $\frac{n_d+1}{2}$ blocks.

The corresponding matrix blocks distribution on processor $p$ is,

if $ p\le \text{ceiling}\left(\frac{n_d+1}{2}\right)\left[ n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]-1$ then \[ \begin{split}

J &= \text{floor}\left(\frac{p}{\text{ceiling}\left(\frac{n_d+1}{2}\right)}\right)+1; \\ I &= \mod\left( p, \text{ceiling}\left(\frac{n_d+1}{2}\right)\right)+J;

\end{split} \]

else,

\[ \begin{split} J &= n_d+2-\text{ceiling}\left(\frac{n_d+1}{2}\right)+\text{floor}\left( \frac{q}{\text{floor}\left(\frac{n_d+1}{2}\right)}\right); \\ I &= \mod\left( \mod\left(q, \text{floor}\left(\frac{n_d+1}{2}\right)\right) +J-1, n_d\right)+1; \end{split} \] where $q = p-\text{ceiling}\left(\frac{n_d+1}{2}\right)\left[n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]$.

The nodeinfo.mod can be modified accordingly.

#### The Parallel SpMV Algorithm Based on the 2D cyclic Distribution

$A = L + D + L^t$, where $L$ is a strict lower triangular matrix and $D$ is a diagonal matrix. $T = L+D$.\[ A x + y \to y \Rightarrow T x + (x^t L)^t + y \to y \]

subroutine SpMV(nnz, row, col, mat, ncol, x, yt, nrow, y, xt, ierr) use nodeinfo use mpi implicit none ! A*x + y -> y ! A = L + D + L^t, where L is a strict lower triangular matrix, ! D is a diagonal matrix, L^t is the transpose of L. ! Let T = L+D; A = T + L^t. ! A*x = T*x + L^t*x = T*x + (x^t*L)^t ! ! x/y is the local vector. Only x/y on the diagonal are initialized/output ! xt/yt is the local vector of the transpose vector of x/y ! note their local dimensionalities nrow, ncol are not the same as the original vector. integer :: nnz, ncol, nrow, ierr, row(nnz), col(nnz); double precision :: x(ncol), yt(ncol), y(nrow), xt(nrow), mat(nnz); intent(in) :: nnz, ncol, nrow, row, col, mat; ! (L+D)*x + y -> y call mpi_bcast(x, ncol, mpi_real8, 0, col_comm, ierr) ! local SpMV if(idiag.eq.0) then xt = y end if do i = 1, nnz xt( 1+(row(i)-1)/ndiag ) = xt( 1+(row(i)-1)/ndiag ) + mat(i)*x( 1+(col(i)-1)/ndiag ) end do call mpi_reduce(xt, y, nrow, mpi_real8, mpi_sum, 0, row_comm, ierr) ! x^t * L -> y^t if(idiag.eq.0) then xt = x x = y end if call mpi_bcast(xt, nrow, mpi_real8, 0, row_comm, ierr) ! local SpMV do i = 1, nnz if( row(i) .eq. col(i) ) then x( 1+(col(i)-1)/ndiag ) = x( 1+(col(i)-1)/ndiag ) + xt( 1+(row(i)-1)/ndiag )*conjg(mat(i)) end if end do call mpi_reduce(x, yt, ncol, mpi_real8, mpi_sum, 0, col_comm, ierr) if(idiag.eq.0) then y = y + yt x = xt endif end subroutine

## No comments:

## Post a Comment