diff --git a/vecxt/src-js/floatmatrix.scala b/vecxt/src-js/floatmatrix.scala index 7ab25cf7..9d4010f2 100644 --- a/vecxt/src-js/floatmatrix.scala +++ b/vecxt/src-js/floatmatrix.scala @@ -11,6 +11,107 @@ object JsFloatMatrix: extension (m: Matrix[Float]) + @targetName("matmulFloat") + def @@(b: Matrix[Float]): Matrix[Float] = + m.matmul(b, 1.0f, 0.0f) + + @targetName("matmulFloatNonDefault") + def matmul(b: Matrix[Float], alpha: Float, beta: Float): Matrix[Float] = + dimMatCheck(m, b) + val newArr: Array[Float] = Array.ofDim[Float](m.rows * b.cols) + val newmat = Matrix[Float](newArr, m.rows, b.cols) + m.`matmulInPlace!`(b, newmat, alpha, beta) + newmat + end matmul + + /** Writes `alpha * (m @@ b) + beta * c` into `c` in place, via a JS `sgemm` shim. Float counterpart of + * `JsDoubleMatrix.matmulInPlace!`; see there for the reasoning behind the layout choices, which carries over + * unchanged aside from element type. + * + * `c` must already be shaped `(m.rows, b.cols)` and dense column-major — `ldc` is hardcoded to `m.rows` below, and + * `sgemm` also reads `c` when `beta != 0`, so any other shape or layout would be silently written to (or read + * from) incorrectly rather than rejected. Use `matmul`/`@@` instead if you don't already have a conforming `c` to + * write into; they allocate one for you. + */ + @targetName("matmulFloatInPlace") + def `matmulInPlace!`(b: Matrix[Float], c: Matrix[Float], alpha: Float, beta: Float): Unit = + dimMatCheck(m, b) + matmulOutputCheck(m, b, c) + println("PERFORMING WARNING in matmul on JS") + println("THIS method copies into native JS types. Then copies back out. Expect catastrophic performance.") + + if m.hasSimpleContiguousMemoryLayout && b.hasSimpleContiguousMemoryLayout then + val lda = if m.isDenseColMajor then m.rows else m.cols + val ldb = if b.isDenseColMajor then b.rows else b.cols + + val transB = if b.isDenseColMajor then "no-transpose" else "transpose" + val transA = if m.isDenseColMajor then "no-transpose" else "transpose" + + // See JsDoubleMatrix.matmulInPlace! for why `order` is always "column-major" here regardless of m/b's own + // orientation. + val outArr = new Float32Array(c.raw.toJSArray) + sgemm( + "column-major", + transA, + transB, + m.rows, + b.cols, + m.cols, + alpha, + new Float32Array(m.raw.toJSArray), + lda, + new Float32Array(b.raw.toJSArray), + ldb, + beta, + outArr, + m.rows + ) + // copy result back into c.raw (Scala Array[Float]) element-wise + val copyLen = Math.min(outArr.length, c.raw.length) + var ci = 0 + while ci < copyLen do + c.raw(ci) = outArr(ci) + ci += 1 + end while + else if blasLeadingDimensionCheck(m) && blasLeadingDimensionCheck(b) then + val transB = if b.rowStride == 1 then "no-transpose" else "transpose" + val transA = if m.rowStride == 1 then "no-transpose" else "transpose" + + // See the fully-dense branch above: `order` stays "column-major" regardless of m/b's own orientation, + // since transA/transB/lda/ldb already implement the transpose trick for that fixed order. + val outArr = new Float32Array(c.raw.toJSArray) + sgemm( + "column-major", + transA, + transB, + m.rows, + b.cols, + m.cols, + alpha, + // convert backing Scala Array[Float] to Float32Array slice (copies) + new Float32Array(m.raw.toJSArray).subarray(m.offset), + if m.rowStride == 1 then m.colStride else m.rowStride, + new Float32Array(b.raw.toJSArray).subarray(b.offset), + if b.rowStride == 1 then b.colStride else b.rowStride, + beta, + outArr, + m.rows + ) + // copy result back into c.raw (Scala Array[Float]) element-wise + val copyLen2 = Math.min(outArr.length, c.raw.length) + var cj = 0 + while cj < copyLen2 do + c.raw(cj) = outArr(cj) + cj += 1 + end while + else + throw UnsupportedLayoutException( + s"matmulInPlace! does not support this combination of matrix layouts. m: ${m.layoutString}, b: ${b.layoutString}" + ) + end if + + end `matmulInPlace!` + /** Writes `alpha * (m @@ vec) + beta * y` into `y` in place, via the stdlib `sgemv` shim. JS counterpart of * `JvmFloatMatrix.*=`; see `JvmDoubleMatrix.*=` for the reasoning behind the guards. * diff --git a/vecxt/src-js/stdlib.facade.scala b/vecxt/src-js/stdlib.facade.scala index ccfb79af..b7e86416 100644 --- a/vecxt/src-js/stdlib.facade.scala +++ b/vecxt/src-js/stdlib.facade.scala @@ -56,6 +56,31 @@ object dgemm extends js.Object: end dgemm +/** `sgemm( ord, ta, tb, M, N, K, α, A, lda, B, ldb, β, C, ldc )` — `C = α*op(A)*op(B) + β*C`. Single-precision twin of + * [[dgemm]]; see its scaladoc, which carries over unchanged aside from element type. + */ +@js.native +@JSImport("@stdlib/blas/base/sgemm/lib", JSImport.Default) +object sgemm extends js.Object: + def apply( + ord: String, + transA: String, + transB: String, + m: Int, + n: Int, + k: Int, + alpha: Float, + a: Float32Array, + lda: Int, + b: Float32Array, + ldb: Int, + beta: Float, + c: Float32Array, + ldc: Int + ): Unit = js.native + +end sgemm + @js.native @JSImport("@stdlib/blas/base/dgemv/lib", JSImport.Default) object dgemv extends js.Object: diff --git a/vecxt/src-native/floatmatrix.scala b/vecxt/src-native/floatmatrix.scala index 39aa09cc..bdc878b6 100644 --- a/vecxt/src-native/floatmatrix.scala +++ b/vecxt/src-native/floatmatrix.scala @@ -12,6 +12,84 @@ object NativeFloatMatrix: extension (m: Matrix[Float]) + @targetName("matmulFloat") + def @@(b: Matrix[Float]): Matrix[Float] = + m.matmul(b, 1.0f, 0.0f) + + @targetName("matmulFloatNonDefault") + def matmul(b: Matrix[Float], alpha: Float, beta: Float): Matrix[Float] = + dimMatCheck(m, b) + val newArr: Array[Float] = Array.ofDim[Float](m.rows * b.cols) + val newmat = Matrix[Float](newArr, m.rows, b.cols) + m.`matmulInPlace!`(b, newmat, alpha, beta) + newmat + end matmul + + /** Writes `alpha * (m @@ b) + beta * c` into `c` in place, via `cblas_sgemm`. Float counterpart of + * `NativeDoubleMatrix.matmulInPlace!`; see there for the reasoning behind the `trans`/`order` choices, which + * carries over unchanged aside from element type. + * + * `c` must already be shaped `(m.rows, b.cols)` and dense column-major — `ldc` is hardcoded to `m.rows` below, and + * `sgemm` also reads `c` when `beta != 0`, so any other shape or layout would be silently written to (or read + * from) incorrectly rather than rejected. Use `matmul`/`@@` instead if you don't already have a conforming `c` to + * write into; they allocate one for you. + */ + @targetName("matmulFloatInPlace") + def `matmulInPlace!`(b: Matrix[Float], c: Matrix[Float], alpha: Float, beta: Float): Unit = + dimMatCheck(m, b) + matmulOutputCheck(m, b, c) + + if m.hasSimpleContiguousMemoryLayout && b.hasSimpleContiguousMemoryLayout then + val lda = if m.isDenseColMajor then m.rows else m.cols + val ldb = if b.isDenseColMajor then b.rows else b.cols + val transB = if b.isDenseColMajor then blasEnums.CblasNoTrans else blasEnums.CblasTrans + val transA = if m.isDenseColMajor then blasEnums.CblasNoTrans else blasEnums.CblasTrans + + // See NativeDoubleMatrix.matmulInPlace! for why `order` is always CblasColMajor here. + blas.cblas_sgemm( + blasEnums.CblasColMajor, + transA, + transB, + m.rows, + b.cols, + m.cols, + alpha, + m.raw.at(0), + lda, + b.raw.at(0), + ldb, + beta, + c.raw.at(0), + m.rows + ) + else if blasLeadingDimensionCheck(m) && blasLeadingDimensionCheck(b) then + val transB = if b.rowStride == 1 then blasEnums.CblasNoTrans else blasEnums.CblasTrans + val transA = if m.rowStride == 1 then blasEnums.CblasNoTrans else blasEnums.CblasTrans + // See the fully-dense branch above: order stays CblasColMajor regardless of m/b's own orientation. + blas.cblas_sgemm( + blasEnums.CblasColMajor, + transA, + transB, + m.rows, + b.cols, + m.cols, + alpha, + m.raw.at(m.offset), + if m.rowStride == 1 then m.colStride else m.rowStride, + b.raw.at(b.offset), + if b.rowStride == 1 then b.colStride else b.rowStride, + beta, + c.raw.at(c.offset), + m.rows + ) + else + throw UnsupportedLayoutException( + s"matmulInPlace! does not support this combination of matrix layouts. m: ${m.layoutString}, b: ${b.layoutString}" + ) + + end if + end `matmulInPlace!` + /** Writes `alpha * (m @@ vec) + beta * y` into `y` in place, via CBLAS `cblas_sgemv`. Native counterpart of * `JvmFloatMatrix.*=`; see `JvmDoubleMatrix.*=` for the reasoning behind the guards. * diff --git a/vecxt/test/src-jvm/TODO.test.scala b/vecxt/test/src-jvm/TODO.test.scala deleted file mode 100644 index 849f1c4e..00000000 --- a/vecxt/test/src-jvm/TODO.test.scala +++ /dev/null @@ -1,32 +0,0 @@ -package vecxt - -import munit.FunSuite - -import all.* - -class TODO extends FunSuite: - - test("matmulInPlace! throws for unsupported general Float layouts"): - val left = Matrix[Float]( - Array[Float](1.0f, 90.0f, 2.0f, 91.0f, 92.0f, 3.0f, 93.0f, 4.0f), - 2, - 2, - 2, - 5, - 0 - ) - val right = Matrix[Float]( - Array[Float](5.0f, 80.0f, 6.0f, 81.0f, 82.0f, 7.0f, 83.0f, 8.0f), - 2, - 2, - 2, - 5, - 0 - ) - val out = Matrix.zeros[Float]((2, 2)) - - intercept[UnsupportedLayoutException] { - left.`matmulInPlace!`(right, out, alpha = 1.0f, beta = 0.0f) - } - -end TODO diff --git a/vecxt/test/src/floatMatMul.test.scala b/vecxt/test/src/floatMatMul.test.scala new file mode 100644 index 00000000..9acadb86 --- /dev/null +++ b/vecxt/test/src/floatMatMul.test.scala @@ -0,0 +1,217 @@ +package vecxt + +import all.* +import munit.FunSuite + +/** `Matrix[Float] @@ Matrix[Float]`, `matmul`, and `matmulInPlace!` across every memory layout. Shared rather than + * JVM-only because the operation now exists on every platform — netlib `sgemm` on the JVM, `cblas_sgemm` on Native, + * and the stdlib `sgemm` shim on JS — following the same layout rules as [[DifferentMemoryLayoutTests]] does for + * `Double`. Previously `matmulInPlace!`, `matmul`, and `@@` existed only on the JVM for `Float`. + */ +class FloatMatMulSuite extends FunSuite: + + test("scalars in Float matmul") { + def makeMat = Matrix[Float](Array.tabulate[Float](9)(_.toFloat + 1), 3, 3, 3, 1, 0) + val eye = Matrix.eye[Float](3) + + assertMatrixEquals(makeMat @@ eye, makeMat) + + val doubledMat = Matrix[Float](Array.tabulate[Float](9)(idx => (idx.toFloat + 1) * 2.0f), 3, 3, 3, 1, 0) + assertMatrixEquals(eye.matmul(makeMat, 2.0f, 0.0f), doubledMat) + + val outMat = Matrix.eye[Float](3) + + makeMat.`matmulInPlace!`(eye, outMat, 2.0f, 2.0f) + + // hand-computed: alpha * (makeMat @@ eye) + beta * outBefore, outBefore == eye + val expected = Matrix[Float]( + Array.tabulate[Float](9) { idx => + val row = idx % 3 + val col = idx / 3 + doubledMat(row, col) + (if row == col then 2.0f else 0.0f) + }, + 3, + 3, + 1, + 3, + 0 + ) + assertMatrixEquals(outMat, expected) + } + + test("Float matmul col major * row major, all combinations") { + val matRow = Matrix[Float](Array.tabulate[Float](9)(_.toFloat + 1), 3, 3, 3, 1, 0) + val matCol = Matrix[Float](Array.tabulate[Float](9)(_.toFloat + 1), 3, 3, 1, 3, 0) + + val mat = matCol @@ matCol + assertEqualsDouble(mat(0, 0).toDouble, 1 * 1 + 4 * 2 + 3 * 7, 0.0001) + assertEqualsDouble(mat(0, 2).toDouble, 1 * 7 + 4 * 8 + 7 * 9, 0.0001) + assertEqualsDouble(mat(1, 1).toDouble, 4 * 2 + 5 * 5 + 8 * 6, 0.0001) + assertEqualsDouble(mat(2, 0).toDouble, 3 * 1 + 6 * 2 + 9 * 3, 0.0001) + + val mat2 = matCol @@ matRow + assertEqualsDouble(mat2(0, 0).toDouble, 1 * 1 + 4 * 4 + 7 * 7, 0.0001) + assertEqualsDouble(mat2(0, 2).toDouble, 1 * 3 + 4 * 6 + 7 * 9, 0.0001) + assertEqualsDouble(mat2(1, 1).toDouble, 2 * 2 + 5 * 5 + 8 * 8, 0.0001) + assertEqualsDouble(mat2(2, 0).toDouble, 3 * 1 + 6 * 4 + 9 * 7, 0.0001) + + val mat3 = matRow @@ matRow + assertEqualsDouble(mat3(0, 0).toDouble, 1 * 1 + 2 * 4 + 7 * 3, 0.0001) + assertEqualsDouble(mat3(0, 2).toDouble, 1 * 3 + 2 * 6 + 3 * 9, 0.0001) + assertEqualsDouble(mat3(1, 1).toDouble, 4 * 2 + 5 * 5 + 6 * 8, 0.0001) + assertEqualsDouble(mat3(2, 0).toDouble, 7 * 1 + 8 * 4 + 9 * 7, 0.0001) + + val mat4 = matRow @@ matCol + assertEqualsDouble(mat4(0, 0).toDouble, 1 * 1 + 2 * 2 + 3 * 3, 0.0001) + assertEqualsDouble(mat4(0, 2).toDouble, 1 * 7 + 2 * 8 + 3 * 9, 0.0001) + assertEqualsDouble(mat4(1, 1).toDouble, 4 * 4 + 5 * 5 + 6 * 6, 0.0001) + assertEqualsDouble(mat4(2, 0).toDouble, 7 * 1 + 8 * 2 + 9 * 3, 0.0001) + } + + test("Float matmul with offset (submatrix) views") { + val mat1 = Matrix.fromRows( + Array(1.0f, 2.0f, 3.0f, 4.0f), + Array(5.0f, 6.0f, 7.0f, 8.0f), + Array(9.0f, 10.0f, 11.0f, 12.0f), + Array(13.0f, 14.0f, 15.0f, 16.0f) + ) + val mat2 = Matrix.fromRows( + Array(1.0f, 2.0f, 3.0f, 4.0f), + Array(5.0f, 6.0f, 7.0f, 8.0f), + Array(9.0f, 10.0f, 11.0f, 12.0f), + Array(13.0f, 14.0f, 15.0f, 16.0f), + Array(1.0f, 2.0f, 3.0f, 4.0f) + ) + + val subMat = Range.Inclusive(1, 2, 1) + + val zeroCopy = mat1(subMat, subMat) + val zeroCopy2 = mat2(subMat, subMat) + + val newMat = zeroCopy @@ zeroCopy2 + + assertEqualsDouble(newMat(0, 0).toDouble, 6 * 6 + 7 * 10, 0.000001) + assertEqualsDouble(newMat(1, 0).toDouble, 10 * 6 + 11 * 10, 0.000001) + assertEqualsDouble(newMat(1, 1).toDouble, 10 * 7 + 11 * 11, 0.000001) + assertEqualsDouble(newMat(0, 1).toDouble, 6 * 7 + 7 * 11, 0.000001) + + val view1 = mat1(Range.Inclusive(0, 2, 1), Range.Inclusive(1, 2, 1)) + val view2 = mat2(Range.Inclusive(1, 2, 1), Range.Inclusive(0, 3, 1)) + val viewMul = view1 @@ view2 + + assertMatrixEquals( + viewMul, + Matrix.fromRows( + Array(37.0f, 42.0f, 47.0f, 52.0f), + Array(93.0f, 106.0f, 119.0f, 132.0f), + Array(149.0f, 170.0f, 191.0f, 212.0f) + ) + ) + } + + test("Float matmul non-square, all layout combinations") { + // A (3x2), logical: [[1,2],[3,4],[5,6]] + val aColMajor = Matrix[Float](Array(1.0f, 3.0f, 5.0f, 2.0f, 4.0f, 6.0f), 3, 2, 1, 3, 0) + val aRowMajor = Matrix[Float](Array(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f), 3, 2, 2, 1, 0) + + // B (2x4), logical: [[1,2,3,4],[5,6,7,8]] + val bColMajor = Matrix[Float](Array(1.0f, 5.0f, 2.0f, 6.0f, 3.0f, 7.0f, 4.0f, 8.0f), 2, 4, 1, 2, 0) + val bRowMajor = Matrix[Float](Array(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f), 2, 4, 4, 1, 0) + + val expected = Matrix.fromRows( + Array(11.0f, 14.0f, 17.0f, 20.0f), + Array(23.0f, 30.0f, 37.0f, 44.0f), + Array(35.0f, 46.0f, 57.0f, 68.0f) + ) + + assertMatrixEquals(aColMajor @@ bColMajor, expected) + assertMatrixEquals(aColMajor @@ bRowMajor, expected) + assertMatrixEquals(aRowMajor @@ bColMajor, expected) + assertMatrixEquals(aRowMajor @@ bRowMajor, expected) + } + + test("Float matmul non-square, beta accumulation") { + val a = Matrix.fromRows(Array(1.0f, 2.0f), Array(3.0f, 4.0f), Array(5.0f, 6.0f)) + val b = Matrix.fromRows(Array(1.0f, 2.0f, 3.0f, 4.0f), Array(5.0f, 6.0f, 7.0f, 8.0f)) + val out = Matrix.fromRows( + Array(1.0f, 1.0f, 1.0f, 1.0f), + Array(1.0f, 1.0f, 1.0f, 1.0f), + Array(1.0f, 1.0f, 1.0f, 1.0f) + ) + + a.`matmulInPlace!`(b, out, 2.0f, 3.0f) + + // hand-computed: alpha * (a @@ b) + beta * outBefore + assertMatrixEquals( + out, + Matrix.fromRows( + Array(25.0f, 31.0f, 37.0f, 43.0f), + Array(49.0f, 63.0f, 77.0f, 91.0f), + Array(73.0f, 95.0f, 117.0f, 139.0f) + ) + ) + } + + test("Float matmul rejects doubly-strided b with unit-rowStride m") { + val m = Matrix[Float](Array.tabulate[Float](6)(_.toFloat + 1), 2, 2, 1, 3, 0) + val b = Matrix[Float](Array.tabulate[Float](10)(_.toFloat + 1), 2, 2, 2, 5, 0) + + intercept[UnsupportedLayoutException] { + m @@ b + } + } + + test("Float matmul rejects a broadcast operand rather than handing BLAS lda = 0") { + val broadcast = Matrix[Float](Array(1.0f, 2.0f), 2, 2, 1, 0, 0) + val dense = Matrix.fromRows(Array(1.0f, 0.0f), Array(0.0f, 1.0f)) + + assert(broadcast.rowStride == 1 && broadcast.colStride == 0, "fixture must be the half-satisfying case") + intercept[UnsupportedLayoutException](broadcast @@ dense) + intercept[UnsupportedLayoutException](dense @@ broadcast) + } + + test("Float matmulInPlace! throws when c is row-major instead of column-major") { + val a = Matrix.fromRows(Array(1.0f, 2.0f, 3.0f), Array(4.0f, 5.0f, 6.0f)) // 2x3 + val b = Matrix.fromRows(Array(1.0f, 2.0f), Array(3.0f, 4.0f), Array(5.0f, 6.0f)) // 3x2 + val cRowMajor = Matrix[Float](Array.ofDim[Float](4), 2, 2, 2, 1, 0) // correctly shaped, but row-major + + intercept[UnsupportedLayoutException] { + a.`matmulInPlace!`(b, cRowMajor, 1.0f, 0.0f) + } + } + + test("Float matmulInPlace! throws MatrixDimensionMismatch when c is the wrong shape") { + val a = Matrix.fromRows(Array(1.0f, 2.0f, 3.0f), Array(4.0f, 5.0f, 6.0f)) // 2x3 + val b = Matrix.fromRows(Array(1.0f, 2.0f), Array(3.0f, 4.0f), Array(5.0f, 6.0f)) // 3x2 + val wrongSizeC = Matrix.zeros[Float]((3, 3)) // should be (2, 2) == (a.rows, b.cols) + + intercept[MatrixDimensionMismatch] { + a.`matmulInPlace!`(b, wrongSizeC, 1.0f, 0.0f) + } + } + + test("Float matmulInPlace! throws for unsupported general layouts") { + val left = Matrix[Float]( + Array[Float](1.0f, 90.0f, 2.0f, 91.0f, 92.0f, 3.0f, 93.0f, 4.0f), + 2, + 2, + 2, + 5, + 0 + ) + val right = Matrix[Float]( + Array[Float](5.0f, 80.0f, 6.0f, 81.0f, 82.0f, 7.0f, 83.0f, 8.0f), + 2, + 2, + 2, + 5, + 0 + ) + val out = Matrix.zeros[Float]((2, 2)) + + intercept[UnsupportedLayoutException] { + left.`matmulInPlace!`(right, out, alpha = 1.0f, beta = 0.0f) + } + } + +end FloatMatMulSuite