Merge pull request #37 from Quafadas/copilot/fix-7ebc2283-e18a-45e5-b01d-9656d56639cf

Implement Hadamard product for matrices with different memory layouts

Author: Quafadas

vecxt/src-js/array_native.scala

@@ -456,6 +456,17 @@ object JsNativeDoubleArrays:
       res
     end *
 
+    inline def *=(d: NArray[Double])(using inline boundsCheck: BoundsCheck): Unit =
+      dimCheck(vec, d)
+      val n = vec.length
+
+      var i = 0
+      while i < n do
+        vec(i) = vec(i) * d(i)
+        i += 1
+      end while
+    end *=
+
     inline def outer(other: NArray[Double])(using ClassTag[Double]): Matrix[Double] =
       val n = vec.length
       val m = other.length

vecxt/src-jvm/arrays.scala

@@ -1176,6 +1176,23 @@ object arrays:
       out
     end *
 
+    inline def *=(d: Array[Double])(using inline boundsCheck: BoundsCheck): Unit =
+      dimCheck(vec, d)
+      var i = 0
+      while i < spd.loopBound(vec.length) do
+        DoubleVector
+          .fromArray(spd, vec, i)
+          .mul(DoubleVector.fromArray(spd, d, i))
+          .intoArray(vec, i)
+        i += spdl
+      end while
+
+      while i < vec.length do
+        vec(i) = vec(i) * d(i)
+        i = i + 1
+      end while
+    end *=
+
     inline def /(d: Array[Double])(using inline boundsCheck: BoundsCheck): Array[Double] =
       dimCheck(vec, d)
       val out = new Array[Double](vec.length)

vecxt/src-native/array_js_native.scala

@@ -456,6 +456,17 @@ object JsNativeDoubleArrays:
       res
     end *
 
+    inline def *=(d: NArray[Double])(using inline boundsCheck: BoundsCheck): Unit =
+      dimCheck(vec, d)
+      val n = vec.length
+
+      var i = 0
+      while i < n do
+        vec(i) = vec(i) * d(i)
+        i += 1
+      end while
+    end *=
+
     inline def outer(other: NArray[Double])(using ClassTag[Double]): Matrix[Double] =
       val n = vec.length
       val m = other.length

vecxt/src/MatrixInstance.scala

@@ -127,10 +127,23 @@ object MatrixInstance:
       * @return
       */
     def deepCopy(using ct: ClassTag[A]): Matrix[A] =
+      deepCopy(asRowMajor = false)
+    end deepCopy
+
+    /** Returns a deep copy of the matrix with specified layout. Copies elements one by one.
+      *
+      * @param asRowMajor
+      *   If true, returns row-major layout; if false, returns column-major layout
+      * @param ct
+      * @return
+      */
+    def deepCopy(asRowMajor: Boolean)(using ct: ClassTag[A]): Matrix[A] =
       // println(s"Deep copying matrix with shape ${m.shape} and offset ${m.offset}")
       import BoundsCheck.DoBoundsCheck.no
       val newRaw = NArray.ofSize[A](m.numel)
-      val newMat = Matrix(newRaw, m.rows, m.cols, 1, m.rows, 0)
+      val newMat =
+        if asRowMajor then Matrix(newRaw, m.rows, m.cols, m.cols, 1, 0) // row-major: rowStride = cols, colStride = 1
+        else Matrix(newRaw, m.rows, m.cols, 1, m.rows, 0) // column-major: rowStride = 1, colStride = rows
       var i = 0
       for row <- 0 until m.rows do
         for col <- 0 until m.cols do

vecxt/src/doublematrix.scala

@@ -142,9 +142,38 @@ object DoubleMatrix:
       sameDimMatCheck(m, m2)
 
       if sameDenseElementWiseMemoryLayoutCheck(m, m2) then
+        // Fast path: use SIMD-optimized array multiplication
         val newArr = vecxt.arrays.*(m.raw)(m2.raw)
         Matrix[Double](newArr, m.rows, m.cols, m.rowStride, m.colStride, m.offset)(using BoundsCheck.DoBoundsCheck.no)
-      else ???
+      else
+        // Different memory layouts: materialize one matrix to match the other's layout
+        if m.isDenseColMajor then
+          // m is dense column-major, materialize m2 to column-major and multiply in-place
+          val m2Dense = m2.deepCopy(asRowMajor = false)
+          vecxt.arrays.*=(m2Dense.raw)(m.raw)
+          m2Dense
+        else if m.isDenseRowMajor then
+          // m is dense row-major, materialize m2 to row-major and multiply in-place
+          val m2Dense = m2.deepCopy(asRowMajor = true)
+          vecxt.arrays.*=(m2Dense.raw)(m.raw)
+          m2Dense
+        else if m2.isDenseColMajor then
+          // m2 is dense column-major, materialize m to column-major and multiply in-place
+          val mDense = m.deepCopy(asRowMajor = false)
+          vecxt.arrays.*=(mDense.raw)(m2.raw)
+          mDense
+        else if m2.isDenseRowMajor then
+          // m2 is dense row-major, materialize m to row-major and multiply in-place
+          val mDense = m.deepCopy(asRowMajor = true)
+          vecxt.arrays.*=(mDense.raw)(m2.raw)
+          mDense
+        else
+          // Neither is dense, materialize both to column-major and use SIMD multiplication
+          val mDense = m.deepCopy(asRowMajor = false)
+          val m2Dense = m2.deepCopy(asRowMajor = false)
+          val newArr = vecxt.arrays.*(mDense.raw)(m2Dense.raw)
+          Matrix[Double](newArr, m.rows, m.cols)(using BoundsCheck.DoBoundsCheck.no)
+        end if
       end if
     end hadamard
 

