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

  • [x] Vector
  • [x] Matrix
  • [x] eye
  • [x] setDiag

Matrix Operations

  • [x] transpose
  • dot
    • [x] 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
  • [x] outer
  • [x] inv
  • [x] cross
  • [x] trace
  • [ ] lu
    • assumes 0-based non-strided indices
  • [ ] det
    • assumes 0-based non-strided indices (calls lu)
  • [x] 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
  • [x] eigvalsh
  • [ ] eigh
    • always returns 0-based vector instead of inheriting domain
  • [x] 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))
    • [x] 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)
  • [x] isSquare