Next: Performance for sparse matrix Up: The Sparse Matrix Benchmark Previous: I/O considerations

## The Computational Kernel

Once the matrix is read in, it is converted into the internal data structure. Here the popular compacted sparse row format for storing sparse matrices is used. This comprises two integer arrays and one double precision array , where is of dimension and points to the start of each row in the arrays and . The array is of dimension and contains the column indices, while is also of dimension and contains the matrix entries.

The main kernel in the computational part of the benchmark is the code for the product of two matrices. In F90 this is of the following form, where the product is calculated. Each matrix is represented in the compacted sparse row format; for example, matrix is represented by arrays , matrix by arrays , and product matrix by arrays .

```
subroutine spmatmul(n,m,ia,ja,a,ib,jb,b,ic,jc,c)
! C = A*B
integer (kind = myint), intent (in) :: n,m
integer (kind = myint), intent (in) :: ia(*),ib(*),ja(*),jb(*)
integer (kind = myint), intent (out) :: ic(*),jc(*)
real (kind = myreal), intent (in) :: a(*),b(*)
real (kind = myreal), intent (out) :: c(*)
real (kind = myreal) :: aij
! initialise the mask array which is an array that has
! non-zero value if the column index already exist, in which
! case the value is the index of that column
nz = 0
ic(1) = 1
do i=1,n
do j = ia(i),ia(i+1)-1
aij = a(j)
neigh = ja(j)
do k = ib(neigh),ib(neigh+1)-1
if (icol == 0) then
nz = nz + 1
c(nz) =  aij*b(k)
else
c(icol) = c(icol) + aij*b(k)
end if
end do
end do
ic(i+1) = nz+1
end do
end subroutine spmatmul
```

The main characteristic of the above computation is that although the arrays are accessed sequentially, the arrays and may be accessed randomly through indirect addressing.

In Java, the kernel looks very similar:

```
static void spmatmul_double(
int n, int m,
double[] a, int[] ia,int[] ja,
double[] b, int[] ib,int[] jb,
double[] c, int[] ic,int[] jc)
{
int nz;
double aij;
int neighbour;

// extra space for FORTRAN like array indexing

// starting from one for FORTRAN like array index
for (l = 1; l <= m; l++) mask[l] = 0;

ic[0] = 1;
nz = 0;

// starting from one for FORTRAN like array indexing
for (i = 1; i <= n; i++) {
for (j = ia[i]; j < ia[i+1]; j++){
neighbour = ja[j];
aij = a[j];
for (k = ib[neighbour]; k < ib[neighbour+1]; k++){
if (icol ==  0) {
c[nz] =  aij*b[k];
}
else {
c[icol] +=  aij*b[k];
}
}
}
for (j = ic[i]; j < nz + 1; j++) mask[jc[j]] = 0;
ic[i+1] = nz+1;
}

}
```

Notice that because the matrix is stored such that the row and column entries always start from , in Java each array will be given one extra element so that the array addressing can start from . In C however it is possible to utilize the flexibility of pointer manipulation to avoid this extra element, as the following C code shows.

```
void spmatmul_double(
const int* n, const int* m,
const double *a, const int *ia,const int *ja,
const double *b, const int *ib,const int *jb,
double *c, int *ic,int *jc,
const int* ind_base)
{
int nz;
const double *aij;
const int *neighbour;

mask = (int*) malloc( (*m) * sizeof(int));

for (i = 0; i != *m; i++) *pmask++ = 0;
/* shift the index base for fortran like indexing: *ind_base=1 */
ib -= *ind_base;
jb -= *ind_base;
b -= *ind_base;
jc -= *ind_base;
c -= *ind_base;

aij = a;
neighbour = ja;
nz = -1+*ind_base;
*ic = 1;
for (i = 0; i != *n; i++) {
for (j = ia[i]; j != ia[i+1]; j++){
for (k = ib[*neighbour];k != ib[*neighbour+1];k++){
if (icol ==  0)
{
c[nz] =  (*aij)*b[k];
}
else
{
c[icol] +=  (*aij)*b[k];
}
}
aij++;
neighbour++;
}
for (j = *ic; j != nz + 1; j++) mask[jc[j]] = 0;
*(++ic) = nz+1;
}