Dgemm with offset - falsely says index out of bounds? #23
Description
I'm getting this error;
[122] vecxt.SubmatrixTest:
[122] ==> X vecxt.SubmatrixTest.Offset 0.08s java.lang.IndexOutOfBoundsException: Index 9 out of bounds for length 9
[122] at dev.ludovic.netlib.blas.AbstractBLAS.checkIndex(AbstractBLAS.java:51)
[122] at dev.ludovic.netlib.blas.AbstractBLAS.dgemm(AbstractBLAS.java:294)
[122] at vecxt.SubmatrixTest.$init$$$anonfun$3$$anonfun$1(submatrix.test.scala:71)
For what I believe to be a valid call to dgemm. I've minimised it as far as possible, although I'm not sure how best to present that. It's all open so easily shareable if that helps or anyone is interested.
Here's a log which tries to set out the data structures;
[122] 4.0 5.0 6.0
[122] 7.0 8.0 9.0
[122] rows: 2, cols: 3, rowStride: 1, colStride: 3, offset: 1, data length: 9
[122] submat raw data
[122] (0, 1.0),(1, 4.0),(2, 7.0),(3, 2.0),(4, 5.0),(5, 8.0),(6, 3.0),(7, 6.0),(8, 9.0)
[122] submat1 matrix
[122] 1.0
[122] 2.0
[122] 3.0
[122] rows: 3, cols: 1, rowStride: 1, colStride: 3, offset: 0, data length: 3
[122] vecxt.SubmatrixTest:
[122] ==> X vecxt.SubmatrixTest.Offset 0.08s java.lang.IndexOutOfBoundsException: Index 9 out of bounds for length 9
Here
- offset = 1
- (rows - 1) * rowStride = 1
- (cols - 1) * colStride = 6
I think the max required index, should be 8, rather than 9.
I was very surprised to find this, to the extent that I spent a lot of time checking and rechecking. I still couldn't see why this shouldn't have worked? As I'm working in scala, I ultimately ended up compiling the same test into scala native and feeding it to CBLAS. CBLAS seems to be fine with this. I'm filing this here mostly on the strength of that finding as opposed to towering confidence in my own assessment.
Here's the raw code, I can share a runnable example, if that helps.
val mat1 = Matrix.fromRows[Double](
NArray[Double](1.0, 2.0, 3.0),
NArray[Double](1.0, 2.0, 3.0) + 3.0,
NArray[Double](1.0, 2.0, 3.0) + 6.0,
)
val submat1 = Matrix.fromColumns[Double](Array(1.0, 2, 3))
for i <- 0.until(2) do
println("-----")
println(s"i: $i")
val submat = mat1.submatrix( NArray(i, i+1), :: )
println(submat.printMat)
println(submat.layout)
println("submat raw data")
println(submat.raw.zipWithIndex.map { case (v, i) => s"($i, $v)" }.mkString(","))
println("submat1 matrix")
println(submat1.printMat)
println(submat1.layout)
// this inlines to a call to dgemm**
val result = submat @@ submat1
println(result.printMat)
println("-----")first matrix m.
second matrix b
inline def `matmulInPlace!`(b: Matrix[Double], c: Matrix[Double], alpha: Double = 1.0, beta: Double = 0.0)(using
inline boundsCheck: BoundsCheck
): Unit =
dimMatCheck(m, b)
blas.dgemm(
"N",
"N",
m.rows,
b.cols,
m.cols,
alpha,
m.raw,
m.offset,
m.colStride,
b.raw,
b.offset,
b.colStride,
beta,
c.raw,
c.offset,
m.rows
)