diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1fd169b10..e542d210c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -16,13 +16,14 @@
 - ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`.
 - Coroutine-deterministic Monte-Carlo scope with a random number generator.
 - Some minor utilities to `kmath-for-real`.
+- Generic operation result parameter to `MatrixContext`
 
 ### Changed
 - Package changed from `scientifik` to `kscience.kmath`.
-- Gradle version: 6.6 -> 6.7
+- Gradle version: 6.6 -> 6.7.1
 - Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
 - `Polynomial` secondary constructor made function.
-- Kotlin version: 1.3.72 -> 1.4.20-M1
+- Kotlin version: 1.3.72 -> 1.4.20
 - `kmath-ast` doesn't depend on heavy `kotlin-reflect` library.
 - Full autodiff refactoring based on `Symbol`
 - `kmath-prob` renamed to `kmath-stat`
@@ -30,6 +31,7 @@
 - Use `Point<Double>` instead of specialized type in `kmath-for-real`
 - Optimized dot product for buffer matrices moved to `kmath-for-real`
 - EjmlMatrix context is an object
+- Matrix LUP `inverse` renamed to `inverseWithLUP`
 
 ### Deprecated
 
@@ -37,6 +39,7 @@
 - `kmath-koma` module because it doesn't support Kotlin 1.4.
 - Support of `legacy` JS backend (we will support only IR)
 - `toGrid` method.
+- Public visibility of `BufferAccessor2D`
 
 ### Fixed
 - `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
