Chapel: Cannot multiply sparse matrix with derived vector

[buddha314]

Almost certainly related to https://github.com/chapel-lang/chapel/issues/8023 but I don't want presume.

Once again our story begins with a matrix A. It is sparse. I want to multiply this matrix with a 1 vector.

use LinearAlgebra, LayoutCS;

var nv: int = 8,
    D: domain(2) = {1..nv, 1..nv},
    SD: sparse subdomain(D),
    A: [SD] real;

SD += (1,2); A[1,2] = 1;
SD += (1,3); A[1,3] = 1;
SD += (1,4); A[1,4] = 1;
SD += (2,2); A[2,2] = 1;
SD += (2,4); A[2,4] = 1;
SD += (3,4); A[3,4] = 1;
SD += (4,5); A[4,5] = 1;
SD += (5,6); A[5,6] = 1;
SD += (6,7); A[6,7] = 1;
SD += (6,8); A[6,8] = 1;
SD += (7,8); A[7,8] = 1;

writeln(A);
writeln(diag(A));
writeln(diag(diag(A)));

Good so far. But wait, there is a villain in this story as well.

> chpl -I/usr/local/Cellar/openblas/0.2.20/include -L/usr/local/Cellar/openblas/0.2.20/lib -lblas digs.chpl
digs.chpl:26: error: unresolved call '[DefaultSparseDom(2,int(64),domain(2,int(64),false))] real(64).dot([domain(2,int(64),false)] real(64))'
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:435: note: candidates are: dot(A: [?Adom] ?eltType, B: [?Bdom] eltType)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:445: note:                 _array.dot(A: [] )
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:452: note:                 dot(A: [?Adom] ?eltType, b)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:460: note:                 dot(a, B: [] )
$CHPL_HOME/modules/packages/BLAS.chpl:2074: note:                 dot(X: [?xD] ?eltType, Y: [?yD] eltType, incY: c_int = 1, incX: c_int = 1)

I cannot multiply this matrix with it's own diagonal. Why would I want to? Because I'm crazy that way. Next, I can't multiply it with a 1 vector.

 chpl -I/usr/local/Cellar/openblas/0.2.20/include -L/usr/local/Cellar/openblas/0.2.20/lib -lblas digs.chpl
digs.chpl:26: error: unresolved call '[DefaultSparseDom(2,int(64),domain(2,int(64),false))] real(64).dot([domain(2,int(64),false)] real(64))'
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:435: note: candidates are: dot(A: [?Adom] ?eltType, B: [?Bdom] eltType)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:445: note:                 _array.dot(A: [] )
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:452: note:                 dot(A: [?Adom] ?eltType, b)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:460: note:                 dot(a, B: [] )
$CHPL_HOME/modules/packages/BLAS.chpl:2074: note:                 dot(X: [?xD] ?eltType, Y: [?yD] eltType, incY: c_int = 1, incX: c_int = 1)

What is the proper way to do this?

[ben-albrecht]

I see 2 problems with your example:

1. LinearAlgebra.Sparse only works with CSR layout.

var SD: sparse subdomain(D) corresponds to COO layout.

You will need var SD: sparse subdomain(D) dmapped CS(compressRows=true) to specify CSR layout.

1. LinearAlgebra.Sparse only supports a subset of LinearAlgebra procedures. diag() is unfortunately not yet supported as part of that subset.

I'll interpret this issue as a feature request for LinearAlgebra.Sparse.diag().

See LinearAlgebra.Sparse docs for more info.

[buddha314]

Thanks for the update, can you think of a way to fake diag() on a sparse matrix? If I created a new domain as a sparse sub-domain of D and added all the (i,i)'s would I be able to multiply them?

[ben-albrecht]

If I created a new domain as a sparse sub-domain of D and added all the (i,i)'s would I be able to multiply them?

Yes, this should work, assuming you're using CSR layout.

[buddha314]

Building on the example above, I tried this:

var diagDom: sparse subdomain(D)  dmapped CS(compressRows=true);
for i in 1..nv {
  diagDom += (i,i);
}
writeln(A[diagDom]);

And it spat venom during compilation:

 chpl -I/usr/local/Cellar/openblas/0.2.20/include -L/usr/local/Cellar/openblas/0.2.20/lib -lblas digs.chpl
digs.chpl:31: error: unresolved call 'CSDom(2,int(64),domain(2,int(64),false),true,false).dsiGetIndices()'
$CHPL_HOME/modules/internal/DefaultRectangular.chpl:115: note: candidates are: DefaultRectangularDom.dsiGetIndices()
$CHPL_HOME/modules/internal/DefaultOpaque.chpl:95: note:                 DefaultOpaqueDom.dsiGetIndices()

Although, without the writeln it worked okay.

[ben-albrecht]

The venom is from trying to slice a sparse array by a sparse subdomain. That's not supported (at least today).

This is a little more verbose, but gets the job done:

var B: [diagDom] int;

forall i in sD.dim(1) {
  if SD.member((i,i)) {
    B[i, i] = A[i, i];
  }
}

writeln(B);

[buddha314]

Okay, now I can get B, but then....

writeln("\nA.dot(B)\n", A.dot(B));

on compilation this violently expectorates

chpl -I/usr/local/Cellar/openblas/0.2.20/include -L/usr/local/Cellar/openblas/0.2.20/lib -lblas digs.chpl
digs.chpl:41: error: unresolved call '[CSDom(2,int(64),domain(2,int(64),false),true,false)] real(64).dot([CSDom(2,int(64),domain(2,int(64),false),true,false)] real(64))'
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:435: note: candidates are: dot(A: [?Adom] ?eltType, B: [?Bdom] eltType)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:445: note:                 _array.dot(A: [] )
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:452: note:                 dot(A: [?Adom] ?eltType, b)
$CHPL_HOME/modules/packages/LinearAlgebra.chpl:460: note:                 dot(a, B: [] )
$CHPL_HOME/modules/packages/BLAS.chpl:2074: note:                 dot(X: [?xD] ?eltType, Y: [?yD] eltType, incY: c_int = 1, incX: c_int = 1)

[ben-albrecht]

I can't determine the origin of your expectoration without seeing the full example.

Here is a working example derived from your original question.

[buddha314]

Right, that works for me, too. It has a few patterns I haven't been using, so I'll have to work those in. For instance, I've been calling dmapped CS() instead of CSRDomain(). Yours is much cleaner!

[buddha314]

Also, I don't understand how this bit works

const nv = 8;
var SD = CSRDomain(1..nv),
     A = CSRMatrix(SD);

1..nv is one dimensional, so isn't SD one dimensional as well?

[ben-albrecht]

1..nv is one dimensional, so isn't SD one dimensional as well?

See the third overload in the docs.

CSR domains can only be 2D, so the range is interpreted as the range for both dimensions. You can also give it a 2D domain if that makes you feel more comfortable:

var SD = CSRDomain({1..nv, 1..nv});

[buddha314]

All right! This is much cleaner!

What if I don't know the number of rows/columns to begin with?

For instance, I have a Graph class in Chingon that wants a sparse matrix as a field. Currently I'm doing

  class Graph {
    var vdom: domain(2),
        SD: sparse subdomain(vdom) dmapped  CS(compressRows=true),
        A: [SD] real,
        name: string,
        vnames: domain(string),
        vids: [vnames] int,
        nameIndex: [1..0] string;

    proc init(A: []) {
      this.vdom = {A.domain.dim(1), A.domain.dim(2)};
      super.init();
      for ij in A.domain {
        this.SD += ij;
        this.A(ij) = A(ij);
      }
    }

But I'm not sure how to use the CSR* stuff. This seems to throw an error.

  class Graph {
    var vdom: domain(1),
        SD = CSRDomain(1..0),
        A = CSRMatrix(SD),
        name: string,
        vnames: domain(string),
        vids: [vnames] int,
        nameIndex: [1..0] string;

    proc init(A: []) {
      this.vdom = {A.domain.dim(1), A.domain.dim(2)};  // Yes, this will error, I'll sort it
      super.init();
      for ij in A.domain {
        this.SD += ij;
        this.A(ij) = A(ij);
      }
    }

So I don't know how big A will be until the object is constructed.

[ben-albrecht]

This issue seems to boil down to a feature request for supporting dynamically bounded sparse arrays, which is the focus of issue #7544. I'm going to close this issue as duplicate of that one.

Let me know if you disagree with that observation.