diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt
index 8816ba5cf..c2d74379e 100644
--- a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt
+++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt
@@ -1,10 +1,13 @@
 package scientifik.kmath.linear
 
 import koma.matrix.ejml.EJMLMatrixFactory
+import scientifik.kmath.operations.RealField
 import scientifik.kmath.structures.Matrix
+import kotlin.contracts.ExperimentalContracts
 import kotlin.random.Random
 import kotlin.system.measureTimeMillis
 
+@ExperimentalContracts
 fun main() {
     val random = Random(12224)
     val dim = 100
@@ -32,10 +35,8 @@ fun main() {
 
     //commons-math
 
-    val cmContext = CMLUPSolver
-
     val commonsTime = measureTimeMillis {
-        cmContext.run {
+        CMMatrixContext.run {
             val cm = matrix.toCM()             //avoid overhead on conversion
             repeat(n) {
                 val res = inverse(cm)
@@ -48,10 +49,8 @@ fun main() {
 
     //koma-ejml
 
-    val komaContext = KomaMatrixContext(EJMLMatrixFactory())
-
     val komaTime = measureTimeMillis {
-        komaContext.run {
+        KomaMatrixContext(EJMLMatrixFactory(), RealField).run {
             val km = matrix.toKoma()      //avoid overhead on conversion
             repeat(n) {
                 val res = inverse(km)
diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt
index 7b1bd9e7e..31a890867 100644
--- a/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt
+++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt
@@ -1,6 +1,7 @@
 package scientifik.kmath.linear
 
 import koma.matrix.ejml.EJMLMatrixFactory
+import scientifik.kmath.operations.RealField
 import scientifik.kmath.structures.Matrix
 import kotlin.random.Random
 import kotlin.system.measureTimeMillis
@@ -27,7 +28,7 @@ fun main() {
     }
 
 
-    KomaMatrixContext(EJMLMatrixFactory()).run {
+    KomaMatrixContext(EJMLMatrixFactory(), RealField).run {
         val komaMatrix1 = matrix1.toKoma()
         val komaMatrix2 = matrix2.toKoma()
 
diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt
index 46b1b56dd..111b76ff0 100644
--- a/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt
+++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt
@@ -23,7 +23,7 @@ fun main() {
 
     val complexTime = measureTimeMillis {
         complexField.run {
-            var res: ComplexNDElement = one
+            var res = one
             repeat(n) {
                 res += 1.0
             }
diff --git a/build.gradle.kts b/build.gradle.kts
index 964eb4998..3f25a622c 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -1,9 +1,6 @@
 import com.moowork.gradle.node.NodeExtension
 import com.moowork.gradle.node.npm.NpmTask
 import com.moowork.gradle.node.task.NodeTask
-import org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension
-import org.jetbrains.kotlin.gradle.tasks.Kotlin2JsCompile
-import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
 
 buildscript {
     val kotlinVersion: String by rootProject.extra("1.3.21")
@@ -194,6 +191,8 @@ subprojects {
                 sourceSets.all {
                     languageSettings.progressiveMode = true
                     languageSettings.enableLanguageFeature("InlineClasses")
+                    languageSettings.useExperimentalAnnotation("ExperimentalContracts")
+                    //languageSettings.enableLanguageFeature("Contracts")
                 }
             }
 
diff --git a/doc/linear.md b/doc/linear.md
index 6208deaff..bbcc435ba 100644
--- a/doc/linear.md
+++ b/doc/linear.md
@@ -1 +1,19 @@
-**TODO**
\ No newline at end of file
+## Basic linear algebra layout
+
+Kmath support for linear algebra organized in a context-oriented way. Meaning that operations are in most cases declared 
+in context classes, and are not the members of classes that store data. This allows more flexible approach to maintain multiple 
+back-ends. The new operations added as extensions to contexts instead of being member functions of data structures.
+
+Two major contexts used for linear algebra and hyper-geometry:
+
+* `VectorSpace` forms a mathematical space on top of array-like structure (`Buffer` and its typealias `Point` used for geometry). 
+
+* `MatrixContext` forms a space-like context for 2d-structures. It does not store matrix size and therefore does not implement
+`Space` interface (it is not possible to create zero element without knowing the matrix size). 
+
+## Vector spaces
+
+
+## Matrix operations
+
+## Back-end overview
\ No newline at end of file
diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
index 513dfe517..9103f0fb2 100644
--- a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
+++ b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
@@ -27,7 +27,7 @@ fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {
     CMMatrix(Array2DRowRealMatrix(array))
 }
 
-fun RealMatrix.toMatrix() = CMMatrix(this)
+fun RealMatrix.asMatrix() = CMMatrix(this)
 
 class CMVector(val origin: RealVector) : Point<Double> {
     override val size: Int get() = origin.dimension
@@ -47,7 +47,6 @@ fun Point<Double>.toCM(): CMVector = if (this is CMVector) {
 fun RealVector.toPoint() = CMVector(this)
 
 object CMMatrixContext : MatrixContext<Double> {
-
     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))
@@ -59,35 +58,19 @@ object CMMatrixContext : MatrixContext<Double> {
     override fun Matrix<Double>.dot(vector: Point<Double>): CMVector =
         CMVector(this.toCM().origin.preMultiply(vector.toCM().origin))
 
-
     override fun Matrix<Double>.unaryMinus(): CMMatrix =
         produce(rowNum, colNum) { i, j -> -get(i, j) }
 
-    override fun Matrix<Double>.plus(b: Matrix<Double>) =
-        CMMatrix(this.toCM().origin.multiply(b.toCM().origin))
+    override fun add(a: Matrix<Double>, b: Matrix<Double>) =
+        CMMatrix(a.toCM().origin.multiply(b.toCM().origin))
 
     override fun Matrix<Double>.minus(b: Matrix<Double>) =
         CMMatrix(this.toCM().origin.subtract(b.toCM().origin))
 
-    override fun Matrix<Double>.times(value: Double) =
-        CMMatrix(this.toCM().origin.scalarMultiply(value.toDouble()))
-}
+    override fun multiply(a: Matrix<Double>, k: Number) =
+        CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
 
-object CMLUPSolver: LinearSolver<Double>{
-    override fun solve(a: Matrix<Double>, b: Matrix<Double>): CMMatrix {
-        val decomposition = LUDecomposition(a.toCM().origin)
-        return decomposition.solver.solve(b.toCM().origin).toMatrix()
-    }
-
-    override fun solve(a: Matrix<Double>, b: Point<Double>): CMVector {
-        val decomposition = LUDecomposition(a.toCM().origin)
-        return decomposition.solver.solve(b.toCM().origin).toPoint()
-    }
-
-    override fun inverse(a: Matrix<Double>): CMMatrix {
-        val decomposition = LUDecomposition(a.toCM().origin)
-        return decomposition.solver.inverse.toMatrix()
-    }
+    override fun Matrix<Double>.times(value: Double): Matrix<Double>  = produce(rowNum,colNum){i,j-> get(i,j)*value}
 }
 
 operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin))
diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMSolver.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMSolver.kt
new file mode 100644
index 000000000..a56e8c93e
--- /dev/null
+++ b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMSolver.kt
@@ -0,0 +1,39 @@
+package scientifik.kmath.linear
+
+import org.apache.commons.math3.linear.*
+import scientifik.kmath.structures.Matrix
+
+enum class CMDecomposition {
+    LUP,
+    QR,
+    RRQR,
+    EIGEN,
+    CHOLESKY
+}
+
+
+fun CMMatrixContext.solver(a: Matrix<Double>, decomposition: CMDecomposition = CMDecomposition.LUP) =
+    when (decomposition) {
+        CMDecomposition.LUP -> LUDecomposition(a.toCM().origin).solver
+        CMDecomposition.RRQR -> RRQRDecomposition(a.toCM().origin).solver
+        CMDecomposition.QR -> QRDecomposition(a.toCM().origin).solver
+        CMDecomposition.EIGEN -> EigenDecomposition(a.toCM().origin).solver
+        CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM().origin).solver
+    }
+
+fun CMMatrixContext.solve(
+    a: Matrix<Double>,
+    b: Matrix<Double>,
+    decomposition: CMDecomposition = CMDecomposition.LUP
+) = solver(a, decomposition).solve(b.toCM().origin).asMatrix()
+
+fun CMMatrixContext.solve(
+    a: Matrix<Double>,
+    b: Point<Double>,
+    decomposition: CMDecomposition = CMDecomposition.LUP
+) = solver(a, decomposition).solve(b.toCM().origin).toPoint()
+
+fun CMMatrixContext.inverse(
+    a: Matrix<Double>,
+    decomposition: CMDecomposition = CMDecomposition.LUP
+) = solver(a, decomposition).inverse.asMatrix()
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
index fc16c1f8e..739c170d4 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
@@ -1,6 +1,5 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Ring
 import scientifik.kmath.structures.*
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt
index 62c360f26..8ac18ee4a 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt
@@ -7,109 +7,6 @@ import scientifik.kmath.structures.*
 import scientifik.kmath.structures.Buffer.Companion.boxing
 import kotlin.math.sqrt
 
-
-/**
- * Basic operations on matrices. Operates on [Matrix]
- */
-interface MatrixContext<T : Any> {
-    /**
-     * Produce a matrix with this context and given dimensions
-     */
-    fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
-
-    infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
-
-    infix fun Matrix<T>.dot(vector: Point<T>): Point<T>
-
-    operator fun Matrix<T>.unaryMinus(): Matrix<T>
-
-    operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T>
-
-    operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T>
-
-    operator fun Matrix<T>.times(value: T): Matrix<T>
-
-    operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
-
-    companion object {
-        /**
-         * Non-boxing double matrix
-         */
-        val real = BufferMatrixContext(RealField, Buffer.Companion::auto)
-
-        /**
-         * A structured matrix with custom buffer
-         */
-        fun <T : Any, R : Ring<T>> buffered(
-            ring: R,
-            bufferFactory: BufferFactory<T> = ::boxing
-        ): GenericMatrixContext<T, R> =
-            BufferMatrixContext(ring, bufferFactory)
-
-        /**
-         * Automatic buffered matrix, unboxed if it is possible
-         */
-        inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
-            buffered(ring, Buffer.Companion::auto)
-    }
-}
-
-interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
-    /**
-     * The ring context for matrix elements
-     */
-    val elementContext: R
-
-    /**
-     * Produce a point compatible with matrix space
-     */
-    fun point(size: Int, initializer: (Int) -> T): Point<T>
-
-    override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
-        //TODO add typed error
-        if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})")
-        return produce(rowNum, other.colNum) { i, j ->
-            val row = rows[i]
-            val column = other.columns[j]
-            with(elementContext) {
-                sum(row.asSequence().zip(column.asSequence(), ::multiply))
-            }
-        }
-    }
-
-    override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
-        //TODO add typed error
-        if (this.colNum != vector.size) error("Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})")
-        return point(rowNum) { i ->
-            val row = rows[i]
-            with(elementContext) {
-                sum(row.asSequence().zip(vector.asSequence(), ::multiply))
-            }
-        }
-    }
-
-    override operator fun Matrix<T>.unaryMinus() =
-        produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } }
-
-    override operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T> {
-        if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] + [${b.rowNum},${b.colNum}]")
-        return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
-    }
-
-    override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
-        if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]")
-        return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
-    }
-
-    operator fun Matrix<T>.times(number: Number): Matrix<T> =
-        produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * number } }
-
-    operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
-
-    override fun Matrix<T>.times(value: T): Matrix<T> =
-        produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
-}
-
 /**
  * A 2d structure plus optional matrix-specific features
  */
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
index dc46216e8..4bc463333 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
@@ -4,6 +4,9 @@ import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Ring
 import scientifik.kmath.structures.*
+import kotlin.contracts.ExperimentalContracts
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
 import kotlin.reflect.KClass
 
 /**
@@ -63,7 +66,12 @@ class LUPDecomposition<T : Any>(
 
 }
 
-open class BufferAccessor<T : Any>(val type: KClass<T>, val field: Field<T>, val rowNum: Int, val colNum: Int) {
+internal open class BufferAccessor<T : Any>(
+    val type: KClass<T>,
+    val field: Field<T>,
+    val rowNum: Int,
+    val colNum: Int
+) {
     open operator fun MutableBuffer<T>.get(i: Int, j: Int) = get(i + colNum * j)
     open operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
         set(i + colNum * j, value)
@@ -102,7 +110,8 @@ open class BufferAccessor<T : Any>(val type: KClass<T>, val field: Field<T>, val
 /**
  * Specialized LU operations for Doubles
  */
-class RealBufferAccessor(rowNum: Int, colNum: Int) : BufferAccessor<Double>(Double::class, RealField, rowNum, colNum) {
+private class RealBufferAccessor(rowNum: Int, colNum: Int) :
+    BufferAccessor<Double>(Double::class, RealField, rowNum, colNum) {
     override fun MutableBuffer<Double>.get(i: Int, j: Int) = (this as DoubleBuffer).array[i + colNum * j]
     override fun MutableBuffer<Double>.set(i: Int, j: Int, value: Double) {
         (this as DoubleBuffer).array[i + colNum * j] = value
@@ -125,24 +134,33 @@ class RealBufferAccessor(rowNum: Int, colNum: Int) : BufferAccessor<Double>(Doub
     }
 }
 
-fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.buildAccessor(
-    type:KClass<T>,
+@ExperimentalContracts
+private inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.withAccessor(
+    type: KClass<T>,
     rowNum: Int,
-    colNum: Int
-): BufferAccessor<T> {
-    return if (elementContext == RealField) {
+    colNum: Int,
+    block: BufferAccessor<T>.() -> Unit
+) {
+    contract {
+        callsInPlace(block, InvocationKind.EXACTLY_ONCE)
+    }
+    if (elementContext == RealField) {
         @Suppress("UNCHECKED_CAST")
         RealBufferAccessor(rowNum, colNum) as BufferAccessor<T>
     } else {
         BufferAccessor(type, elementContext, rowNum, colNum)
-    }
+    }.run(block)
 }
 
-fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T) =
+private fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T) =
     if (value > elementContext.zero) value else with(elementContext) { -value }
 
 
-fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
+/**
+ * Create a lup decomposition of generic matrix
+ */
+@ExperimentalContracts
+fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
     type: KClass<T>,
     matrix: Matrix<T>,
     checkSingular: (T) -> Boolean
@@ -155,7 +173,7 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
     val m = matrix.colNum
     val pivot = IntArray(matrix.rowNum)
 
-    buildAccessor(type, matrix.rowNum, matrix.colNum).run {
+    withAccessor(type, matrix.rowNum, matrix.colNum) {
 
         val lu = create(matrix)
 
@@ -170,21 +188,12 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
 
             // upper
             for (row in 0 until col) {
-//                var sum = lu[row, col]
-//                for (i in 0 until row) {
-//                    sum -= lu[row, i] * lu[i, col]
-//                }
                 val sum = lu.innerProduct(row, col, row)
                 lu[row, col] = field.run { lu[row, col] - sum }
             }
 
             // lower
             val max = (col until m).maxBy { row ->
-                //                var sum = lu[row, col]
-//                for (i in 0 until col) {
-//                    sum -= lu[row, i] * lu[i, col]
-//                }
-//                lu[row, col] = sum
                 val sum = lu.innerProduct(row, col, col)
                 lu[row, col] = field.run { lu[row, col] - sum }
                 abs(sum)
@@ -214,14 +223,17 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
                 //lu[row, col] = lu[row, col] / luDiag
             }
         }
-        return scientifik.kmath.linear.LUPDecomposition(elementContext, lu.collect(), pivot, even)
-
+        return LUPDecomposition(elementContext, lu.collect(), pivot, even)
     }
 }
 
+@ExperimentalContracts
+fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>) = lup(Double::class, matrix) { it < 1e-11 }
+
 /**
  * Solve a linear equation **a*x = b**
  */
+@ExperimentalContracts
 fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
     type: KClass<T>,
     a: Matrix<T>,
@@ -233,9 +245,9 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
     }
 
     // Use existing decomposition if it is provided by matrix
-    val decomposition = a.getFeature() ?: lupDecompose(type, a, checkSingular)
+    val decomposition = a.getFeature() ?: lup(type, a, checkSingular)
 
-    buildAccessor(type, a.rowNum, a.colNum).run {
+    withAccessor(type, a.rowNum, a.colNum) {
 
         val lu = create(decomposition.lu)
 
@@ -271,14 +283,19 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
 
         return produce(a.rowNum, a.colNum) { i, j -> bp[i, j] }
     }
-
 }
 
+@ExperimentalContracts
+fun GenericMatrixContext<Double, RealField>.solve(a: Matrix<Double>, b: Matrix<Double>) =
+    solve(Double::class, a, b) { it < 1e-11 }
+
+@ExperimentalContracts
 inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.inverse(
     matrix: Matrix<T>,
     noinline checkSingular: (T) -> Boolean
 ) =
     solve(T::class, matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
 
+@ExperimentalContracts
 fun GenericMatrixContext<Double, RealField>.inverse(matrix: Matrix<Double>) =
     inverse(matrix) { it < 1e-11 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 4254fd7ea..7e631f5fc 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -13,7 +13,7 @@ import scientifik.kmath.structures.asSequence
  */
 interface LinearSolver<T : Any> {
     fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
-    fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.toMatrix()).asPoint()
+    fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.asMatrix()).asPoint()
     fun inverse(a: Matrix<T>): Matrix<T>
 }
 
@@ -43,4 +43,4 @@ fun <T : Any> Matrix<T>.asPoint(): Point<T> =
         error("Can't convert matrix with more than one column to vector")
     }
 