index 0ac6ea5f0..ec8714617 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
@@ -29,7 +29,7 @@ class LinearAlgebraBenchmark {
 
     @Benchmark
     fun kmathLUPInversion() {
-        MatrixContext.real.inverse(matrix)
+        MatrixContext.real.inverseWithLUP(matrix)
     }
 
     @Benchmark
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 377d167bc..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
@@ -5,8 +5,7 @@ import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.NDStructure
 import org.apache.commons.math3.linear.*
 
-public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
-    FeaturedMatrix<Double> {
+public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
     public override val rowNum: Int get() = origin.rowDimension
     public override val colNum: Int get() = origin.columnDimension
 
@@ -55,7 +54,7 @@ public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
 
 public fun RealVector.toPoint(): CMVector = CMVector(this)
 
-public object CMMatrixContext : MatrixContext<Double> {
+public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
         val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
         return CMMatrix(Array2DRowRealMatrix(array))
@@ -79,7 +78,7 @@ public object CMMatrixContext : MatrixContext<Double> {
     public override fun multiply(a: Matrix<Double>, k: Number): CMMatrix =
         CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
 
-    public override operator fun Matrix<Double>.times(value: Double): Matrix<Double> =
+    public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
         produce(rowNum, colNum) { i, j -> get(i, j) * value }
 }
 
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 f6ed34be0..8b50bbe33 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
@@ -10,7 +10,7 @@ import kscience.kmath.structures.*
 public class BufferMatrixContext<T : Any, R : Ring<T>>(
     public override val elementContext: R,
     private val bufferFactory: BufferFactory<T>,
-) : GenericMatrixContext<T, R> {
+) : GenericMatrixContext<T, R, BufferMatrix<T>> {
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix<T> {
         val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
         return BufferMatrix(rows, columns, buffer)
@@ -22,7 +22,7 @@ public class BufferMatrixContext<T : Any, R : Ring<T>>(
 }
 
 @Suppress("OVERRIDE_BY_INLINE")
-public object RealMatrixContext : GenericMatrixContext<Double, RealField> {
+public object RealMatrixContext : GenericMatrixContext<Double, RealField, BufferMatrix<Double>> {
     public override val elementContext: RealField
         get() = RealField
 
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
index 350379717..68272203c 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
@@ -57,7 +57,7 @@ public inline fun <reified T : Any> Matrix<*>.getFeature(): T? =
 /**
  * Diagonal matrix of ones. The matrix is virtual no actual matrix is created
  */
-public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
+public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
     VirtualMatrix(rows, columns, DiagonalFeature) { i, j ->
         if (i == j) elementContext.one else elementContext.zero
     }
@@ -66,7 +66,7 @@ public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, colu
 /**
  * A virtual matrix of zeroes
  */
-public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
+public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
     VirtualMatrix(rows, columns) { _, _ -> elementContext.zero }
 
 public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature
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 bb80bcafd..099fa1909 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
@@ -1,25 +1,18 @@
 package kscience.kmath.linear
 
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.Ring
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.BufferAccessor2D
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.Structure2D
-import kotlin.reflect.KClass
+import kscience.kmath.operations.*
+import kscience.kmath.structures.*
 
 /**
  * Common implementation of [LUPDecompositionFeature]
  */
 public class LUPDecomposition<T : Any>(
-    public val context: GenericMatrixContext<T, out Field<T>>,
+    public val context: MatrixContext<T, FeaturedMatrix<T>>,
+    public val elementContext: Field<T>,
     public val lu: Structure2D<T>,
     public val pivot: IntArray,
-    private val even: Boolean
+    private val even: Boolean,
 ) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
-    public val elementContext: Field<T>
-        get() = context.elementContext
 
     /**
      * Returns the matrix L of the decomposition.
@@ -64,23 +57,25 @@ public class LUPDecomposition<T : Any>(
 
 }
 
-public fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T): T =
+@PublishedApi
+internal fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, *>.abs(value: T): T =
     if (value > elementContext.zero) value else elementContext { -value }
 
 /**
- * Create a lup decomposition of generic matrix
+ * Create a lup decomposition of generic matrix.
  */
-public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
-    type: KClass<T>,
+public fun <T : Comparable<T>> MatrixContext<T, FeaturedMatrix<T>>.lup(
+    factory: MutableBufferFactory<T>,
+    elementContext: Field<T>,
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
+    checkSingular: (T) -> Boolean,
 ): LUPDecomposition<T> {
     require(matrix.rowNum == matrix.colNum) { "LU decomposition supports only square matrices" }
     val m = matrix.colNum
     val pivot = IntArray(matrix.rowNum)
 
     //TODO just waits for KEEP-176
-    BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
+    BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
         elementContext {
             val lu = create(matrix)
 
@@ -112,14 +107,14 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
                     luRow[col] = sum
 
                     // maintain best permutation choice
-                    if (this@lup.abs(sum) > largest) {
-                        largest = this@lup.abs(sum)
+                    if (abs(sum) > largest) {
+                        largest = abs(sum)
                         max = row
                     }
                 }
 
                 // Singularity check
-                check(!checkSingular(this@lup.abs(lu[max, col]))) { "The matrix is singular" }
+                check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
 
                 // Pivot if necessary
                 if (max != col) {
@@ -143,23 +138,23 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
                 for (row in col + 1 until m) lu[row, col] /= luDiag
             }
 
-            return LUPDecomposition(this@lup, lu.collect(), pivot, even)
+            return LUPDecomposition(this@lup, elementContext, lu.collect(), pivot, even)
         }
     }
 }
 
-public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.lup(
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): LUPDecomposition<T> = lup(T::class, matrix, checkSingular)
+    noinline checkSingular: (T) -> Boolean,
+): LUPDecomposition<T> = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular)
 
-public fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
-    lup(Double::class, matrix) { it < 1e-11 }
+public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
+    lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
 
-public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T>): Matrix<T> {
+public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> {
     require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
 
-    BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
+    BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
         elementContext {
             // Apply permutations to b
             val bp = create { _, _ -> zero }
@@ -201,27 +196,34 @@ public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T
     }
 }
 
-public inline fun <reified T : Any> LUPDecomposition<T>.solve(matrix: Matrix<T>): Matrix<T> = solve(T::class, matrix)
+public inline fun <reified T : Any> LUPDecomposition<T>.solveWithLUP(matrix: Matrix<T>): Matrix<T> =
+    solveWithLUP(MutableBuffer.Companion::auto, matrix)
 
 /**
- * Solve a linear equation **a*x = b**
+ * Solve a linear equation **a*x = b** using LUP decomposition
  */
-public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.solveWithLUP(
     a: Matrix<T>,
     b: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): Matrix<T> {
+    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+    noinline checkSingular: (T) -> Boolean,
+): FeaturedMatrix<T> {
     // Use existing decomposition if it is provided by matrix
-    val decomposition = a.getFeature() ?: lup(T::class, a, checkSingular)
-    return decomposition.solve(T::class, b)
+    val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular)
+    return decomposition.solveWithLUP(bufferFactory, b)
 }
 
-public fun RealMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b) { it < 1e-11 }
+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>.inverse(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.inverseWithLUP(
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
+    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+    noinline checkSingular: (T) -> Boolean,
+): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
 
-public fun RealMatrixContext.inverse(matrix: Matrix<Double>): Matrix<Double> =
-    solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 }
+/**
+ * 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 f4dbce89a..d9dc57b0f 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
@@ -12,15 +12,16 @@ import kscience.kmath.structures.asSequence
 /**
  * Basic operations on matrices. Operates on [Matrix]
  */
-public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
+public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Matrix<T>> {
     /**
      * Produce a matrix with this context and given dimensions
      */
-    public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
+    public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
 
-    public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): Matrix<T> = when (operation) {
+    @Suppress("UNCHECKED_CAST")
+    public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): M = when (operation) {
         "dot" -> left dot right
-        else -> super.binaryOperation(operation, left, right)
+        else -> super.binaryOperation(operation, left, right) as M
     }
 
     /**
@@ -30,7 +31,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param other the multiplier.
      * @return the dot product.
      */
-    public infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
+    public infix fun Matrix<T>.dot(other: Matrix<T>): M
 
     /**
      * Computes the dot product of this matrix and a vector.
@@ -48,7 +49,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param value the multiplier.
      * @receiver the product.
      */
-    public operator fun Matrix<T>.times(value: T): Matrix<T>
+    public operator fun Matrix<T>.times(value: T): M
 
     /**
      * Multiplies an element by a matrix of it.
@@ -57,7 +58,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param value the multiplier.
      * @receiver the product.
      */
-    public operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
+    public operator fun T.times(m: Matrix<T>): M = m * this
 
     public companion object {
         /**
@@ -70,18 +71,18 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
          */
         public fun <T : Any, R : Ring<T>> buffered(
             ring: R,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
-        ): GenericMatrixContext<T, R> = BufferMatrixContext(ring, bufferFactory)
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+        ): GenericMatrixContext<T, R, BufferMatrix<T>> = BufferMatrixContext(ring, bufferFactory)
 
         /**
          * Automatic buffered matrix, unboxed if it is possible
          */
-        public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
+        public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R, BufferMatrix<T>> =
             buffered(ring, Buffer.Companion::auto)
     }
 }
 
