Constructing a Sparse Tropical Limit Function in Chapel

Question

Given matrices A and B the tropical product is defined to be the usual matrix product with multiplication traded out for addition and addition traded out for minimum. That is, it returns a new matrix C such that,

C_ij = minimum(A_ij, B_ij, A_i1 + B_1j, A_i2 + B_12,..., A_im + B_mj)

Given the underlying adjacency matrix A_g of a graph g, the nth "power" with respect to the tropical product represents the connections between nodes reachable in at most n steps. That is, C_ij = (A**n)_ij has value m if nodes i and j are separated by m<=n edges.

In general, given some graph with N nodes. The diameter of the graph can only be at most N; and, given a graph with diameter k, A**n = A**k for all n>k and the matrix D_ij = A**k is called the "distance matrix" entries representing the distances between all nodes in the graph.

I have written a tropical product function in chapel and I want to write a function that takes an adjacency matrix and returns the resulting distance matrix. I have tried the following approaches to no avail. Guidance in getting past these errors would be greatly appreciated!

proc tropicLimit(A:[] real,B:[] real) {
 var R = tropic(A,B);
 if A == R {
   return A;
 } else {
   tropicLimit(R,B);
 }
}

which threw a domain mismatch error so I made the following edit:

proc tropicLimit(A:[] real,B:[] real) {
 var R = tropic(A,B);
 if A.domain == R.domain {
   if && reduce (A == R) {
     return R;
   } else {
     tropicLimit(R,B);
   }
 } else {
   tropicLimit(R,B);
 }
}

which throws

src/MatrixOps.chpl:602: error: control reaches end of function that returns a value

proc tropicLimit(A:[] real,B:[] real) {
 var R = tropic(A,B);
 if A.domain == R.domain {
   if && reduce (A == R) {  // Line 605 is this one
   } else {
     tropicLimit(R,B);
   }
 } else {
   tropicLimit(R,B);
 }
return R;
}

Brings me back to this error

src/MatrixOps.chpl:605: error: halt reached - Sparse arrays can't be zippered with anything other than their domains and sibling arrays (CS layout)

I also tried using a for loop with a break condition but that didn't work either

proc tropicLimit(B:[] real) {
 var R = tropic(B,B);
 for n in B.domain.dim(2) {
   var S = tropic(R,B);
   if S.domain != R.domain {
    R = S; // Intended to just reassign the handle "R" to the contents of "S" i.o.w. destructive update of R
   } else {
     break;
   }   
 }
 return R;
}

Any suggestions?

Answer 1

src/MatrixOps.chpl:605: error: halt reached - Sparse arrays can't be zippered with anything other than their domains and sibling arrays (CS layout)

I believe you are encountering a limitation of zippering sparse arrays in the current implementation, documented in #6577.

Removing some unknowns from the equation, I believe this distilled code snippet demonstrates the issue you are encountering:

use LayoutCS;

var dom = {1..10, 1..10};

var Adom: sparse subdomain(dom) dmapped CS();
var Bdom: sparse subdomain(dom) dmapped CS();

var A: [Adom] real;
var B: [Bdom] real;

Adom += (1,1);
Bdom += (1,1);

A[1,1] = 1.0;
B[1,1] = 2.0;


writeln(A.domain == B.domain); // true
var willThisWork = && reduce (A == B);
// dang.chpl:19: error: halt reached - Sparse arrays can't be zippered with
//           anything other than their domains and sibling arrays (CS layout)

As a work-around, I would suggest looping over the sparse indices after confirming the domains are equal and performing a && reduce. This is something you could wrap in a helper function, e.g.

proc main() { 
  var dom = {1..10, 1..10};

  var Adom: sparse subdomain(dom) dmapped CS();
  var Bdom: sparse subdomain(dom) dmapped CS();

  var A: [Adom] real;
  var B: [Bdom] real;

  Adom += (1,1);
  Bdom += (1,1);

  A[1,1] = 1.0;
  B[1,1] = 2.0;

  if A.domain == B.domain {
    writeln(equal(A, B));
  }
}

/* Some day, this should be A.equals(B) ! */
proc equal(A: [], B: []) {
  // You could also return 'false' if domains do not match
  assert(A.domain == B.domain);

  var s = true;

  forall (i,j) in A.domain with (&& reduce s) {
    s &&= (A[i,j] == B[i,j]);
  }

  return s;
}

src/MatrixOps.chpl:602: error: control reaches end of function that returns a value

This error is a result of not returning something in every condition. I believe you intended to do:

proc tropicLimit(A:[] real,B:[] real) {
 var R = tropic(A,B);
 if A.domain == R.domain {
   if && reduce (A == R) {
     return R;
   } else {
     return tropicLimit(R,B);
   }
 } else {
   return tropicLimit(R,B);
 }
}