Indirect addressing

Indirect addressing (pointers) are very common in codes for sparse matrices (PDE-problems) for example. In this assignment you will study the performance of a sparse daxpy-operation.

Let x and y be two double precision vectors with n = 1000000 elements each. Let p an integer vector with indices, a is a double precision number. We would like to form the following sparse daxpy-operation:

do k = 1, n
  j = p(k)
  y(j) = y(j) + a * x(j)
end do

or in C

for(k = 0; k < n; k++) {
  j = p[k];
  y[j] += a * x[j];
}

In our test the p-vector will have the same length as x and y, but in an application it could be shorter. You will test two different values of p and compare with using no indirect addressing.

Implement sparse_daxpy(n, a, x, y, p) (the code above) and daxpy(n, a, x, y) (write your own daxpy; if you are using the standard daxpy you must use the increment parameters in the call, daxpy(n, a, x, 1, y, 1)) and write the following main program:

• initialize x, y and a.
   
• initialize p the first way and call sparse_daxpy 100 times; measure the elapsed time.
   
• initialize p the second way and call sparse_daxpy 100 times; measure the elapsed time.
   
• call daxpy 100 times; measure the elapsed time.

Use full optimization in all cases.

You may want to read about large arrays in C first. See the C-intro. under the Diary.

You should be able to fill out the following table when you ready:

Comments?