-fun <T : Any> Point<T>.toMatrix() = VirtualMatrix(size, 1) { i, _ -> get(i) }
\ No newline at end of file
+fun <T : Any> Point<T>.asMatrix() = VirtualMatrix(size, 1) { i, _ -> get(i) }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt
new file mode 100644
index 000000000..f97c55cc4
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt
@@ -0,0 +1,106 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.operations.SpaceOperations
+import scientifik.kmath.operations.sum
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.BufferFactory
+import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.asSequence
+
+/**
+ * Basic operations on matrices. Operates on [Matrix]
+ */
+interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
+    /**
+     * Produce a matrix with this context and given dimensions
+     */
+    fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
+
+    infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
+
+    infix fun Matrix<T>.dot(vector: Point<T>): Point<T>
+
+    operator fun Matrix<T>.times(value: T): Matrix<T>
+
+    operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
+
+    companion object {
+        /**
+         * Non-boxing double matrix
+         */
+        val real = BufferMatrixContext(RealField, Buffer.Companion::auto)
+
+        /**
+         * A structured matrix with custom buffer
+         */
+        fun <T : Any, R : Ring<T>> buffered(
+            ring: R,
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
+        ): GenericMatrixContext<T, R> =
+            BufferMatrixContext(ring, bufferFactory)
+
+        /**
+         * Automatic buffered matrix, unboxed if it is possible
+         */
+        inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
+            buffered(ring, Buffer.Companion::auto)
+    }
+}
+
+interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
+    /**
+     * The ring context for matrix elements
+     */
+    val elementContext: R
+
+    /**
+     * Produce a point compatible with matrix space
+     */
+    fun point(size: Int, initializer: (Int) -> T): Point<T>
+
+    override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
+        //TODO add typed error
+        if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})")
+        return produce(rowNum, other.colNum) { i, j ->
+            val row = rows[i]
+            val column = other.columns[j]
+            with(elementContext) {
+                sum(row.asSequence().zip(column.asSequence(), ::multiply))
+            }
+        }
+    }
+
+    override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
+        //TODO add typed error
+        if (this.colNum != vector.size) error("Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})")
+        return point(rowNum) { i ->
+            val row = rows[i]
+            with(elementContext) {
+                sum(row.asSequence().zip(vector.asSequence(), ::multiply))
+            }
+        }
+    }
+
+    override operator fun Matrix<T>.unaryMinus() =
+        produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } }
+
+    override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
+        if (a.rowNum != b.rowNum || a.colNum != b.colNum) error("Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]")
+        return produce(a.rowNum, a.colNum) { i, j -> elementContext.run { a.get(i, j) + b[i, j] } }
+    }
+
+    override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
+        if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]")
+        return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
+    }
+
+    override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
+        produce(a.rowNum, a.colNum) { i, j -> elementContext.run { a.get(i, j) * k } }
+
+    operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
+
+    override fun Matrix<T>.times(value: T): Matrix<T> =
+        produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
+}
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index d74d0ed4f..18a0021c3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -4,68 +4,23 @@ import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Space
 import scientifik.kmath.operations.SpaceElement
 import scientifik.kmath.structures.Buffer
-import scientifik.kmath.structures.BufferFactory
 import scientifik.kmath.structures.asSequence
+import kotlin.jvm.JvmName
 
 typealias Point<T> = Buffer<T>
 
-/**
- * A linear space for vectors.
- * Could be used on any point-like structure
- */
-interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
 
-    val size: Int
-
-    val space: S
-
-    fun produce(initializer: (Int) -> T): Point<T>
-
-    /**
-     * Produce a space-element of this vector space for expressions
-     */
-    fun produceElement(initializer: (Int) -> T): Vector<T, S>
-
-    override val zero: Point<T> get() = produce { space.zero }
-
-    override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { with(space) { a[it] + b[it] } }
-
-    override fun multiply(a: Point<T>, k: Number): Point<T> = produce { with(space) { a[it] * k } }
-
-    //TODO add basis
-
-    companion object {
-
-        private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
-
-        /**
-         * Non-boxing double vector space
-         */
-        fun real(size: Int): BufferVectorSpace<Double, RealField> {
-            return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, Buffer.Companion::auto) }
-        }
-
-        /**
-         * A structured vector space with custom buffer
-         */
-        fun <T : Any, S : Space<T>> buffered(
-            size: Int,
-            space: S,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
-        ): VectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
-
-        /**
-         * Automatic buffered vector, unboxed if it is possible
-         */
-        inline fun <reified T : Any, S : Space<T>> smart(size: Int, space: S): VectorSpace<T, S> =
-            buffered(size, space, Buffer.Companion::auto)
-    }
-}
+fun <T : Any, S : Space<T>> BufferVectorSpace<T, S>.produceElement(initializer: (Int) -> T): Vector<T, S> =
+    BufferVector(this, produce(initializer))
 
+@JvmName("produceRealElement")
+fun BufferVectorSpace<Double, RealField>.produceElement(initializer: (Int) -> Double): Vector<Double, RealField> =
+    BufferVector(this, produce(initializer))
 
 /**
  * A point coupled to the linear space
  */
+@Deprecated("Use VectorContext instead")
 interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, Vector<T, S>, VectorSpace<T, S>>, Point<T> {
     override val size: Int get() = context.size
 
@@ -90,16 +45,7 @@ interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, Vector<T, S>, V
     }
 }
 
-data class BufferVectorSpace<T : Any, S : Space<T>>(
-    override val size: Int,
-    override val space: S,
-    val bufferFactory: BufferFactory<T>
-) : VectorSpace<T, S> {
-    override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer)
-    override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
-}
-
-
+@Deprecated("Use VectorContext instead")
 data class BufferVector<T : Any, S : Space<T>>(override val context: VectorSpace<T, S>, val buffer: Buffer<T>) :
     Vector<T, S> {
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VectorSpace.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VectorSpace.kt
new file mode 100644
index 000000000..8e14e2882
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VectorSpace.kt
@@ -0,0 +1,75 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.operations.Space
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.BufferFactory
+
+/**
+ * A linear space for vectors.
+ * Could be used on any point-like structure
+ */
+interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
+
+    val size: Int
+
+    val space: S
+
+    fun produce(initializer: (Int) -> T): Point<T>
+
+    /**
+     * Produce a space-element of this vector space for expressions
+     */
+    //fun produceElement(initializer: (Int) -> T): Vector<T, S>
+
+    override val zero: Point<T> get() = produce { space.zero }
+
+    override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { with(space) { a[it] + b[it] } }
+
+    override fun multiply(a: Point<T>, k: Number): Point<T> = produce { with(space) { a[it] * k } }
+
+    //TODO add basis
+
+    companion object {
+
+        private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
+
+        /**
+         * Non-boxing double vector space
+         */
+        fun real(size: Int): BufferVectorSpace<Double, RealField> {
+            return realSpaceCache.getOrPut(size) {
+                BufferVectorSpace(
+                    size,
+                    RealField,
+                    Buffer.Companion::auto
+                )
+            }
+        }
+
+        /**
+         * A structured vector space with custom buffer
+         */
+        fun <T : Any, S : Space<T>> buffered(
+            size: Int,
+            space: S,
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
+        ) = BufferVectorSpace(size, space, bufferFactory)
+
+        /**
+         * Automatic buffered vector, unboxed if it is possible
+         */
+        inline fun <reified T : Any, S : Space<T>> auto(size: Int, space: S): VectorSpace<T, S> =
+            buffered(size, space, Buffer.Companion::auto)
+    }
+}
+
+
+class BufferVectorSpace<T : Any, S : Space<T>>(
+    override val size: Int,
+    override val space: S,
+    val bufferFactory: BufferFactory<T>
+) : VectorSpace<T, S> {
+    override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer)
+    //override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
index 887a480bb..c6f0f7a48 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
@@ -73,6 +73,7 @@ interface NDSpace<T, S : Space<T>, N : NDStructure<T>> : Space<N>, NDAlgebra<T,
      */
     override fun multiply(a: N, k: Number): N = map(a) { multiply(it, k) }
 
+    //TODO move to extensions after KEEP-176
     operator fun N.plus(arg: T) = map(this) { value -> add(arg, value) }
     operator fun N.minus(arg: T) = map(this) { value -> add(arg, -value) }
 
@@ -90,6 +91,7 @@ interface NDRing<T, R : Ring<T>, N : NDStructure<T>> : Ring<N>, NDSpace<T, R, N>
      */
     override fun multiply(a: N, b: N): N = combine(a, b) { aValue, bValue -> multiply(aValue, bValue) }
 
+    //TODO move to extensions after KEEP-176
     operator fun N.times(arg: T) = map(this) { value -> multiply(arg, value) }
     operator fun T.times(arg: N) = map(arg) { value -> multiply(this@times, value) }
 }
@@ -109,6 +111,7 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
      */
     override fun divide(a: N, b: N): N = combine(a, b) { aValue, bValue -> divide(aValue, bValue) }
 
+    //TODO move to extensions after KEEP-176
     operator fun N.div(arg: T) = map(this) { value -> divide(arg, value) }
     operator fun T.div(arg: N) = map(arg) { divide(it, this@div) }
 
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
index 8d09b6d7b..f7d603dab 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
@@ -1,5 +1,6 @@
 package scientifik.kmath.linear
 
+import scientifik.kmath.structures.Matrix
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
@@ -16,7 +17,7 @@ class MatrixTest {
     @Test
     fun testVectorToMatrix() {
         val vector = Vector.real(5) { it.toDouble() }
-        val matrix = vector.toMatrix()
+        val matrix = vector.asMatrix()
         assertEquals(4.0, matrix[4, 0])
     }
 
@@ -33,8 +34,8 @@ class MatrixTest {
         val vector1 = Vector.real(5) { it.toDouble() }
         val vector2 = Vector.real(5) { 5 - it.toDouble() }
 
-        val matrix1 = vector1.toMatrix()
-        val matrix2 = vector2.toMatrix().transpose()
+        val matrix1 = vector1.asMatrix()
+        val matrix2 = vector2.asMatrix().transpose()
         val product = MatrixContext.real.run { matrix1 dot matrix2 }
 
 
@@ -44,7 +45,7 @@ class MatrixTest {
 
     @Test
     fun testBuilder() {
-        val matrix = FeaturedMatrix.build<Double>(2, 3)(
+        val matrix = Matrix.build<Double>(2, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
index 8887d3a32..89f08ae5a 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
@@ -1,25 +1,28 @@
 package scientifik.kmath.linear
 
+import scientifik.kmath.structures.Matrix
+import kotlin.contracts.ExperimentalContracts
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
+@ExperimentalContracts
 class RealLUSolverTest {
+
     @Test
     fun testInvertOne() {
         val matrix = MatrixContext.real.one(2, 2)
-        val inverted = LUSolver.real.inverse(matrix)
+        val inverted = MatrixContext.real.inverse(matrix)
         assertEquals(matrix, inverted)
     }
 
     @Test
     fun testInvert() {
-        val matrix = FeaturedMatrix.square(
+        val matrix = Matrix.square(
             3.0, 1.0,
             1.0, 3.0
         )
 
-        val decomposed = LUSolver.real.decompose(matrix)
-        val decomposition = decomposed.getFeature<LUPDecomposition<Double>>()!!
+        val decomposition = MatrixContext.real.lup(matrix)
 
         //Check determinant
         assertEquals(8.0, decomposition.determinant)
@@ -29,9 +32,9 @@ class RealLUSolverTest {
             assertEquals(decomposition.p dot matrix, decomposition.l dot decomposition.u)
         }
 
-        val inverted = LUSolver.real.inverse(decomposed)
+        val inverted = MatrixContext.real.inverse(matrix)
 
-        val expected = FeaturedMatrix.square(
+        val expected = Matrix.square(
             0.375, -0.125,
             -0.125, 0.375
         )
diff --git a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
index b95f04887..74681ac48 100644
--- a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
+++ b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
@@ -2,10 +2,14 @@ package scientifik.kmath.linear
 
 import koma.extensions.fill
 import koma.matrix.MatrixFactory
+import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.Matrix
 
-class KomaMatrixContext<T : Any>(val factory: MatrixFactory<koma.matrix.Matrix<T>>) : MatrixContext<T>,
-    LinearSolver<T> {
+class KomaMatrixContext<T : Any>(
+    private val factory: MatrixFactory<koma.matrix.Matrix<T>>,
+    private val space: Space<T>
+) :
+    MatrixContext<T> {
 
     override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T) =
         KomaMatrix(factory.zeros(rows, columns).fill(initializer))
@@ -32,28 +36,38 @@ class KomaMatrixContext<T : Any>(val factory: MatrixFactory<koma.matrix.Matrix<T
     override fun Matrix<T>.unaryMinus() =
         KomaMatrix(this.toKoma().origin.unaryMinus())
 
-    override fun Matrix<T>.plus(b: Matrix<T>) =
-        KomaMatrix(this.toKoma().origin + b.toKoma().origin)
+    override fun add(a: Matrix<T>, b: Matrix<T>) =
+        KomaMatrix(a.toKoma().origin + b.toKoma().origin)
 
     override fun Matrix<T>.minus(b: Matrix<T>) =
         KomaMatrix(this.toKoma().origin - b.toKoma().origin)
 
+    override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
+        produce(a.rowNum, a.colNum) { i, j -> space.run { a[i, j] * k } }
+
     override fun Matrix<T>.times(value: T) =
         KomaMatrix(this.toKoma().origin * value)
 
+    companion object {
 
-    override fun solve(a: Matrix<T>, b: Matrix<T>) =
-        KomaMatrix(a.toKoma().origin.solve(b.toKoma().origin))
+    }
 
-    override fun inverse(a: Matrix<T>) =
-        KomaMatrix(a.toKoma().origin.inv())
 }
 
+fun <T : Any> KomaMatrixContext<T>.solve(a: Matrix<T>, b: Matrix<T>) =
+    KomaMatrix(a.toKoma().origin.solve(b.toKoma().origin))
+
+fun <T : Any> KomaMatrixContext<T>.solve(a: Matrix<T>, b: Point<T>) =
+    KomaVector(a.toKoma().origin.solve(b.toKoma().origin))
+
+fun <T : Any> KomaMatrixContext<T>.inverse(a: Matrix<T>) =
+    KomaMatrix(a.toKoma().origin.inv())
+
 class KomaMatrix<T : Any>(val origin: koma.matrix.Matrix<T>, features: Set<MatrixFeature>? = null) : FeaturedMatrix<T> {
     override val rowNum: Int get() = origin.numRows()
     override val colNum: Int get() = origin.numCols()
 
-    override val shape: IntArray get() = intArrayOf(origin.numRows(),origin.numCols())
+    override val shape: IntArray get() = intArrayOf(origin.numRows(), origin.numCols())
 
     override val features: Set<MatrixFeature> = features ?: setOf(
         object : DeterminantFeature<T> {