vecxt/test/src/matrix.test.scala

@@ -565,6 +565,75 @@ class MatrixExtensionSuite extends FunSuite:
     assertVecEquals[Double](result.raw, NArray[Double](2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0))
   }
 
+  test("hadamard product with non-simple layout") {
+    // Create base matrices
+    val base1 = Matrix[Double](NArray.tabulate[Double](9)(i => (i + 1).toDouble), 3, 3)
+    val base2 = Matrix[Double](NArray.tabulate[Double](9)(i => (i + 10).toDouble), 3, 3)
+
+    // Create views with non-simple layouts
+    val view1 = base1(::, NArray(1, 2)) // columns 1 and 2
+    val view2 = base2(::, NArray(1, 2)) // columns 1 and 2
+
+    // Compute hadamard product
+    val result = view1.hadamard(view2)
+
+    // Expected values (column-major):
+    // view1 columns 1,2 of base1 = [[4,7], [5,8], [6,9]]
+    // view2 columns 1,2 of base2 = [[13,16], [14,17], [15,18]]
+    // hadamard should give element-wise multiplication
+    assert(result.rows == 3)
+    assert(result.cols == 2)
+    assertEqualsDouble(result(0, 0), 4.0 * 13.0, 0.0001) // 52
+    assertEqualsDouble(result(1, 0), 5.0 * 14.0, 0.0001) // 70
+    assertEqualsDouble(result(2, 0), 6.0 * 15.0, 0.0001) // 90
+    assertEqualsDouble(result(0, 1), 7.0 * 16.0, 0.0001) // 112
+    assertEqualsDouble(result(1, 1), 8.0 * 17.0, 0.0001) // 136
+    assertEqualsDouble(result(2, 1), 9.0 * 18.0, 0.0001) // 162
+  }
+
+  test("hadamard product with mixed layouts") {
+    // One matrix with simple layout, one with non-simple
+    val simple = Matrix[Double](NArray(1.0, 2.0, 3.0, 4.0, 5.0, 6.0), 3, 2)
+    val base = Matrix[Double](NArray.tabulate[Double](9)(i => (i + 10).toDouble), 3, 3)
+    val view = base(::, NArray(0, 2)) // columns 0 and 2
+
+    val result = simple.hadamard(view)
+
+    assert(result.rows == 3)
+    assert(result.cols == 2)
+    assertEqualsDouble(result(0, 0), 1.0 * 10.0, 0.0001)
+    assertEqualsDouble(result(1, 0), 2.0 * 11.0, 0.0001)
+    assertEqualsDouble(result(2, 0), 3.0 * 12.0, 0.0001)
+    assertEqualsDouble(result(0, 1), 4.0 * 16.0, 0.0001)
+    assertEqualsDouble(result(1, 1), 5.0 * 17.0, 0.0001)
+    assertEqualsDouble(result(2, 1), 6.0 * 18.0, 0.0001)
+  }
+
+  test("hadamard product with transposed matrix") {
+    val mat1 = Matrix[Double](NArray(1.0, 2.0, 3.0, 4.0, 5.0, 6.0), 2, 3)
+    val mat2 = Matrix[Double](NArray(10.0, 20.0, 30.0, 40.0, 50.0, 60.0), 3, 2)
+
+    // Transpose mat2 to have the same shape as mat1
+    val mat2T = mat2.transpose
+
+    assert(mat1.rows == mat2T.rows)
+    assert(mat1.cols == mat2T.cols)
+
+    val result = mat1.hadamard(mat2T)
+
+    assert(result.rows == 2)
+    assert(result.cols == 3)
+    // mat1 (col-major 2x3): [[1,3,5], [2,4,6]]
+    // mat2 (col-major 3x2): [[10,40], [20,50], [30,60]]
+    // mat2T (row-major 2x3): [[10,20,30], [40,50,60]]
+    assertEqualsDouble(result(0, 0), 1.0 * 10.0, 0.0001) // 10
+    assertEqualsDouble(result(1, 0), 2.0 * 40.0, 0.0001) // 80
+    assertEqualsDouble(result(0, 1), 3.0 * 20.0, 0.0001) // 60
+    assertEqualsDouble(result(1, 1), 4.0 * 50.0, 0.0001) // 200
+    assertEqualsDouble(result(0, 2), 5.0 * 30.0, 0.0001) // 150
+    assertEqualsDouble(result(1, 2), 6.0 * 60.0, 0.0001) // 360
+  }
+
   test("map rows") {
     val mapped = mat1to9.mapRows[Double](row => row * 2)
     assertVecEquals[Double](mapped.raw, mat1to9.raw * 2)

vecxt/test/src/notImpl.test.scala

@@ -22,9 +22,6 @@ class NotImplTest extends FunSuite:
     intercept[NotImplementedError](m / 2.0)
   }
 
-  test("DoubleMatrix.hadamard with non-simple layout throws") {
-    intercept[NotImplementedError](m.hadamard(m))
-  }
   test("DoubleMatrix./:/ with non-simple layout throws") {
     intercept[NotImplementedError](m./:/(m))
   }