CODEWITHSUNDEEP
sparse-matrixchapel
Brian Dolan
Brian Dolan

Reputation: 3136

Custom rowSums of a Matrix in Chapel

Following up on this question. I have a Matrix (yes, I do) which will be large and sparse.

A = [
     [0,   0, 0, 1.2, 0]
     [0,   0, 0, 0,   0]
     [3.5, 0, 0, 0,   0]
     [0    7, 0, 0,   0]
]

And I want to create a vector v that has the sum v[j] = v[j,] * log(v[j,]) for each row in A

I believe there is an iterator like [x * log(x) for x in row] do... but I'm having a hard time finding the syntax. One particular bugaboo is to avoid taking log(0), so maybe an if statement in the iterator?

Upvotes: 1

Views: 53

Answers (1)

ben-albrecht
ben-albrecht

Reputation: 1865

I believe there is an iterator like [x * log(x) for x in row] do... but I'm having a hard time finding the syntax.

Instead of creating an iterator, we can create a function to compute x * log(x) and just pass an array (or array slice) to it, allowing promotion to take care of the rest.

Instead of doing a + reduce over the array slices like we did before,

forall i in indices {
  rowsums[i] = + reduce(A[i, ..]);
}

We can do a + reduce over a promoted operation on the array slices, like this:

forall i in indices {
  rowsums[i] = + reduce(logProduct(A[i, ..]));
}

where logProduct(x) can include an if-statement handling the special case of 0, as you mentioned above.

Putting this all together looks something like this:

config const n = 10;

proc main() {
  const indices = 1..n;
  const Adom = {indices, indices};
  var A: [Adom] real;

  populate(A);

  var v = rowSums(A);

  writeln('matrix:');
  writeln(A);

  writeln('row sums:');
  writeln(v);
}


/* Populate A, leaving many zeros */
proc populate(A: [?Adom]) {
  forall i in Adom.dim(1) by 2 do // every other row
    forall j in Adom.dim(2) by 3 do // every third column
      A[i, j] = i*j;
}

/* Compute row sums with custom function (logProduct) */
proc rowSums(A: [?Adom] ?t) {
  var v: [Adom.dim(1)] t;

  [i in v.domain] v[i] = + reduce(logProduct(A[i, ..]));

  return v;
}

/* Custom function to handle log(0) case */
proc logProduct(x: real) {
  if x == 0 then return 0;
  return x * log(x);
}

Upvotes: 1

Related Questions