diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt
similarity index 73%
rename from examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt
rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt
index 8823e86db..9d2b02245 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt
@@ -2,22 +2,19 @@ package kscience.kmath.benchmarks
 
 import kotlinx.benchmark.Benchmark
 import kscience.kmath.commons.linear.CMMatrixContext
+import kscience.kmath.commons.linear.CMMatrixContext.dot
 import kscience.kmath.commons.linear.toCM
 import kscience.kmath.ejml.EjmlMatrixContext
 import kscience.kmath.ejml.toEjml
-import kscience.kmath.linear.BufferMatrixContext
-import kscience.kmath.linear.RealMatrixContext
 import kscience.kmath.linear.real
-import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Buffer
 import kscience.kmath.structures.Matrix
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
 import kotlin.random.Random
 
 @State(Scope.Benchmark)
-class DotBenchmark {
+class MultiplicationBenchmark {
     companion object {
         val random = Random(12224)
         val dim = 1000
@@ -35,14 +32,14 @@ class DotBenchmark {
 
     @Benchmark
     fun commonsMathMultiplication() {
-        CMMatrixContext {
+        CMMatrixContext.invoke {
             cmMatrix1 dot cmMatrix2
         }
     }
 
     @Benchmark
     fun ejmlMultiplication() {
-        EjmlMatrixContext {
+        EjmlMatrixContext.invoke {
             ejmlMatrix1 dot ejmlMatrix2
         }
     }
@@ -51,22 +48,13 @@ class DotBenchmark {
     fun ejmlMultiplicationwithConversion() {
         val ejmlMatrix1 = matrix1.toEjml()
         val ejmlMatrix2 = matrix2.toEjml()
-        EjmlMatrixContext {
+        EjmlMatrixContext.invoke {
             ejmlMatrix1 dot ejmlMatrix2
         }
     }
 
     @Benchmark
     fun bufferedMultiplication() {
-        BufferMatrixContext(RealField, Buffer.Companion::real).invoke{
-            matrix1 dot matrix2
-        }
-    }
-
-    @Benchmark
-    fun realMultiplication(){
-        RealMatrixContext {
-            matrix1 dot matrix2
-        }
+        matrix1 dot matrix2
     }
 }
\ No newline at end of file
diff --git a/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt b/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt
index 96684f7dc..987eea16f 100644
--- a/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt
+++ b/examples/src/main/kotlin/kscience/kmath/structures/typeSafeDimensions.kt
@@ -4,8 +4,9 @@ import kscience.kmath.dimensions.D2
 import kscience.kmath.dimensions.D3
 import kscience.kmath.dimensions.DMatrixContext
 import kscience.kmath.dimensions.Dimension
+import kscience.kmath.operations.RealField
 
-private fun DMatrixContext<Double>.simple() {
+private fun DMatrixContext<Double, RealField>.simple() {
     val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
     val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
 
@@ -17,7 +18,7 @@ private object D5 : Dimension {
     override val dim: UInt = 5u
 }
 
-private fun DMatrixContext<Double>.custom() {
+private fun DMatrixContext<Double, RealField>.custom() {
     val m1 = produce<D2, D5> { i, j -> (i + j).toDouble() }
     val m2 = produce<D5, D2> { i, j -> (i - j).toDouble() }
     val m3 = produce<D2, D2> { i, j -> (i - j).toDouble() }
diff --git a/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MST.kt b/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MST.kt
index 6cf746722..212fd0d0b 100644
--- a/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MST.kt
+++ b/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MST.kt
@@ -2,10 +2,9 @@ package kscience.kmath.ast
 
 import kscience.kmath.operations.Algebra
 import kscience.kmath.operations.NumericAlgebra
-import kscience.kmath.operations.RealField
 
 /**
- * A Mathematical Syntax Tree node for mathematical expressions.
+ * A Mathematical Syntax Tree (MST) node for mathematical expressions.
  *
  * @author Alexander Nozik
  */
@@ -57,21 +56,22 @@ public fun <T> Algebra<T>.evaluate(node: MST): T = when (node) {
         ?: error("Numeric nodes are not supported by $this")
 
     is MST.Symbolic -> symbol(node.value)
-    is MST.Unary -> unaryOperationFunction(node.operation)(evaluate(node.value))
+
+    is MST.Unary -> when {
+        this is NumericAlgebra && node.value is MST.Numeric -> unaryOperationFunction(node.operation)(number(node.value.value))
+        else -> unaryOperationFunction(node.operation)(evaluate(node.value))
+    }
 
     is MST.Binary -> when {
-        this !is NumericAlgebra -> binaryOperationFunction(node.operation)(evaluate(node.left), evaluate(node.right))
+        this is NumericAlgebra && node.left is MST.Numeric && node.right is MST.Numeric ->
+            binaryOperationFunction(node.operation)(number(node.left.value), number(node.right.value))
 
-        node.left is MST.Numeric && node.right is MST.Numeric -> {
-            val number = RealField
-                .binaryOperationFunction(node.operation)
-                .invoke(node.left.value.toDouble(), node.right.value.toDouble())
+        this is NumericAlgebra && node.left is MST.Numeric ->
+            leftSideNumberOperationFunction(node.operation)(node.left.value, evaluate(node.right))
 
-            number(number)
-        }
+        this is NumericAlgebra && node.right is MST.Numeric ->
+            rightSideNumberOperationFunction(node.operation)(evaluate(node.left), node.right.value)
 
-        node.left is MST.Numeric -> leftSideNumberOperationFunction(node.operation)(node.left.value, evaluate(node.right))
-        node.right is MST.Numeric -> rightSideNumberOperationFunction(node.operation)(evaluate(node.left), node.right.value)
         else -> binaryOperationFunction(node.operation)(evaluate(node.left), evaluate(node.right))
     }
 }
diff --git a/kmath-ast/src/jsMain/kotlin/kscience/kmath/estree/estree.kt b/kmath-ast/src/jsMain/kotlin/kscience/kmath/estree/estree.kt
index 159c5d5ec..5c08ada31 100644
--- a/kmath-ast/src/jsMain/kotlin/kscience/kmath/estree/estree.kt
+++ b/kmath-ast/src/jsMain/kotlin/kscience/kmath/estree/estree.kt
@@ -1,18 +1,18 @@
 package kscience.kmath.estree
 
 import kscience.kmath.ast.MST
+import kscience.kmath.ast.MST.*
 import kscience.kmath.ast.MstExpression
 import kscience.kmath.estree.internal.ESTreeBuilder
 import kscience.kmath.estree.internal.estree.BaseExpression
 import kscience.kmath.expressions.Expression
 import kscience.kmath.operations.Algebra
 import kscience.kmath.operations.NumericAlgebra
-import kscience.kmath.operations.RealField
 
 @PublishedApi
 internal fun <T> MST.compileWith(algebra: Algebra<T>): Expression<T> {
     fun ESTreeBuilder<T>.visit(node: MST): BaseExpression = when (node) {
-        is MST.Symbolic -> {
+        is Symbolic -> {
             val symbol = try {
                 algebra.symbol(node.value)
             } catch (ignored: IllegalStateException) {
@@ -25,25 +25,29 @@ internal fun <T> MST.compileWith(algebra: Algebra<T>): Expression<T> {
                 variable(node.value)
         }
 
-        is MST.Numeric -> constant(node.value)
-        is MST.Unary -> call(algebra.unaryOperationFunction(node.operation), visit(node.value))
+        is Numeric -> constant(node.value)
 
-        is MST.Binary -> when {
-            algebra is NumericAlgebra<T> && node.left is MST.Numeric && node.right is MST.Numeric -> constant(
-                algebra.number(
-                    RealField
-                        .binaryOperationFunction(node.operation)
-                        .invoke(node.left.value.toDouble(), node.right.value.toDouble())
-                )
+        is Unary -> when {
+            algebra is NumericAlgebra && node.value is Numeric -> constant(
+                algebra.unaryOperationFunction(node.operation)(algebra.number(node.value.value)))
+
+            else -> call(algebra.unaryOperationFunction(node.operation), visit(node.value))
+        }
+
+        is Binary -> when {
+            algebra is NumericAlgebra && node.left is Numeric && node.right is Numeric -> constant(
+                algebra
+                    .binaryOperationFunction(node.operation)
+                    .invoke(algebra.number(node.left.value), algebra.number(node.right.value))
             )
 
-            algebra is NumericAlgebra<T> && node.left is MST.Numeric -> call(
+            algebra is NumericAlgebra && node.left is Numeric -> call(
                 algebra.leftSideNumberOperationFunction(node.operation),
                 visit(node.left),
                 visit(node.right),
             )
 
-            algebra is NumericAlgebra<T> && node.right is MST.Numeric -> call(
+            algebra is NumericAlgebra && node.right is Numeric -> call(
                 algebra.rightSideNumberOperationFunction(node.operation),
                 visit(node.left),
                 visit(node.right),
diff --git a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt
index b98c0bfce..55cdec243 100644
--- a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt
+++ b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt
@@ -3,11 +3,11 @@ package kscience.kmath.asm
 import kscience.kmath.asm.internal.AsmBuilder
 import kscience.kmath.asm.internal.buildName
 import kscience.kmath.ast.MST
+import kscience.kmath.ast.MST.*
 import kscience.kmath.ast.MstExpression
 import kscience.kmath.expressions.Expression
 import kscience.kmath.operations.Algebra
 import kscience.kmath.operations.NumericAlgebra
-import kscience.kmath.operations.RealField
 
 /**
  * Compiles given MST to an Expression using AST compiler.
@@ -20,7 +20,7 @@ import kscience.kmath.operations.RealField
 @PublishedApi
 internal fun <T : Any> MST.compileWith(type: Class<T>, algebra: Algebra<T>): Expression<T> {
     fun AsmBuilder<T>.visit(node: MST): Unit = when (node) {
-        is MST.Symbolic -> {
+        is Symbolic -> {
             val symbol = try {
                 algebra.symbol(node.value)
             } catch (ignored: IllegalStateException) {
@@ -33,24 +33,29 @@ internal fun <T : Any> MST.compileWith(type: Class<T>, algebra: Algebra<T>): Exp
                 loadVariable(node.value)
         }
 
-        is MST.Numeric -> loadNumberConstant(node.value)
-        is MST.Unary -> buildCall(algebra.unaryOperationFunction(node.operation)) { visit(node.value) }
+        is Numeric -> loadNumberConstant(node.value)
 
-        is MST.Binary -> when {
-            algebra is NumericAlgebra<T> && node.left is MST.Numeric && node.right is MST.Numeric -> loadObjectConstant(
-                algebra.number(
-                    RealField
-                        .binaryOperationFunction(node.operation)
-                        .invoke(node.left.value.toDouble(), node.right.value.toDouble())
-                )
+        is Unary -> when {
+            algebra is NumericAlgebra && node.value is Numeric -> loadObjectConstant(
+                algebra.unaryOperationFunction(node.operation)(algebra.number(node.value.value)))
+
+            else -> buildCall(algebra.unaryOperationFunction(node.operation)) { visit(node.value) }
+        }
+
+        is Binary -> when {
+            algebra is NumericAlgebra && node.left is Numeric && node.right is Numeric -> loadObjectConstant(
+                algebra.binaryOperationFunction(node.operation)
+                    .invoke(algebra.number(node.left.value), algebra.number(node.right.value))
             )
 
-            algebra is NumericAlgebra<T> && node.left is MST.Numeric -> buildCall(algebra.leftSideNumberOperationFunction(node.operation)) {
+            algebra is NumericAlgebra && node.left is Numeric -> buildCall(
+                algebra.leftSideNumberOperationFunction(node.operation)) {
                 visit(node.left)
                 visit(node.right)
             }
 
-            algebra is NumericAlgebra<T> && node.right is MST.Numeric -> buildCall(algebra.rightSideNumberOperationFunction(node.operation)) {
+            algebra is NumericAlgebra && node.right is Numeric -> buildCall(
+                algebra.rightSideNumberOperationFunction(node.operation)) {
                 visit(node.left)
                 visit(node.right)
             }
diff --git a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt
index 1edbed28d..93d8d1143 100644
--- a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt
+++ b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt
@@ -191,7 +191,7 @@ internal class AsmBuilder<T>(
         }
 
         val cls = classLoader.defineClass(className, classWriter.toByteArray())
-        java.io.File("dump.class").writeBytes(classWriter.toByteArray())
+        // java.io.File("dump.class").writeBytes(classWriter.toByteArray())
         val l = MethodHandles.publicLookup()
 
         if (hasConstants)
diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
index 49888f8d6..712927400 100644
--- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
+++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
@@ -29,7 +29,6 @@ public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature
     }
 }
 
-//TODO move inside context
 public fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {
     this
 } else {
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
index 402161207..8b50bbe33 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
@@ -1,10 +1,8 @@
 package kscience.kmath.linear
 
+import kscience.kmath.operations.RealField
 import kscience.kmath.operations.Ring
-import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.BufferFactory
-import kscience.kmath.structures.NDStructure
-import kscience.kmath.structures.asSequence
+import kscience.kmath.structures.*
 
 /**
  * Basic implementation of Matrix space based on [NDStructure]
@@ -23,6 +21,24 @@ public class BufferMatrixContext<T : Any, R : Ring<T>>(
     public companion object
 }
 
+@Suppress("OVERRIDE_BY_INLINE")
+public object RealMatrixContext : GenericMatrixContext<Double, RealField, BufferMatrix<Double>> {
+    public override val elementContext: RealField
+        get() = RealField
+
+    public override inline fun produce(
+        rows: Int,
+        columns: Int,
+        initializer: (i: Int, j: Int) -> Double,
+    ): BufferMatrix<Double> {
+        val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
+        return BufferMatrix(rows, columns, buffer)
+    }
+
+    public override inline fun point(size: Int, initializer: (Int) -> Double): Point<Double> =
+        RealBuffer(size, initializer)
+}
+
 public class BufferMatrix<T : Any>(
     public override val rowNum: Int,
     public override val colNum: Int,
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
index bf2a9f59e..099fa1909 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
@@ -213,8 +213,17 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
     return decomposition.solveWithLUP(bufferFactory, b)
 }
 
+public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): FeaturedMatrix<Double> =
+    solveWithLUP(a, b) { it < 1e-11 }
+
 public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.inverseWithLUP(
     matrix: Matrix<T>,
     noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
     noinline checkSingular: (T) -> Boolean,
-): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
\ No newline at end of file
+): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
+
+/**
+ * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
+ */
+public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): FeaturedMatrix<Double> =
+    solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
index 8c28a240f..648af605c 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
@@ -18,11 +18,6 @@ public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Mat
      */
     public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
 
-    /**
-     * Produce a point compatible with matrix space (and possibly optimized for it)
-     */
-    public fun point(size: Int, initializer: (Int) -> T): Point<T> = Buffer.boxing(size, initializer)
-
     @Suppress("UNCHECKED_CAST")
     public override fun binaryOperationFunction(operation: String): (left: Matrix<T>, right: Matrix<T>) -> M =
         when (operation) {
@@ -67,6 +62,10 @@ public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Mat
     public operator fun T.times(m: Matrix<T>): M = m * this
 
     public companion object {
+        /**
+         * Non-boxing double matrix
+         */
+        public val real: RealMatrixContext = RealMatrixContext
 
         /**
          * A structured matrix with custom buffer
@@ -90,6 +89,11 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> :
      */
     public val elementContext: R
 
+    /**
+     * Produce a point compatible with matrix space
+     */
+    public fun point(size: Int, initializer: (Int) -> T): Point<T>
+
     public override infix fun Matrix<T>.dot(other: Matrix<T>): M {
         //TODO add typed error
         require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt
deleted file mode 100644
index 772b20f3b..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt
+++ /dev/null
@@ -1,84 +0,0 @@
-package kscience.kmath.linear
-
-import kscience.kmath.operations.RealField
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.MutableBuffer
-import kscience.kmath.structures.MutableBufferFactory
-import kscience.kmath.structures.RealBuffer
-
-@Suppress("OVERRIDE_BY_INLINE")
-public object RealMatrixContext : MatrixContext<Double, BufferMatrix<Double>> {
-
-    public override inline fun produce(
-        rows: Int,
-        columns: Int,
-        initializer: (i: Int, j: Int) -> Double,
-    ): BufferMatrix<Double> {
-        val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
-        return BufferMatrix(rows, columns, buffer)
-    }
-
-    private fun Matrix<Double>.wrap(): BufferMatrix<Double> = if (this is BufferMatrix) this else {
-        produce(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    public fun one(rows: Int, columns: Int): FeaturedMatrix<Double> = VirtualMatrix(rows, columns, DiagonalFeature) { i, j ->
-        if (i == j) 1.0 else 0.0
-    }
-
-    public override infix fun Matrix<Double>.dot(other: Matrix<Double>): BufferMatrix<Double> {
-        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
-        return produce(rowNum, other.colNum) { i, j ->
-            var res = 0.0
-            for (l in 0 until colNum) {
-                res += get(i, l) * other.get(l, j)
-            }
-            res
-        }
-    }
-
-    public override infix fun Matrix<Double>.dot(vector: Point<Double>): Point<Double> {
-        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
-        return RealBuffer(rowNum) { i ->
-            var res = 0.0
-            for (j in 0 until colNum) {
-                res += get(i, j) * vector[j]
-            }
-            res
-        }
-    }
-
-    override fun add(a: Matrix<Double>, b: Matrix<Double>): BufferMatrix<Double> {
-        require(a.rowNum == b.rowNum) { "Row number mismatch in matrix addition. Left side: ${a.rowNum}, right side: ${b.rowNum}" }
-        require(a.colNum == b.colNum) { "Column number mismatch in matrix addition. Left side: ${a.colNum}, right side: ${b.colNum}" }
-        return produce(a.rowNum, a.colNum) { i, j ->
-            a[i, j] + b[i, j]
-        }
-    }
-
-    override fun Matrix<Double>.times(value: Double): BufferMatrix<Double> =
-        produce(rowNum, colNum) { i, j -> get(i, j) * value }
-
-
-    override fun multiply(a: Matrix<Double>, k: Number): BufferMatrix<Double> =
-        produce(a.rowNum, a.colNum) { i, j -> a.get(i, j) * k.toDouble() }
-}
-
-
-/**
- * Partially optimized real-valued matrix
- */
-public val MatrixContext.Companion.real: RealMatrixContext get() = RealMatrixContext
-
-public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): FeaturedMatrix<Double> {
-    // Use existing decomposition if it is provided by matrix
-    val bufferFactory: MutableBufferFactory<Double> = MutableBuffer.Companion::real
-    val decomposition = a.getFeature() ?: lup(bufferFactory, RealField, a) { it < 1e-11 }
-    return decomposition.solveWithLUP(bufferFactory, b)
-}
-
-/**
- * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
- */
-public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): FeaturedMatrix<Double> =
-    solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum))
diff --git a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
index 0422d11b2..68a5dc262 100644
--- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
+++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
@@ -1,6 +1,11 @@
 package kscience.kmath.dimensions
 
-import kscience.kmath.linear.*
+import kscience.kmath.linear.GenericMatrixContext
+import kscience.kmath.linear.MatrixContext
+import kscience.kmath.linear.Point
+import kscience.kmath.linear.transpose
+import kscience.kmath.operations.RealField
+import kscience.kmath.operations.Ring
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.Structure2D
@@ -37,7 +42,7 @@ public interface DMatrix<T, R : Dimension, C : Dimension> : Structure2D<T> {
  * An inline wrapper for a Matrix
  */
 public inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
-    private val structure: Structure2D<T>,
+    private val structure: Structure2D<T>
 ) : DMatrix<T, R, C> {
     override val shape: IntArray get() = structure.shape
     override operator fun get(i: Int, j: Int): T = structure[i, j]
@@ -76,7 +81,7 @@ public inline class DPointWrapper<T, D : Dimension>(public val point: Point<T>)
 /**
  * Basic operations on dimension-safe matrices. Operates on [Matrix]
  */
-public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T, Matrix<T>>) {
+public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri, Matrix<T>>) {
     public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
         require(rowNum == Dimension.dim<R>().toInt()) {
             "Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum"
@@ -110,7 +115,7 @@ public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T,
     }
 
     public inline infix fun <reified R1 : Dimension, reified C1 : Dimension, reified C2 : Dimension> DMatrix<T, R1, C1>.dot(
-        other: DMatrix<T, C1, C2>,
+        other: DMatrix<T, C1, C2>
     ): DMatrix<T, R1, C2> = context { this@dot dot other }.coerce()
 
     public inline infix fun <reified R : Dimension, reified C : Dimension> DMatrix<T, R, C>.dot(vector: DPoint<T, C>): DPoint<T, R> =
@@ -134,19 +139,18 @@ public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T,
     public inline fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.transpose(): DMatrix<T, R, C> =
         context { (this@transpose as Matrix<T>).transpose() }.coerce()
 
+    /**
+     * A square unit matrix
+     */
+    public inline fun <reified D : Dimension> one(): DMatrix<T, D, D> = produce { i, j ->
+        if (i == j) context.elementContext.one else context.elementContext.zero
+    }
+
+    public inline fun <reified R : Dimension, reified C : Dimension> zero(): DMatrix<T, R, C> = produce { _, _ ->
+        context.elementContext.zero
+    }
+
     public companion object {
-        public val real: DMatrixContext<Double> = DMatrixContext(MatrixContext.real)
+        public val real: DMatrixContext<Double, RealField> = DMatrixContext(MatrixContext.real)
     }
 }
-
-
-/**
- * A square unit matrix
- */
-public inline fun <reified D : Dimension> DMatrixContext<Double>.one(): DMatrix<Double, D, D> = produce { i, j ->
-    if (i == j) 1.0 else 0.0
-}
-
-public inline fun <reified R : Dimension, reified C : Dimension> DMatrixContext<Double>.zero(): DMatrix<Double, R, C> = produce { _, _ ->
-    0.0
-}
\ No newline at end of file
diff --git a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt
index 5b330fcce..f44b16753 100644
--- a/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt
+++ b/kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt
@@ -3,7 +3,6 @@ package kscience.dimensions
 import kscience.kmath.dimensions.D2
 import kscience.kmath.dimensions.D3
 import kscience.kmath.dimensions.DMatrixContext
-import kscience.kmath.dimensions.one
 import kotlin.test.Test
 
 internal class DMatrixContextTest {
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
index 772abfbed..e8ad835e5 100644
--- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
@@ -1,7 +1,13 @@
 package kscience.kmath.real
 
-import kscience.kmath.linear.*
+import kscience.kmath.linear.FeaturedMatrix
+import kscience.kmath.linear.MatrixContext
+import kscience.kmath.linear.RealMatrixContext.elementContext
+import kscience.kmath.linear.VirtualMatrix
+import kscience.kmath.linear.inverseWithLUP
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.invoke
+import kscience.kmath.operations.sum
 import kscience.kmath.structures.Buffer
 import kscience.kmath.structures.RealBuffer
 import kscience.kmath.structures.asIterable
@@ -116,7 +122,8 @@ public fun RealMatrix.extractColumn(columnIndex: Int): RealMatrix =
     extractColumns(columnIndex..columnIndex)
 
 public fun RealMatrix.sumByColumn(): RealBuffer = RealBuffer(colNum) { j ->
-    columns[j].asIterable().sum()
+    val column = columns[j]
+    elementContext { sum(column.asIterable()) }
 }
 
 public fun RealMatrix.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->