diff --git a/CHANGELOG.md b/CHANGELOG.md
index 74a11fc3d..89b06e3cf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### Added
 - Convenient matrix builders for rows, columns, vstacks and hstacks
+- Spreadsheet matrix builder
 
 ### Changed
 
diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt
index 64efaf640..824a7af77 100644
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt
+++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt
@@ -9,8 +9,8 @@ import kotlinx.benchmark.Benchmark
 import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
+import space.kscience.kmath.linear.MatrixBuilder
 import space.kscience.kmath.linear.linearSpace
-import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.linear.symmetric
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.tensors.core.symEigJacobi
@@ -24,7 +24,7 @@ internal class TensorAlgebraBenchmark {
         private val random = Random(12224)
         private const val dim = 30
 
-        private val matrix = Float64Field.linearSpace.matrix(dim, dim).symmetric { _, _ -> random.nextDouble() }
+        private val matrix = Float64Field.linearSpace.MatrixBuilder(dim, dim).symmetric { _, _ -> random.nextDouble() }
     }
 
     @Benchmark
diff --git a/build.gradle.kts b/build.gradle.kts
index 987dc9764..781f5822c 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -14,7 +14,7 @@ allprojects {
     }
 
     group = "space.kscience"
-    version = "0.4.1"
+    version = "0.4.2-dev"
 }
 
 dependencies{
diff --git a/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt b/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
index 6c87491eb..0505fc1e4 100644
--- a/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
@@ -6,22 +6,21 @@
 package space.kscience.kmath.operations
 
 import space.kscience.kmath.commons.linear.CMLinearSpace
-import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.linear.MatrixBuilder
+import space.kscience.kmath.linear.fill
 import space.kscience.kmath.nd.Float64BufferND
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.mutableStructureND
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.viktor.viktorAlgebra
-import kotlin.collections.component1
-import kotlin.collections.component2
 
 fun main() {
     val viktorStructure = Float64Field.viktorAlgebra.mutableStructureND(2, 2) { (i, j) ->
         if (i == j) 2.0 else 0.0
     }
 
-    val cmMatrix: Structure2D<Float64> = CMLinearSpace.matrix(2, 2)(0.0, 1.0, 0.0, 3.0)
+    val cmMatrix: Structure2D<Float64> = CMLinearSpace.MatrixBuilder(2, 2).fill(0.0, 1.0, 0.0, 3.0)
 
     val res: Float64BufferND = Float64Field.ndAlgebra {
         exp(viktorStructure) + 2.0 * cmMatrix
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt
index 0b3ac56fd..ff77c366e 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt
@@ -5,7 +5,7 @@
 
 package space.kscience.kmath.linear
 
-import space.kscience.attributes.FlagAttribute
+import space.kscience.attributes.Attributes
 import space.kscience.attributes.SafeType
 import space.kscience.attributes.WithType
 import space.kscience.kmath.UnstableKMathAPI
@@ -13,62 +13,51 @@ import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.structures.BufferAccessor2D
 import space.kscience.kmath.structures.MutableBufferFactory
 
-public class MatrixBuilder<T : Any, out A : Ring<T>>(
+/**
+ * A builder for matrix with fixed size
+ */
+@UnstableKMathAPI
+public class MatrixBuilder<T, out A : Ring<T>>(
     public val linearSpace: LinearSpace<T, A>,
-    public val rows: Int,
-    public val columns: Int,
+    public val rowNum: Int,
+    public val colNum: Int,
 ) : WithType<T> {
 
     override val type: SafeType<T> get() = linearSpace.type
 
-    public operator fun invoke(vararg elements: T): Matrix<T> {
-        require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
-        return linearSpace.buildMatrix(rows, columns) { i, j -> elements[i * columns + j] }
-    }
+}
 
-    //TODO add specific matrix builder functions like diagonal, etc
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> MatrixBuilder<T, A>.sparse(): SparseMatrix<T> =
+    SparseMatrix(rowNum, colNum, linearSpace.elementAlgebra.zero)
+
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> MatrixBuilder<T, A>.fill(vararg elements: T): Matrix<T> {
+    require(rowNum * colNum == elements.size) { "The number of elements ${elements.size} is not equal $rowNum * $colNum" }
+    return linearSpace.buildMatrix(rowNum, colNum) { i, j -> elements[i * colNum + j] }
 }
 
 /**
  * Create a matrix builder with given number of rows and columns
  */
 @UnstableKMathAPI
-public fun <T : Any, A : Ring<T>> LinearSpace<T, A>.matrix(rows: Int, columns: Int): MatrixBuilder<T, A> =
+public fun <T, A : Ring<T>> LinearSpace<T, A>.MatrixBuilder(rows: Int, columns: Int): MatrixBuilder<T, A> =
     MatrixBuilder(this, rows, columns)
 
-@UnstableKMathAPI
-public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Point<T> {
-    return buildVector(elements.size) { elements[it] }
-}
-
-public inline fun <T : Any> LinearSpace<T, Ring<T>>.row(
-    size: Int,
-    crossinline builder: (Int) -> T,
-): Matrix<T> = buildMatrix(1, size) { _, j -> builder(j) }
-
-public fun <T : Any> LinearSpace<T, Ring<T>>.row(vararg values: T): Matrix<T> = row(values.size, values::get)
-
-public inline fun <T : Any> LinearSpace<T, Ring<T>>.column(
-    size: Int,
-    crossinline builder: (Int) -> T,
-): Matrix<T> = buildMatrix(size, 1) { i, _ -> builder(i) }
-
-public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)
-
-public object Symmetric : MatrixAttribute<Unit>, FlagAttribute
-
 /**
- * Naive implementation of a symmetric matrix builder, that adds a [Symmetric] tag. The resulting matrix contains
- * full `size^2` number of elements, but caches elements during calls to save [builder] calls. [builder] is always called in the
- * upper triangle region meaning that `i <= j`
+ * Naive implementation of a symmetric matrix builder, that adds a [Symmetric] tag.
+ * The resulting matrix contains full `size^2` number of elements,
+ * but caches elements during calls to save [builder] calls.
+ * Always called in the upper triangle region meaning that `i <= j`
  */
-public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
-    builder: (i: Int, j: Int) -> T,
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
+    builder: A.(i: Int, j: Int) -> T,
 ): Matrix<T> {
-    require(columns == rows) { "In order to build symmetric matrix, number of rows $rows should be equal to number of columns $columns" }
-    return with(BufferAccessor2D<T?>(rows, rows, MutableBufferFactory(type))) {
+    require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
+    return with(BufferAccessor2D<T?>(rowNum, rowNum, MutableBufferFactory(type))) {
         val cache = HashMap<IntArray, T>()
-        linearSpace.buildMatrix(rows, rows) { i, j ->
+        linearSpace.buildMatrix(this@symmetric.rowNum, this@symmetric.rowNum) { i, j ->
             val index = intArrayOf(i, j)
             val cached = cache[index]
             if (cached == null) {
@@ -81,4 +70,50 @@ public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
             }
         }.withAttribute(Symmetric)
     }
-}
\ No newline at end of file
+}
+
+/**
+ * Create a diagonal matrix with given factory.
+ */
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> MatrixBuilder<T, A>.diagonal(
+    builder: A.(Int) -> T
+): Matrix<T> = with(linearSpace.elementAlgebra) {
+    require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
+    return VirtualMatrix(rowNum, colNum, attributes = Attributes(IsDiagonal)) { i, j ->
+        check(i in 0 until rowNum) { "$i out of bounds: 0..<$rowNum" }
+        check(j in 0 until colNum) { "$j out of bounds: 0..<$colNum" }
+        if (i == j) {
+            builder(i)
+        } else {
+            zero
+        }
+    }
+}
+
+/**
+ * Create a diagonal matrix from elements
+ */
+@UnstableKMathAPI
+public fun <T> MatrixBuilder<T, Ring<T>>.diagonal(vararg elements: T): Matrix<T> {
+    require(colNum == rowNum) { "In order to build symmetric matrix, number of rows $rowNum should be equal to number of columns $colNum" }
+    return return VirtualMatrix(rowNum, colNum, attributes = Attributes(IsDiagonal)) { i, j ->
+        check(i in 0 until rowNum) { "$i out of bounds: 0..<$rowNum" }
+        check(j in 0 until colNum) { "$j out of bounds: 0..<$colNum" }
+        if (i == j) {
+            elements[i]
+        } else {
+            linearSpace.elementAlgebra.zero
+        }
+    }
+}
+
+/**
+ * Create a lazily evaluated virtual matrix with a given size
+ */
+@UnstableKMathAPI
+public fun <T : Any> MatrixBuilder<T, *>.virtual(
+    attributes: Attributes = Attributes.EMPTY,
+    generator: (i: Int, j: Int) -> T,
+): VirtualMatrix<T> = VirtualMatrix(rowNum, colNum, attributes, generator)
+
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/SparseMatrix.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/SparseMatrix.kt
new file mode 100644
index 000000000..601924ecd
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/SparseMatrix.kt
@@ -0,0 +1,69 @@
+/*
+ * Copyright 2018-2025 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.attributes.Attributes
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.operations.Ring
+
+/**
+ * Mutable sparse matrix that stores values only for non-zero cells ([DOK format](https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_(DOK))).
+ *
+ * [SparseMatrix] is ineffective, but does not depend on particular [LinearSpace]
+ *
+ * Using this class is almost always a [PerformancePitfall]. It should be used only for special cases like manual matrix building.
+ */
+@UnstableKMathAPI
+public class SparseMatrix<T>(
+    override val rowNum: Int,
+    override val colNum: Int,
+    private val zero: T,
+    cells: Map<Pair<Int, Int>, T> = emptyMap(),
+    override val attributes: Attributes = Attributes.EMPTY,
+) : MutableMatrix<T> {
+
+    private val cells = cells.toMutableMap()
+
+    override fun get(i: Int, j: Int): T {
+        if (i !in 0 until rowNum) throw IndexOutOfBoundsException("Row index $i out of row range 0..$rowNum")
+        if (j !in 0 until colNum) throw IndexOutOfBoundsException("Column index $j out of column range 0..$colNum")
+        return cells[i to j] ?: zero
+    }
+
+    override fun set(i: Int, j: Int, value: T) {
+        require(i in 0 until rowNum) { "Row index $i is out of bounds: 0..<$rowNum" }
+        require(j in 0 until colNum) { "Row index $j is out of bounds: 0..<$colNum" }
+        val coordinates = i to j
+        if (cells[coordinates] != null || value != zero) {
+            cells[coordinates] = value
+        }
+    }
+
+    override fun set(index: IntArray, value: T) {
+        require(index.size == 2) { "Index array must contain two elements." }
+        set(index[0], index[1], value)
+    }
+}
+
+@UnstableKMathAPI
+public fun <T> SparseMatrix<T>.fill(vararg elements: T): SparseMatrix<T> {
+    require(rowNum * colNum == elements.size) { "The number of elements ${elements.size} is not equal $rowNum * $colNum" }
+    for (i in 0 until rowNum) {
+        for (j in 0 until colNum) {
+            set(i, j, elements[i * rowNum + j])
+        }
+    }
+    return this
+}
+
+@UnstableKMathAPI
+public fun <T> LinearSpace<T, Ring<T>>.sparse(
+    rows: Int,
+    columns: Int,
+    attributes: Attributes = Attributes.EMPTY,
+    builder: SparseMatrix<T>.() -> Unit = {}
+): SparseMatrix<T> = SparseMatrix(rows, columns, elementAlgebra.zero, attributes = attributes).apply(builder)
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt
index 2b031de06..2809db9fd 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt
@@ -6,7 +6,6 @@
 package space.kscience.kmath.linear
 
 import space.kscience.attributes.Attributes
-import space.kscience.kmath.nd.ShapeND
 
 
 /**
@@ -20,13 +19,5 @@ public class VirtualMatrix<out T>(
     override val attributes: Attributes = Attributes.EMPTY,
     public val generator: (i: Int, j: Int) -> T,
 ) : Matrix<T> {
-
-    override val shape: ShapeND get() = ShapeND(rowNum, colNum)
-
     override operator fun get(i: Int, j: Int): T = generator(i, j)
-}
-
-public fun <T : Any> MatrixBuilder<T, *>.virtual(
-    attributes: Attributes = Attributes.EMPTY,
-    generator: (i: Int, j: Int) -> T,
-): VirtualMatrix<T> = VirtualMatrix(rows, columns, attributes, generator)
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
index 59ab009ec..a1dcc9c46 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
@@ -23,10 +23,17 @@ public interface MatrixScope<T> : WithType<T>
  */
 public interface MatrixAttribute<T> : StructureAttribute<T>
 
+/**
+ * Matrices with this feature are symmetric, meaning `matrix[i,j] == matrix[j,i]`
+ */
+public interface Symmetric : MatrixAttribute<Unit>, FlagAttribute{
+    public companion object: Symmetric
+}
+
 /**
  * Matrices with this feature are considered to have only diagonal non-zero elements.
  */
-public interface IsDiagonal : MatrixAttribute<Unit>, FlagAttribute {
+public interface IsDiagonal : Symmetric {
     public companion object : IsDiagonal
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixBuilders.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixBuilders.kt
new file mode 100644
index 000000000..76f95fb86
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixBuilders.kt
@@ -0,0 +1,81 @@
+/*
+ * Copyright 2018-2025 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.operations.Ring
+
+
+/**
+ * Create a vector from elements
+ */
+public fun <T> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Point<T> =
+    buildVector(elements.size) { elements[it] }
+
+/**
+ * Create a single row matrix
+ */
+public inline fun <T, A : Ring<T>> LinearSpace<T, A>.row(
+    size: Int,
+    crossinline builder: A.(Int) -> T,
+): Matrix<T> = buildMatrix(1, size) { _, j -> builder(j) }
+
+/**
+ * Create a single row matrix from elements
+ */
+public fun <T> LinearSpace<T, Ring<T>>.row(vararg elements: T): Matrix<T> = row(elements.size) { elements[it] }
+
+/**
+ * Create a single column matrix
+ */
+public inline fun <T, A : Ring<T>> LinearSpace<T, A>.column(
+    size: Int,
+    crossinline builder: A.(Int) -> T,
+): Matrix<T> = buildMatrix(size, 1) { i, _ -> builder(i) }
+
+/**
+ * Create a single column matrix from elements
+ */
+public fun <T> LinearSpace<T, Ring<T>>.column(vararg elements: T): Matrix<T> = column(elements.size) { elements[it] }
+
+/**
+ * Stack vertically several matrices with the same number of columns.
+ *
+ * Resulting matrix number of rows is the sum of rows in all [matrices]
+ */
+@PerformancePitfall
+@UnstableKMathAPI
+public fun <T> LinearSpace<T, Ring<T>>.vstack(vararg matrices: Matrix<T>): Matrix<T> {
+    require(matrices.isNotEmpty()) { "No matrices" }
+    val colNum = matrices.first().colNum
+    require(matrices.all { it.colNum == colNum }) { "All matrices must have the same number of columns: $colNum" }
+
+    val rows = matrices.flatMap { it.rows }
+
+    return buildMatrix(matrices.sumOf { it.rowNum }, colNum) { row, column->
+        rows[row][column]
+    }
+}
+
+/**
+ * Stack horizontally several matrices with the same number of rows.
+ *
+ * Resulting matrix number of co is the sum of rows in all [matrices]
+ */
+@PerformancePitfall
+@UnstableKMathAPI
+public fun <T> LinearSpace<T, Ring<T>>.hstack(vararg matrices: Matrix<T>): Matrix<T> {
+    require(matrices.isNotEmpty()) { "No matrices" }
+    val rowNum = matrices.first().rowNum
+    require(matrices.all { it.rowNum == rowNum }) { "All matrices must have the same number of rows: $rowNum" }
+
+    val columns = matrices.flatMap { it.columns }
+
+    return buildMatrix(rowNum, matrices.sumOf { it.colNum }) { row, column->
+        columns[column][row]
+    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
index be2856c3b..a2031dbad 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
@@ -228,7 +228,6 @@ public interface MutableStructureND<T> : StructureND<T> {
      * @param index the indices.
      * @param value the value.
      */
-    @PerformancePitfall
     public operator fun set(index: IntArray, value: T)
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/LogicAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/LogicAlgebra.kt
index c4a26e242..bf1661479 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/LogicAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/LogicAlgebra.kt
@@ -10,7 +10,7 @@ import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.structures.MutableBufferFactory
 
 /**
- * An algebra for generic boolean logic
+ * Algebra for generic boolean logic
  */
 @UnstableKMathAPI
 public interface LogicAlgebra<T : Any> : Algebra<T> {
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt
index 32990c707..d899e3075 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt
@@ -22,7 +22,7 @@ internal class BufferAccessor2D<T>(
     }
 
     inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer<T> =
-        factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
+        factory(this@BufferAccessor2D.rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
 
     fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
 
diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt
index 037fdf3e5..1c9a9ad29 100644
--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt
@@ -29,7 +29,7 @@ class DoubleLUSolverTest {
 
     @Test
     fun testDecomposition() = with(Double.algebra.linearSpace) {
-        val matrix = matrix(2, 2)(
+        val matrix = MatrixBuilder(2, 2).fill(
             3.0, 1.0,
             2.0, 3.0
         )
@@ -44,14 +44,14 @@ class DoubleLUSolverTest {
 
     @Test
     fun testInvert() = Double.algebra.linearSpace.run {
-        val matrix = matrix(2, 2)(
+        val matrix = MatrixBuilder(2, 2).fill(
             3.0, 1.0,
             1.0, 3.0
         )
 
         val inverted = lupSolver().inverse(matrix)
 
-        val expected = matrix(2, 2)(
+        val expected = MatrixBuilder(2, 2).fill(
             0.375, -0.125,
             -0.125, 0.375
         )
diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixBuilderTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixBuilderTest.kt
new file mode 100644
index 000000000..3b53f4361
--- /dev/null
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixBuilderTest.kt
@@ -0,0 +1,38 @@
+/*
+ * Copyright 2018-2025 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+@UnstableKMathAPI
+class MatrixBuilderTest {
+
+    @Test
+    fun buildCompositeMatrix() = with(Float64.algebra.linearSpace) {
+
+        val matrix = vstack(
+            sparse(1, 5) { set(0, 4, 1.0) },
+            hstack(
+                sparse(4, 4).fill(
+                    1.0, 1.0, 0.0, 0.0,
+                    0.0, 1.0, 1.0, 0.0,
+                    0.0, 0.0, 1.0, 1.0,
+                    0.0, 0.0, 0.0, 1.0
+                ),
+                sparse(4, 1)
+            )
+        )
+
+        println(StructureND.toString(matrix))
+
+        assertEquals(1.0, matrix[0, 4])
+    }
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt
index 3489e92a4..bd3486767 100644
--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt
@@ -29,7 +29,7 @@ class MatrixTest {
 
     @Test
     fun testBuilder() = Double.algebra.linearSpace.run {
-        val matrix = matrix(2, 3)(
+        val matrix = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
diff --git a/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt b/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt
index 19a37c9fb..4725d606e 100644
--- a/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt
+++ b/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt
@@ -28,7 +28,7 @@ class ParallelMatrixTest {
 
     @Test
     fun testBuilder() = Float64Field.linearSpace.parallel {
-        val matrix = matrix(2, 3)(
+        val matrix = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt
index 39933076a..c8b26369f 100644
--- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt
+++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt
@@ -37,8 +37,8 @@ public fun realMatrix(rowNum: Int, colNum: Int, initializer: Float64Field.(i: In
     Double.algebra.linearSpace.buildMatrix(rowNum, colNum, initializer)
 
 @OptIn(UnstableKMathAPI::class)
-public fun realMatrix(rowNum: Int, colNum: Int): MatrixBuilder<Double, Float64Field> =
-    Double.algebra.linearSpace.matrix(rowNum, colNum)
+public fun MatrixBuilder(rowNum: Int, colNum: Int): MatrixBuilder<Double, Float64Field> =
+    Double.algebra.linearSpace.MatrixBuilder(rowNum, colNum)
 
 public fun Array<DoubleArray>.toMatrix(): RealMatrix {
     return Double.algebra.linearSpace.buildMatrix(size, this[0].size) { row, col -> this@toMatrix[row][col] }
diff --git a/kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleMatrixTest.kt b/kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleMatrixTest.kt
index 69ed9574d..921e36767 100644
--- a/kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleMatrixTest.kt
+++ b/kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleMatrixTest.kt
@@ -7,8 +7,9 @@ package space.kscience.kmath.real
 
 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.linear.MatrixBuilder
+import space.kscience.kmath.linear.fill
 import space.kscience.kmath.linear.linearSpace
-import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.testutils.contentEquals
@@ -43,11 +44,11 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testRepeatStackVertical() {
-        val matrix1 = realMatrix(2, 3)(
+        val matrix1 = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
-        val matrix2 = realMatrix(6, 3)(
+        val matrix2 = MatrixBuilder(6, 3).fill(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0,
             1.0, 0.0, 0.0,
@@ -60,12 +61,12 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testMatrixAndDouble() = Double.algebra.linearSpace.run {
-        val matrix1 = realMatrix(2, 3)(
+        val matrix1 = space.kscience.kmath.real.MatrixBuilder(2, 3).fill(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = (matrix1 * 2.5 + 1.0 - 2.0) / 2.0
-        val expectedResult = matrix(2, 3)(
+        val expectedResult = this.MatrixBuilder(2, 3).fill(
             0.75, -0.5, 3.25,
             4.5, 7.0, 2.0
         )
@@ -74,13 +75,13 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testDoubleAndMatrix() {
-        val matrix1 = realMatrix(2, 3)(
+        val matrix1 = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = 20.0 - (10.0 + (5.0 * matrix1))
         //val matrix2 = 10.0 + (5.0 * matrix1)
-        val expectedResult = realMatrix(2, 3)(
+        val expectedResult = MatrixBuilder(2, 3).fill(
             5.0, 10.0, -5.0,
             -10.0, -20.0, 0.0
         )
@@ -89,15 +90,15 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testSquareAndPower() {
-        val matrix1 = realMatrix(2, 3)(
+        val matrix1 = MatrixBuilder(2, 3).fill(
             -1.0, 0.0, 3.0,
             4.0, -6.0, -2.0
         )
-        val matrix2 = realMatrix(2, 3)(
+        val matrix2 = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 9.0,
             16.0, 36.0, 4.0
         )
-        val matrix3 = realMatrix(2, 3)(
+        val matrix3 = MatrixBuilder(2, 3).fill(
             -1.0, 0.0, 27.0,
             64.0, -216.0, -8.0
         )
@@ -108,16 +109,16 @@ internal class DoubleMatrixTest {
     @OptIn(UnstableKMathAPI::class)
     @Test
     fun testTwoMatrixOperations() {
-        val matrix1 = realMatrix(2, 3)(
+        val matrix1 = MatrixBuilder(2, 3).fill(
             -1.0, 0.0, 3.0,
             4.0, -6.0, 7.0
         )
-        val matrix2 = realMatrix(2, 3)(
+        val matrix2 = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 3.0,
             4.0, 6.0, -2.0
         )
         val result = matrix1 * matrix2 + matrix1 - matrix2
-        val expectedResult = realMatrix(2, 3)(
+        val expectedResult = MatrixBuilder(2, 3).fill(
             -3.0, 0.0, 9.0,
             16.0, -48.0, -5.0
         )
@@ -126,16 +127,16 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testColumnOperations() {
-        val matrix1 = realMatrix(2, 4)(
+        val matrix1 = MatrixBuilder(2, 4).fill(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )
-        val matrix2 = realMatrix(2, 5)(
+        val matrix2 = MatrixBuilder(2, 5).fill(
             -1.0, 0.0, 3.0, 15.0, -1.0,
             4.0, -6.0, 7.0, -11.0, 4.0
         )
-        val col1 = realMatrix(2, 1)(0.0, -6.0)
-        val cols1to2 = realMatrix(2, 2)(
+        val col1 = MatrixBuilder(2, 1).fill(0.0, -6.0)
+        val cols1to2 = MatrixBuilder(2, 2).fill(
             0.0, 3.0,
             -6.0, 7.0
         )
@@ -160,7 +161,7 @@ internal class DoubleMatrixTest {
 
     @Test
     fun testAllElementOperations() = Double.algebra.linearSpace.run {
-        val matrix1 = matrix(2, 4)(
+        val matrix1 = this.MatrixBuilder(2, 4).fill(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )
diff --git a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/rotations3D.kt b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/rotations3D.kt
index 78fe01733..027632b24 100644
--- a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/rotations3D.kt
+++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/rotations3D.kt
@@ -8,10 +8,7 @@ package space.kscience.kmath.geometry.euclidean3d
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.complex.*
 import space.kscience.kmath.geometry.*
-import space.kscience.kmath.linear.LinearSpace
-import space.kscience.kmath.linear.Matrix
-import space.kscience.kmath.linear.linearSpace
-import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.linear.*
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.structures.Float64
 import kotlin.math.*
@@ -103,7 +100,7 @@ public fun Quaternion.toRotationMatrix(
     linearSpace: LinearSpace<Double, *> = Float64Field.linearSpace,
 ): Matrix<Float64> {
     val s = QuaternionAlgebra.norm(this).pow(-2)
-    return linearSpace.matrix(3, 3)(
+    return linearSpace.MatrixBuilder(3, 3).fill(
         1.0 - 2 * s * (y * y + z * z), 2 * s * (x * y - z * w), 2 * s * (x * z + y * w),
         2 * s * (x * y + z * w), 1.0 - 2 * s * (x * x + z * z), 2 * s * (y * z - x * w),
         2 * s * (x * z - y * w), 2 * s * (y * z + x * w), 1.0 - 2 * s * (x * x + y * y)
diff --git a/kmath-ojalgo/build.gradle.kts b/kmath-ojalgo/build.gradle.kts
index c8c8a4abb..f4ab0e269 100644
--- a/kmath-ojalgo/build.gradle.kts
+++ b/kmath-ojalgo/build.gradle.kts
@@ -15,6 +15,10 @@ kscience {
 //        api(projects.kmathFunctions)
         api(libs.ojalgo)
     }
+
+    jvmTest{
+        implementation(projects.testUtils)
+    }
 }
 
 readme {
diff --git a/kmath-ojalgo/src/jvmTest/kotlin/space/kscience/kmath/ojalgo/OjalgoMatrixTest.kt b/kmath-ojalgo/src/jvmTest/kotlin/space/kscience/kmath/ojalgo/OjalgoMatrixTest.kt
index 803064cad..ac614e672 100644
--- a/kmath-ojalgo/src/jvmTest/kotlin/space/kscience/kmath/ojalgo/OjalgoMatrixTest.kt
+++ b/kmath-ojalgo/src/jvmTest/kotlin/space/kscience/kmath/ojalgo/OjalgoMatrixTest.kt
@@ -30,7 +30,7 @@ class OjalgoMatrixTest {
 
     @Test
     fun testBuilder() = Ojalgo.Companion.R064.linearSpace {
-        val matrix = buildMatrix(2, 3)(
+        val matrix = MatrixBuilder(2, 3).fill(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
@@ -77,7 +77,7 @@ class OjalgoMatrixTest {
 
     @Test
     fun testCholesky() = with(Ojalgo.Companion.R064.linearSpace) {
-        val l = buildMatrix(4, 4)(
+        val l = MatrixBuilder(4, 4).fill(
             1.0, 0.0, 0.0, 0.0,
             1.0, 1.0, 0.0, 0.0,
             1.0, 1.0, 1.0, 0.0,
diff --git a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt
index 337904923..d561ad776 100644
--- a/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt
+++ b/kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/QowOptimizer.kt
@@ -91,7 +91,7 @@ public object QowOptimizer : Optimizer<Double, XYFit> {
      * D(\phi)=E(\phi_k(\theta_0) \phi_l(\theta_0))= disDeriv_k * disDeriv_l /sigma^2
      */
     private fun QoWeight.covarF(): Matrix<Float64> =
-        linearSpace.matrix(size, size).symmetric { s1, s2 ->
+        linearSpace.MatrixBuilder(size, size).symmetric { s1, s2 ->
             (0 until data.size).sumOf { d -> derivs[d, s1] * derivs[d, s2] / dispersion[d] }
         }