-public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
+public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> : MatrixContext<T, M> {
     /**
      * The ring context for matrix elements
      */
@@ -92,7 +93,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
      */
     public fun point(size: Int, initializer: (Int) -> T): Point<T>
 
-    public override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<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})" }
 
@@ -113,10 +114,10 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         }
     }
 
-    public override operator fun Matrix<T>.unaryMinus(): Matrix<T> =
+    public override operator fun Matrix<T>.unaryMinus(): M =
         produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } }
 
-    public override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
+    public override fun add(a: Matrix<T>, b: Matrix<T>): M {
         require(a.rowNum == b.rowNum && a.colNum == b.colNum) {
             "Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]"
         }
@@ -124,7 +125,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } }
     }
 
-    public override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
+    public override operator fun Matrix<T>.minus(b: Matrix<T>): M {
         require(rowNum == b.rowNum && colNum == b.colNum) {
             "Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]"
         }
@@ -132,11 +133,11 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } }
     }
 
-    public override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
+    public override fun multiply(a: Matrix<T>, k: Number): M =
         produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } }
 
-    public operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
+    public operator fun Number.times(matrix: FeaturedMatrix<T>): M = multiply(matrix, this)
 
-    public override operator fun Matrix<T>.times(value: T): Matrix<T> =
+    public override operator fun Matrix<T>.times(value: T): M =
         produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } }
 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
index 657da6b4e..4527a2a42 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
@@ -38,6 +38,11 @@ public fun <T> Space<T>.average(data: Iterable<T>): T = sum(data) / data.count()
  */
 public fun <T> Space<T>.average(data: Sequence<T>): T = sum(data) / data.count()
 
