Chapel: Ensuring input given to Matrix() is accurate when given arrays as arguments

[SamuelHoward]

Summary of Issue

Note: This issue is base on a TODO in chapel/modules/packages/LinearAlgebra.chpl on line 383

The current Matrix() function in the LinearAlgebra package is very versatile. One of the most intuitive uses of Matrix() is when giving it multiple 1-dimensional arrays of equal size and a type:

var matrix1 = Matrix([1,2,3], [4,5,6], [7,8,9], eltType=int);
// Returns matrix:
// 1 2 3
// 4 5 6
// 7 8 9

Naturally, when given arrays that are not uniform, Matrix() fails:

var matrix1 = Matrix([1,2,3,4], [4,5,6], [7,8,9], eltType=int);
// Gives the error:
// error: halt reached - array index out of bounds

It would be worthwhile to add a simple check in Matrix() such that a more precise error is given, such as:

var matrix1 = Matrix([1,2,3,4], [4,5,6], [7,8,9], eltType=int);
// Gives the error:
// error: halt reached - Matrix() expected arrays of equal length

Proposed Solution

The current implementation of Matrix() when given arrays as arguments, in chapel/modules/packages/LinearAlgebra.chpl, is as follows:

proc Matrix(const Arrays: ?t ...?n, type eltType) where isArrayType(t) && t.rank == 1 {
  // TODO -- assert all array domains are same length
  //         Can this be done via type query?

  if Arrays(1).domain.rank != 1 then compilerError("Matrix() expected 1D arrays");

  const dim2 = 1..Arrays(1).domain.dim(1).size,
        dim1 = 1..n;

  var M: [{dim1, dim2}] eltType;

  for i in dim1 do
    M[i, ..] = Arrays(i)[..]: eltType;

  return M;
}

Adding a forall loop would resolve this issue, albeit not optimally. As the TODO above says, using type query would be preferred. Below is a simple solution using a forall loop:

proc Matrix(const Arrays: ?t ...?n, type eltType) where isArrayType(t) /* && t.rank == 1 */ {

    if Arrays(1).domain.rank != 1 then compilerError("Matrix() expected 1D arrays");

    // Added forall loop:
    forall i in 2..n do
      if Arrays(i).size != Arrays(1).size then halt("Matrix() expected arrays of equal length");

    const dim2 = 1..Arrays(1).domain.dim(1).size,
    dim1 = 1..n;

    var M: [{dim1, dim2}] eltType;

    for i in dim1 do
      M[i, ..] = Arrays(i)[..]: eltType;

    return M;
  }

Final Note: This is my first attempt at creating an issue for Chapel. I tried to follow the guidelines, but please let me know if I made a misstep anywhere. Thank you.

[bradcray]

This looks like a job for @ben-albrecht (already done for today, but he should be back tomorrow).

[ben-albrecht]

@SamuelHoward - Thanks for tackling this. A few thoughts:

As the TODO above says, using type query would be preferred.

As the comment implies, I don't think this is possible in Chapel today, so your solution is reasonable.

Adding a forall loop would resolve this issue, albeit not optimally.

Note that there is some overhead when creating tasks for a forall loop, such that it's generally not worth creating the forall unless the inner operation is expensive and/or the array is large (e.g. >1MB). To amortize the cost of the forall creation, we could also utilize the same loop for the assignments below:

    for i in dim1 do
      M[i, ..] = Arrays(i)[..]: eltType;

e.g.

    forall i in dim1 {
      if i > 1 {  /* size check */ }
      M[i, ..] = Arrays(i)[..]: eltType;
    }

We could return to this in the future and try to optimize the matrix creation time further by adding some heuristically determined checks on the matrix size, but let's keep it simple for now.

Another optimization we may include is to only check these sizes if bounds checking is enabled (optimizations disabled), e.g.

    forall i in dim1 {
      boundsChecking {
        if i > 1 {  /* size check */ }
      }
      M[i, ..] = Arrays(i)[..]: eltType;
    }

Would you like to open a PR with these changes and a test demonstrating the error?

[SamuelHoward]

@ben-albrecht - Thank you for your comprehensive response.

What you said makes sense. I will make a corresponding PR shortly. I will comment here with a link to the PR.

Note: I submitted a CLA so that I may contribute. It hasn't been processed quite yet though.

[SamuelHoward]

Here is the PR #15322

Also, my CLA has fortunately been processed as of now.

[ben-albrecht]

Closing now that #15322 is merged. Thanks @SamuelHoward !