New Issue: Index-neutrality of dense LinearAlgebra functions

Notifications Issues
1

17624, "ben-albrecht", "Index-neutrality of dense LinearAlgebra functions", "2021-05-03T15:19:14Z"

The LinearAlgebra library aims to be as index-neutral as possible such that
users can use any offsets or strides in their matrices and vectors (Chapel
arrays) that they wish. This issue summarizes the index-neutrality status of
all the current dense LinearAlgebra functions as of chpl version 1.25.0
pre-release (1125df1d2b).

Methodology

This survey was completed by a combination of reading the source code and
writing short tests to check cases that were not clear in the source. Ideally, we would
like to have index-neutrality tests for every LinearAlgebra function. However,
this approach was not taken for due to time-constraints. That said, we should
add tests for the cases that are fixed in this effort in order to lock in the
new behavior.

Terminology

A function is considered index-neutral if it meets the following criteria:
- Any arrays accepted as arguments can have any offset or stride in their domains
- Any arrays returned will inherit the domain from the argument arrays if possible

For conciseness, I call arrays with an offset of n as n-based arrays, e.g.
an offset of 1 (var A: [1..3] real;) is called a 1-based array.

When I say a function assumes indices should match across ranks, that means an array with mismatched indices such as var A: [1..10, 4..12] int; will cause an error. This kind of index-neutrality is not very important for most users, but we ought to support it for completeness.

Index-neutrality survey

Matrix and Vector Factory Functions

  • Vector
  • Matrix
  • eye
  • setDiag

Matrix Operations

  • transpose
  • dot
    • with BLAS
    • BLAS off
      • assumes 0-based indices (matrix-matrix multiplication)
  • matPow
    • uses `dot, so has same issues when BLAS is off
  • inner
    • assumes non-strided indices
  • outer
  • inv
  • cross
  • trace
  • lu
    • assumes 0-based non-strided indices
  • det
    • assumes 0-based non-strided indices (calls lu)
  • norm
  • solve_tril
    • assumes 0-based non-strided indices
  • solve_triu
    • assumes 0-based non-strided indices
  • solve
    • assumes 0-based non-strided indices (calls lu)
  • cholesky
    • assumes same indices across ranks
  • leastSquares
    • always returns 0-based array instead of inheriting domain
  • eigvalsh
  • eigh
    • always returns 0-based vector instead of inheriting domain
  • eigvals
  • eig
    • always returns 0-based vector instead of inheriting domain
  • svd
    • always returns 0-based array instead of inheriting domain
  • jacobi
    • assumes same indices across ranks
  • kron
    • always returns 0-based array instead of inheriting domain

Matrix Structure

  • diag
    • with a matrix argument
      • always returns 0-based vector (we should probably inherit from A.dim(0))
    • with a vector argument
  • tril
    • assumes non-strided / same indices across ranks (implementation should reindex)
  • triu
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isDiag
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isHermitian
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isSymmetric
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isTril
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isTriu
    • assumes non-strided / same indices across ranks (implementationshould reindex)
  • isSquare