+/**
+ * Absolute of the comparable [value]
+ */
+public fun <T : Comparable<T>> Space<T>.abs(value: T): T = if (value > zero) value else -value
+
 /**
  * Returns the sum of all elements in the iterable in provided space.
  *
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
index 00fc4aef0..5d7ba611f 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
@@ -1,25 +1,31 @@
 package kscience.kmath.structures
 
-import kotlin.reflect.KClass
-
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
  */
-public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val rowNum: Int, public val colNum: Int) {
+internal class BufferAccessor2D<T : Any>(
+    public val rowNum: Int,
+    public val colNum: Int,
+    val factory: MutableBufferFactory<T>,
+) {
     public operator fun Buffer<T>.get(i: Int, j: Int): T = get(i + colNum * j)
 
     public operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
         set(i + colNum * j, value)
     }
 
-    public inline fun create(init: (i: Int, j: Int) -> T): MutableBuffer<T> =
-        MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
+    public inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer<T> =
+        factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
 
     public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
 
     //TODO optimize wrapper
-    public fun MutableBuffer<T>.collect(): Structure2D<T> =
-        NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D()
+    public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.build(
+        DefaultStrides(intArrayOf(rowNum, colNum)),
+        factory
+    ) { (i, j) ->
+        get(i, j)
+    }.as2D()
 
     public inner class Row(public val buffer: MutableBuffer<T>, public val rowIndex: Int) : MutableBuffer<T> {
         override val size: Int get() = colNum
@@ -30,7 +36,7 @@ public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val ro
             buffer[rowIndex, index] = value
         }
 
-        override fun copy(): MutableBuffer<T> = MutableBuffer.auto(type, colNum) { get(it) }
+        override fun copy(): MutableBuffer<T> = factory(colNum) { get(it) }
         override operator fun iterator(): Iterator<T> = (0 until colNum).map(::get).iterator()
 
     }
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
index de07d3639..28dfe46ec 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
@@ -9,7 +9,7 @@ class RealLUSolverTest {
     @Test
     fun testInvertOne() {
         val matrix = MatrixContext.real.one(2, 2)
-        val inverted = MatrixContext.real.inverse(matrix)
+        val inverted = MatrixContext.real.inverseWithLUP(matrix)
         assertEquals(matrix, inverted)
     }
 
@@ -37,7 +37,7 @@ class RealLUSolverTest {
             1.0, 3.0
         )
 
-        val inverted = MatrixContext.real.inverse(matrix)
+        val inverted = MatrixContext.real.inverseWithLUP(matrix)
 
         val expected = Matrix.square(
             0.375, -0.125,
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 4111cb78d..68a5dc262 100644
--- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
+++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
@@ -42,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>(
-    public 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]
@@ -81,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, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri>) {
+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"
diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
index 1ed57fcb3..31792e39c 100644
--- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
+++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
@@ -1,11 +1,9 @@
 package kscience.kmath.ejml
 
-import org.ejml.simple.SimpleMatrix
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.Point
-import kscience.kmath.operations.Space
-import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
+import org.ejml.simple.SimpleMatrix
 
 /**
  * Converts this matrix to EJML one.
@@ -18,7 +16,7 @@ public fun Matrix<Double>.toEjml(): EjmlMatrix =
  *
  * @author Iaroslav Postovalov
  */
-public object EjmlMatrixContext : MatrixContext<Double> {
+public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
 
     /**
      * Converts this vector to EJML one.
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 c0657c0bf..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,12 +1,16 @@
 package kscience.kmath.real
 
+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.*
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.RealBuffer
+import kscience.kmath.structures.asIterable
 import kotlin.math.pow
 
 /*
@@ -21,7 +25,7 @@ import kotlin.math.pow
  *  Functions that help create a real (Double) matrix
  */
 
-public typealias RealMatrix = Matrix<Double>
+public typealias RealMatrix = FeaturedMatrix<Double>
 
 public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum, initializer)
@@ -148,6 +152,11 @@ public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix =
         transform(get(i, j))
     }
 
+/**
+ * Inverse a square real matrix using LUP decomposition
+ */
+public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this)
+
 //extended operations
 
 public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }