From 36609827ddb045f603a98d156308850339c142d4 Mon Sep 17 00:00:00 2001
From: Mikhail Zelenyy <mihail.zelenyy@phystech.edu>
Date: Wed, 23 Jan 2019 04:41:26 +0300
Subject: [PATCH 1/6] Create object with factory method for generating
 RealNDElement.

---
 .../kmath/structures/CreationRoutines.kt      | 169 ++++++++++++++++++
 1 file changed, 169 insertions(+)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt
new file mode 100644
index 000000000..69698a627
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt
@@ -0,0 +1,169 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.operations.RealField.power
+import kotlin.math.*
+
+
+object RealFactories{
+    /**
+     *  Create a NDArray filled with ones
+     */
+    fun ones(vararg shape: Int) = NDElement.real(shape){1.0}
+
+    /**
+     * Create a 2D NDArray, with ones on the diagonal and zeros elsewhere.
+     *
+     * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
+     */
+    fun eye(dim1: Int, dim2: Int, offset : Int = 0) = NDElement.real2D(dim1, dim2){i, j -> if (i == j + offset) 1.0 else 0.0}
+
+    /**
+     * An array with ones at and below the given diagonal and zeros elsewhere.
+     * T[i,j] == 1 for i <= j + offset
+     *
+     * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
+     */
+    fun triangle(dim1: Int, dim2: Int, offset : Int = 0) = NDElement.real2D(dim1, dim2){i, j -> if (i <= j + offset) 1.0 else 0.0}
+
+    /**
+     * Return evenly spaced values within a given interval.
+     *
+     * Values are generated within the half-open interval [start, stop) (in other words, the interval including start but excluding stop).
+     *  @param range use it like:
+     *           (start..stop) to step
+     */
+    fun range(range : Pair<ClosedFloatingPointRange<Double>,Double>) = NDElement.real1D(ceil((range.first.endInclusive - range.first.start)/range.second).toInt()){i-> range.first.start + i*range.second}
+
+    /**
+     *  Return evenly spaced numbers over a specified interval.
+     *  @param range use it like:
+     *           (start..stop) to number
+     *          start is starting value, finaly value depend from endPoint parameter
+     *  @param endPoint If True, right boundary of range is the last sample. Otherwise, it is not included.
+     */
+    fun linspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true): Pair<RealNDElement, Double> {
+        val div = if (endPoint) (range.second - 1) else range.second
+        val delta = range.first.start - range.first.endInclusive
+        if (range.second > 1){
+            val step = delta/div
+            if (step == 0.0){ error("Bad ranges: step = $step")}
+            val result = NDElement.real1D(range.second){
+                if ( endPoint and (it == range.second - 1) ){ range.first.endInclusive}
+                range.first.start + it*step
+            }
+            return  result to step
+        }
+        else{
+            val step = Double.NaN
+            return NDElement.real1D(1){range.first.start} to step
+        }
+
+    }
+
+    /**
+     * Return numbers spaced evenly on a log scale.
+     * @param range use it like:
+     *          (start..stop) to number
+     *          power(base,start) is starting value, endvalue depend from endPoint parameter
+     * @param endPoint If True, power(base,stop) is the last sample. Otherwise, it is not included.
+     * @param base - The base of the log space.
+     */
+    fun logspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true, base : Double = 10.0)  : RealNDElement {
+        val lin = linspace(range, endPoint).first
+        val fun_ = {x: Double -> power(base, x)}
+        return fun_(lin) // FIXME: RealNDElement.map return not suitable type ( `linspace(range, endPoint).first.map{power(base, it}`)
+    }
+    /**
+     *  Return numbers spaced evenly on a log scale (a geometric progression).
+     *
+     *  This is similar to [logspace], but with endpoints specified directly. Each output sample is a constant multiple of the previous.
+     *  @param range use it like:
+     *           (start..stop) to number
+     *          start is starting value, finaly value depend from endPoint parameter
+     *  @param endPoint If True, right boundary of range is the last sample. Otherwise, it is not included.
+     */
+    fun geomspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true) : RealNDElement{
+        var start = range.first.start
+        var stop = range.first.endInclusive
+        val num = range.second
+        if ( start == 0.0 || stop == 0.0){
+            error("Geometric sequence cannot include zero")
+        }
+        var outSign = 1.0
+        if (sign(start) == -1.0 && sign(stop) == -1.0){
+            start = -start
+            stop = -stop
+            outSign = -outSign
+        }
+
+        val log_ = logspace((log(start, 10.0)..log(stop, 10.0) to num), endPoint=endPoint)
+        val fun_ = {x:Double -> outSign*x}
+        return fun_(log_) // FIXME: `outSign*log_`  --- don't define times operator
+
+    }
+
+    /**
+     * Return specified diagonals of 2D NDArray.
+     *
+     * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
+     */
+    fun extractDiagonal(array : RealNDElement, offset: Int = 0): RealNDElement{
+        if (array.dimension != 2){
+            error("Input must be 2D NDArray")}
+        val size = min(array.shape[0], array.shape[0])
+        if (offset>=0){
+            return NDElement.real1D(size){i -> array[i, i+offset]}
+        }
+        else{
+            return NDElement.real1D(size){i -> array[i-offset, i]}
+        }
+
+    }
+
+    /**
+     *   Return a 2-D array with [array] on the [offset] diagonal.
+     *
+     *   @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
+     */
+    fun fromDiagonal(array : RealNDElement, offset: Int = 0): RealNDElement{
+        if (array.dimension != 1){
+            error("Input must be 1D NDArray")}
+        val size = array.shape[0]
+        if (offset>=0){
+            return NDElement.real2D(size, size+offset){
+                i, j -> if (i == j+offset) array[i] else 0.0
+            }
+        }
+        else{
+            return NDElement.real2D(size-offset, size){
+                i, j -> if (i-offset == j) array[j] else 0.0
+            }
+        }
+    }
+
+    /**
+     * Generate a [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix).
+     *
+     *  @param nCols --- number of columns, as default using length of [array]
+     *  @param increasing --- Order of the powers of the columns. If True, the powers increase from left to right, if False (the default) they are reversed. FIXME: Default order like numpy
+     */
+    fun vandermonde(array : RealNDElement, nCols: Int = 0, increasing: Boolean =false): RealNDElement{
+        if (array.dimension != 1){
+            error("Input must be 1D NDArray")}
+        var size = if (nCols ==0) array.shape[0] else nCols
+        if (increasing){
+            return NDElement.real2D(array.shape[0], size){
+                i, j -> power(array[i], j)
+            }
+        }else{
+            return NDElement.real2D(array.shape[0], size){
+                i, j -> power(array[i], size - j - 1)
+            }
+        }
+
+    }
+
+}
+
+
-- 
2.34.1


From 1b411872ef88ed3d96ed1f1aa1282276e1a072d3 Mon Sep 17 00:00:00 2001
From: Mikhail Zelenyy <mihail.zelenyy@phystech.edu>
Date: Wed, 23 Jan 2019 11:52:09 +0300
Subject: [PATCH 2/6] Reformat code and change some name

---
 .../kmath/structures/CreationRoutines.kt      | 128 ++++++++++--------
 1 file changed, 71 insertions(+), 57 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt
index 69698a627..51432bc5b 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/CreationRoutines.kt
@@ -1,22 +1,22 @@
 package scientifik.kmath.structures
 
-import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.RealField.power
 import kotlin.math.*
 
 
-object RealFactories{
+object RealFactory {
     /**
      *  Create a NDArray filled with ones
      */
-    fun ones(vararg shape: Int) = NDElement.real(shape){1.0}
+    fun ones(vararg shape: Int) = NDElement.real(shape) { 1.0 }
 
     /**
      * Create a 2D NDArray, with ones on the diagonal and zeros elsewhere.
      *
      * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
      */
-    fun eye(dim1: Int, dim2: Int, offset : Int = 0) = NDElement.real2D(dim1, dim2){i, j -> if (i == j + offset) 1.0 else 0.0}
+    fun eye(dim1: Int, dim2: Int, offset: Int = 0) =
+        NDElement.real2D(dim1, dim2) { i, j -> if (i == j + offset) 1.0 else 0.0 }
 
     /**
      * An array with ones at and below the given diagonal and zeros elsewhere.
@@ -24,7 +24,8 @@ object RealFactories{
      *
      * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
      */
-    fun triangle(dim1: Int, dim2: Int, offset : Int = 0) = NDElement.real2D(dim1, dim2){i, j -> if (i <= j + offset) 1.0 else 0.0}
+    fun triangle(dim1: Int, dim2: Int, offset: Int = 0) =
+        NDElement.real2D(dim1, dim2) { i, j -> if (i <= j + offset) 1.0 else 0.0 }
 
     /**
      * Return evenly spaced values within a given interval.
@@ -33,7 +34,8 @@ object RealFactories{
      *  @param range use it like:
      *           (start..stop) to step
      */
-    fun range(range : Pair<ClosedFloatingPointRange<Double>,Double>) = NDElement.real1D(ceil((range.first.endInclusive - range.first.start)/range.second).toInt()){i-> range.first.start + i*range.second}
+    fun range(range: Pair<ClosedFloatingPointRange<Double>, Double>) =
+        NDElement.real1D(ceil((range.first.endInclusive - range.first.start) / range.second).toInt()) { i -> range.first.start + i * range.second }
 
     /**
      *  Return evenly spaced numbers over a specified interval.
@@ -42,21 +44,27 @@ object RealFactories{
      *          start is starting value, finaly value depend from endPoint parameter
      *  @param endPoint If True, right boundary of range is the last sample. Otherwise, it is not included.
      */
-    fun linspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true): Pair<RealNDElement, Double> {
+    fun linSpace(
+        range: Pair<ClosedFloatingPointRange<Double>, Int>,
+        endPoint: Boolean = true
+    ): Pair<RealNDElement, Double> {
         val div = if (endPoint) (range.second - 1) else range.second
         val delta = range.first.start - range.first.endInclusive
-        if (range.second > 1){
-            val step = delta/div
-            if (step == 0.0){ error("Bad ranges: step = $step")}
-            val result = NDElement.real1D(range.second){
-                if ( endPoint and (it == range.second - 1) ){ range.first.endInclusive}
-                range.first.start + it*step
+        if (range.second > 1) {
+            val step = delta / div
+            if (step == 0.0) {
+                error("Bad ranges: step = $step")
             }
-            return  result to step
-        }
-        else{
+            val result = NDElement.real1D(range.second) {
+                if (endPoint and (it == range.second - 1)) {
+                    range.first.endInclusive
+                }
+                range.first.start + it * step
+            }
+            return result to step
+        } else {
             val step = Double.NaN
-            return NDElement.real1D(1){range.first.start} to step
+            return NDElement.real1D(1) { range.first.start } to step
         }
 
     }
@@ -69,37 +77,42 @@ object RealFactories{
      * @param endPoint If True, power(base,stop) is the last sample. Otherwise, it is not included.
      * @param base - The base of the log space.
      */
-    fun logspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true, base : Double = 10.0)  : RealNDElement {
-        val lin = linspace(range, endPoint).first
-        val fun_ = {x: Double -> power(base, x)}
-        return fun_(lin) // FIXME: RealNDElement.map return not suitable type ( `linspace(range, endPoint).first.map{power(base, it}`)
+    fun logSpace(
+        range: Pair<ClosedFloatingPointRange<Double>, Int>,
+        endPoint: Boolean = true,
+        base: Double = 10.0
+    ): RealNDElement {
+        val lin = linSpace(range, endPoint).first
+        val tempFun = { x: Double -> power(base, x) }
+        return tempFun(lin) // FIXME: RealNDElement.map return not suitable type ( `linSpace(range, endPoint).first.map{power(base, it}`)
     }
+
     /**
      *  Return numbers spaced evenly on a log scale (a geometric progression).
      *
-     *  This is similar to [logspace], but with endpoints specified directly. Each output sample is a constant multiple of the previous.
+     *  This is similar to [logSpace], but with endpoints specified directly. Each output sample is a constant multiple of the previous.
      *  @param range use it like:
      *           (start..stop) to number
      *          start is starting value, finaly value depend from endPoint parameter
      *  @param endPoint If True, right boundary of range is the last sample. Otherwise, it is not included.
      */
-    fun geomspace(range : Pair<ClosedFloatingPointRange<Double>,Int>, endPoint: Boolean = true) : RealNDElement{
+    fun geomSpace(range: Pair<ClosedFloatingPointRange<Double>, Int>, endPoint: Boolean = true): RealNDElement {
         var start = range.first.start
         var stop = range.first.endInclusive
         val num = range.second
-        if ( start == 0.0 || stop == 0.0){
+        if (start == 0.0 || stop == 0.0) {
             error("Geometric sequence cannot include zero")
         }
         var outSign = 1.0
-        if (sign(start) == -1.0 && sign(stop) == -1.0){
+        if (sign(start) == -1.0 && sign(stop) == -1.0) {
             start = -start
             stop = -stop
             outSign = -outSign
         }
 
-        val log_ = logspace((log(start, 10.0)..log(stop, 10.0) to num), endPoint=endPoint)
-        val fun_ = {x:Double -> outSign*x}
-        return fun_(log_) // FIXME: `outSign*log_`  --- don't define times operator
+        val logRange = logSpace((log(start, 10.0)..log(stop, 10.0) to num), endPoint = endPoint)
+        val function = { x: Double -> outSign * x }
+        return function(logRange) // FIXME: `outSign*log_`  --- don't define times operator
 
     }
 
@@ -108,15 +121,15 @@ object RealFactories{
      *
      * @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
      */
-    fun extractDiagonal(array : RealNDElement, offset: Int = 0): RealNDElement{
-        if (array.dimension != 2){
-            error("Input must be 2D NDArray")}
-        val size = min(array.shape[0], array.shape[0])
-        if (offset>=0){
-            return NDElement.real1D(size){i -> array[i, i+offset]}
+    fun extractDiagonal(array: RealNDElement, offset: Int = 0): RealNDElement {
+        if (array.dimension != 2) {
+            error("Input must be 2D NDArray")
         }
-        else{
-            return NDElement.real1D(size){i -> array[i-offset, i]}
+        val size = min(array.shape[0], array.shape[0])
+        if (offset >= 0) {
+            return NDElement.real1D(size) { i -> array[i, i + offset] }
+        } else {
+            return NDElement.real1D(size) { i -> array[i - offset, i] }
         }
 
     }
@@ -126,18 +139,18 @@ object RealFactories{
      *
      *   @param offset Index of the diagonal: 0 (the default) refers to the main diagonal, a positive value refers to an upper diagonal, and a negative value to a lower diagonal.
      */
-    fun fromDiagonal(array : RealNDElement, offset: Int = 0): RealNDElement{
-        if (array.dimension != 1){
-            error("Input must be 1D NDArray")}
-        val size = array.shape[0]
-        if (offset>=0){
-            return NDElement.real2D(size, size+offset){
-                i, j -> if (i == j+offset) array[i] else 0.0
-            }
+    fun fromDiagonal(array: RealNDElement, offset: Int = 0): RealNDElement {
+        if (array.dimension != 1) {
+            error("Input must be 1D NDArray")
         }
-        else{
-            return NDElement.real2D(size-offset, size){
-                i, j -> if (i-offset == j) array[j] else 0.0
+        val size = array.shape[0]
+        if (offset >= 0) {
+            return NDElement.real2D(size, size + offset) { i, j ->
+                if (i == j + offset) array[i] else 0.0
+            }
+        } else {
+            return NDElement.real2D(size - offset, size) { i, j ->
+                if (i - offset == j) array[j] else 0.0
             }
         }
     }
@@ -148,17 +161,18 @@ object RealFactories{
      *  @param nCols --- number of columns, as default using length of [array]
      *  @param increasing --- Order of the powers of the columns. If True, the powers increase from left to right, if False (the default) they are reversed. FIXME: Default order like numpy
      */
-    fun vandermonde(array : RealNDElement, nCols: Int = 0, increasing: Boolean =false): RealNDElement{
-        if (array.dimension != 1){
-            error("Input must be 1D NDArray")}
-        var size = if (nCols ==0) array.shape[0] else nCols
-        if (increasing){
-            return NDElement.real2D(array.shape[0], size){
-                i, j -> power(array[i], j)
+    fun vandermonde(array: RealNDElement, nCols: Int = 0, increasing: Boolean = false): RealNDElement {
+        if (array.dimension != 1) {
+            error("Input must be 1D NDArray")
+        }
+        val size = if (nCols == 0) array.shape[0] else nCols
+        if (increasing) {
+            return NDElement.real2D(array.shape[0], size) { i, j ->
+                power(array[i], j)
             }
-        }else{
-            return NDElement.real2D(array.shape[0], size){
-                i, j -> power(array[i], size - j - 1)
+        } else {
+            return NDElement.real2D(array.shape[0], size) { i, j ->
+                power(array[i], size - j - 1)
             }
         }
 
-- 
2.34.1


From 77bf8de4f1a5f8f686b3f2c0c71805eb641c67c3 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Wed, 23 Jan 2019 13:30:26 +0300
Subject: [PATCH 3/6] Optimized mapping functions for NDElements

---
 .../scientifik/kmath/structures/NDAlgebra.kt  |  9 +++++---
 .../scientifik/kmath/structures/NDElement.kt  | 23 ++++++++++++-------
 .../kmath/structures/RealNDField.kt           | 20 +++++++++++++---
 settings.gradle.kts                           |  3 +--
 4 files changed, 39 insertions(+), 16 deletions(-)

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 097d52723..7ea768c63 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
@@ -14,7 +14,7 @@ class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : Run
 /**
  * The base interface for all nd-algebra implementations
  * @param T the type of nd-structure element
- * @param C the type of the context
+ * @param C the type of the element context
  * @param N the type of the structure
  */
 interface NDAlgebra<T, C, N : NDStructure<T>> {
@@ -112,10 +112,13 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
     operator fun T.div(arg: N) = map(arg) { divide(it, this@div) }
 
     companion object {
+
+        private val realNDFieldCache = HashMap<IntArray, RealNDField>()
+
         /**
-         * Create a nd-field for [Double] values
+         * Create a nd-field for [Double] values or pull it from cache if it was created previously
          */
-        fun real(shape: IntArray) = RealNDField(shape)
+        fun real(shape: IntArray) = realNDFieldCache.getOrPut(shape){RealNDField(shape)}
 
         /**
          * Create a nd-field with boxing generic buffer
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
index 0fdb53f07..c97f959f3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
@@ -7,6 +7,9 @@ import scientifik.kmath.operations.Space
 
 /**
  * The root for all [NDStructure] based algebra elements. Does not implement algebra element root because of problems with recursive self-types
+ * @param T  the type of the element of the structure
+ * @param C the type of the context for the element
+ * @param N the type of the underlying [NDStructure]
  */
 interface NDElement<T, C, N : NDStructure<T>> : NDStructure<T> {
 
@@ -16,9 +19,6 @@ interface NDElement<T, C, N : NDStructure<T>> : NDStructure<T> {
 
     fun N.wrap(): NDElement<T, C, N>
 
-    fun mapIndexed(transform: C.(index: IntArray, T) -> T) = context.mapIndexed(unwrap(), transform).wrap()
-    fun map(transform: C.(T) -> T) = context.map(unwrap(), transform).wrap()
-
     companion object {
         /**
          * Create a optimized NDArray of doubles
@@ -61,10 +61,17 @@ interface NDElement<T, C, N : NDStructure<T>> : NDStructure<T> {
     }
 }
 
+
+fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.mapIndexed(transform: C.(index: IntArray, T) -> T) =
+    context.mapIndexed(unwrap(), transform).wrap()
+
+fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.map(transform: C.(T) -> T) = context.map(unwrap(), transform).wrap()
+
+
 /**
  * Element by element application of any operation on elements to the whole [NDElement]
  */
-operator fun <T, C> Function1<T, T>.invoke(ndElement: NDElement<T, C, *>) =
+operator fun <T, C, N : NDStructure<T>> Function1<T, T>.invoke(ndElement: NDElement<T, C, N>) =
     ndElement.map { value -> this@invoke(value) }
 
 /* plus and minus */
@@ -72,13 +79,13 @@ operator fun <T, C> Function1<T, T>.invoke(ndElement: NDElement<T, C, *>) =
 /**
  * Summation operation for [NDElement] and single element
  */
-operator fun <T, S : Space<T>> NDElement<T, S, *>.plus(arg: T) =
+operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.plus(arg: T) =
     map { value -> arg + value }
 
 /**
  * Subtraction operation between [NDElement] and single element
  */
-operator fun <T, S : Space<T>> NDElement<T, S, *>.minus(arg: T) =
+operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.minus(arg: T) =
     map { value -> arg - value }
 
 /* prod and div */
@@ -86,13 +93,13 @@ operator fun <T, S : Space<T>> NDElement<T, S, *>.minus(arg: T) =
 /**
  * Product operation for [NDElement] and single element
  */
-operator fun <T, R : Ring<T>> NDElement<T, R, *>.times(arg: T) =
+operator fun <T, R : Ring<T>, N : NDStructure<T>> NDElement<T, R, N>.times(arg: T) =
     map { value -> arg * value }
 
 /**
  * Division operation between [NDElement] and single element
  */
-operator fun <T, F : Field<T>> NDElement<T, F, *>.div(arg: T) =
+operator fun <T, F : Field<T>, N : NDStructure<T>> NDElement<T, F, N>.div(arg: T) =
     map { value -> arg / value }
 
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
index bc5832e1c..d652bb8a8 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
@@ -86,11 +86,25 @@ inline fun BufferedNDField<Double, RealField>.produceInline(crossinline initiali
     return BufferedNDFieldElement(this, DoubleBuffer(array))
 }
 
+/**
+ * Map one [RealNDElement] using function with indexes
+ */
+inline fun RealNDElement.mapIndexed(crossinline transform: RealField.(index: IntArray, Double) -> Double) =
+    context.produceInline { offset -> transform(strides.index(offset), buffer[offset]) }
+
+/**
+ * Map one [RealNDElement] using function without indexes
+ */
+inline fun RealNDElement.map(crossinline transform: RealField.(Double) -> Double): RealNDElement {
+    val array = DoubleArray(strides.linearSize) { offset -> RealField.transform(buffer[offset]) }
+    return BufferedNDFieldElement(context, DoubleBuffer(array))
+}
+
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
 operator fun Function1<Double, Double>.invoke(ndElement: RealNDElement) =
-    ndElement.context.produceInline { i -> invoke(ndElement.buffer[i]) }
+    ndElement.map { this@invoke(it) }
 
 
 /* plus and minus */
@@ -99,10 +113,10 @@ operator fun Function1<Double, Double>.invoke(ndElement: RealNDElement) =
  * Summation operation for [BufferedNDElement] and single element
  */
 operator fun RealNDElement.plus(arg: Double) =
-    context.produceInline { i -> buffer[i] + arg }
+    map { it + arg }
 
 /**
  * Subtraction operation between [BufferedNDElement] and single element
  */
 operator fun RealNDElement.minus(arg: Double) =
-    context.produceInline { i -> buffer[i] - arg }
+    map { it - arg }
diff --git a/settings.gradle.kts b/settings.gradle.kts
index a4464d01f..ca738647e 100644
--- a/settings.gradle.kts
+++ b/settings.gradle.kts
@@ -2,8 +2,7 @@ pluginManagement {
     repositories {
         mavenCentral()
         maven("https://plugins.gradle.org/m2/")
-        maven { setUrl("https://dl.bintray.com/kotlin/kotlin-eap") }
-        maven { setUrl("https://plugins.gradle.org/m2/") }
+        maven ("https://dl.bintray.com/kotlin/kotlin-eap")
     }
 }
 
-- 
2.34.1


From ce29e0e26c99cf25351fc1c1652a6567454ceabe Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 26 Jan 2019 19:38:18 +0300
Subject: [PATCH 4/6] Migrate to stable 1.3.20

---
 build.gradle.kts                              | 14 ++---
 kmath-core/build.gradle                       | 50 -----------------
 kmath-core/build.gradle.kts                   | 55 +++++++++++++++++++
 .../kotlin/scientifik/kmath/linear/Matrix.kt  |  2 +-
 .../kmath/linear/RealLUSolverTest.kt          |  4 +-
 kmath-koma/build.gradle.kts                   |  2 +-
 settings.gradle.kts                           |  2 +-
 7 files changed, 65 insertions(+), 64 deletions(-)
 delete mode 100644 kmath-core/build.gradle
 create mode 100644 kmath-core/build.gradle.kts

diff --git a/build.gradle.kts b/build.gradle.kts
index e5c4fc39b..6f0478fc2 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -1,14 +1,10 @@
 buildscript {
-    extra["kotlinVersion"] = "1.3.20-eap-100"
-    extra["ioVersion"] = "0.1.2"
-    extra["coroutinesVersion"] = "1.1.0"
-
-    val kotlinVersion: String by extra
-    val ioVersion: String by extra
-    val coroutinesVersion: String by extra
+    val kotlinVersion: String by rootProject.extra("1.3.20")
+    val ioVersion: String by rootProject.extra("0.1.2")
+    val coroutinesVersion: String by rootProject.extra("1.1.1")
 
     repositories {
-        maven("https://dl.bintray.com/kotlin/kotlin-eap")
+        //maven("https://dl.bintray.com/kotlin/kotlin-eap")
         jcenter()
     }
 
@@ -28,7 +24,7 @@ allprojects {
     apply(plugin = "com.jfrog.artifactory")
 
     group = "scientifik"
-    version = "0.0.3-dev-3"
+    version = "0.0.3-dev-4"
 
     repositories {
         maven("https://dl.bintray.com/kotlin/kotlin-eap")
diff --git a/kmath-core/build.gradle b/kmath-core/build.gradle
deleted file mode 100644
index b1a6cccb5..000000000
--- a/kmath-core/build.gradle
+++ /dev/null
@@ -1,50 +0,0 @@
-plugins {
-    id "org.jetbrains.kotlin.multiplatform"
-}
-
-kotlin {
-    jvm {
-        compilations["main"].kotlinOptions.jvmTarget = "1.8"
-        compilations["test"].kotlinOptions.jvmTarget = "1.8"
-    }
-    js()
-
-    sourceSets {
-        commonMain {
-            dependencies {
-                api 'org.jetbrains.kotlin:kotlin-stdlib-common'
-            }
-        }
-        commonTest {
-            dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-test-common'
-                implementation 'org.jetbrains.kotlin:kotlin-test-annotations-common'
-            }
-        }
-        jvmMain {
-            dependencies {
-                api 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
-            }
-        }
-        jvmTest {
-            dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-test'
-                implementation 'org.jetbrains.kotlin:kotlin-test-junit'
-            }
-        }
-        jsMain {
-            dependencies {
-                api 'org.jetbrains.kotlin:kotlin-stdlib-js'
-            }
-        }
-        jsTest {
-            dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-test-js'
-            }
-        }
-//        mingwMain {
-//        }
-//        mingwTest {
-//        }
-    }
-}
diff --git a/kmath-core/build.gradle.kts b/kmath-core/build.gradle.kts
new file mode 100644
index 000000000..2e753f624
--- /dev/null
+++ b/kmath-core/build.gradle.kts
@@ -0,0 +1,55 @@
+plugins {
+    kotlin("multiplatform")
+}
+
+
+kotlin {
+    jvm {
+        compilations.all {
+            kotlinOptions {
+                jvmTarget = "1.8"
+                freeCompilerArgs += "-progressive"
+            }
+        }
+    }
+    js()
+
+    sourceSets {
+        val commonMain by getting {
+            dependencies {
+                api(kotlin("stdlib"))
+            }
+        }
+        val commonTest by getting {
+            dependencies {
+                implementation(kotlin("test-common"))
+                implementation(kotlin("test-annotations-common"))
+            }
+        }
+        val jvmMain by getting {
+            dependencies {
+                api(kotlin("stdlib-jdk8"))
+            }
+        }
+        val jvmTest by getting {
+            dependencies {
+                implementation(kotlin("test"))
+                implementation(kotlin("test-junit"))
+            }
+        }
+        val jsMain by getting {
+            dependencies {
+                api(kotlin("stdlib-js"))
+            }
+        }
+        val jsTest by getting {
+            dependencies {
+                implementation(kotlin("test-js"))
+            }
+        }
+//        mingwMain {
+//        }
+//        mingwTest {
+//        }
+    }
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
index 5b4f3f12d..abe44e81a 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -167,7 +167,7 @@ interface Matrix<T : Any> : NDStructure<T> {
         /**
          * Build a square matrix from given elements.
          */
-        fun <T : Any> build(vararg elements: T): Matrix<T> {
+        fun <T : Any> square(vararg elements: T): Matrix<T> {
             val size: Int = sqrt(elements.size.toDouble()).toInt()
             if (size * size != elements.size) error("The number of elements ${elements.size} is not a full square")
             val buffer = elements.asBuffer()
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 ea104355c..bfa720369 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
@@ -13,7 +13,7 @@ class RealLUSolverTest {
 
     @Test
     fun testInvert() {
-        val matrix = Matrix.build(
+        val matrix = Matrix.square(
             3.0, 1.0,
             1.0, 3.0
         )
@@ -31,7 +31,7 @@ class RealLUSolverTest {
 
         val inverted = LUSolver.real.inverse(decomposed)
 
-        val expected = Matrix.build(
+        val expected = Matrix.square(
             0.375, -0.125,
             -0.125, 0.375
         )
diff --git a/kmath-koma/build.gradle.kts b/kmath-koma/build.gradle.kts
index 20ad02fbe..0e527e263 100644
--- a/kmath-koma/build.gradle.kts
+++ b/kmath-koma/build.gradle.kts
@@ -1,5 +1,5 @@
 plugins {
-    id("kotlin-multiplatform")
+    kotlin("multiplatform")
 }
 
 repositories {
diff --git a/settings.gradle.kts b/settings.gradle.kts
index ca738647e..cee985432 100644
--- a/settings.gradle.kts
+++ b/settings.gradle.kts
@@ -2,7 +2,7 @@ pluginManagement {
     repositories {
         mavenCentral()
         maven("https://plugins.gradle.org/m2/")
-        maven ("https://dl.bintray.com/kotlin/kotlin-eap")
+        //maven ("https://dl.bintray.com/kotlin/kotlin-eap")
     }
 }
 
-- 
2.34.1


From 569ff6357b42fdbf90616127c59a40a7564608a5 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 26 Jan 2019 22:04:28 +0300
Subject: [PATCH 5/6] Added matrix builder for small matrices

---
 .../kotlin/scientifik/kmath/linear/Matrix.kt  | 30 +++++++------------
 .../scientifik/kmath/linear/MatrixFeatures.kt | 30 +++++++++++++++++++
 .../scientifik/kmath/linear/MatrixTest.kt     | 10 +++++++
 kmath-koma/build.gradle.kts                   |  8 +++--
 4 files changed, 56 insertions(+), 22 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
index abe44e81a..68a81f077 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -107,26 +107,6 @@ interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
 }
 
-/**
- * A marker interface representing some matrix feature like diagonal, sparce, zero, etc. Features used to optimize matrix
- * operations performance in some cases.
- */
-interface MatrixFeature
-
-object DiagonalFeature : MatrixFeature
-
-object ZeroFeature : MatrixFeature
-
-object UnitFeature : MatrixFeature
-
-interface InverseMatrixFeature<T : Any> : MatrixFeature {
-    val inverse: Matrix<T>
-}
-
-interface DeterminantFeature<T : Any> : MatrixFeature {
-    val determinant: T
-}
-
 /**
  * Specialized 2-d structure
  */
@@ -173,6 +153,16 @@ interface Matrix<T : Any> : NDStructure<T> {
             val buffer = elements.asBuffer()
             return BufferMatrix(size, size, buffer)
         }
+
+        fun <T : Any> build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)
+    }
+}
+
+class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
+    operator fun invoke(vararg elements: T): Matrix<T> {
+        if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns")
+        val buffer = elements.asBuffer()
+        return BufferMatrix(rows, columns, buffer)
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt
new file mode 100644
index 000000000..d9f4e58ca
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt
@@ -0,0 +1,30 @@
+package scientifik.kmath.linear
+
+/**
+ * A marker interface representing some matrix feature like diagonal, sparce, zero, etc. Features used to optimize matrix
+ * operations performance in some cases.
+ */
+interface MatrixFeature
+
+/**
+ * The matrix with this feature is considered to have only diagonal non-null elements
+ */
+object DiagonalFeature : MatrixFeature
+
+/**
+ * Matix with this feature has all zero elements
+ */
+object ZeroFeature : MatrixFeature
+
+/**
+ * Matrix with this feature have unit elements on diagonal and zero elements in all other places
+ */
+object UnitFeature : MatrixFeature
+
+interface InverseMatrixFeature<T : Any> : MatrixFeature {
+    val inverse: Matrix<T>
+}
+
+interface DeterminantFeature<T : Any> : MatrixFeature {
+    val determinant: T
+}
\ No newline at end of file
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 2a00b84cd..61aa506c4 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
@@ -41,4 +41,14 @@ class MatrixTest {
         assertEquals(5.0, product[1, 0])
         assertEquals(6.0, product[2, 2])
     }
+
+    @Test
+    fun testBuilder() {
+        val matrix = Matrix.build<Double>(2, 3)(
+            1.0, 0.0, 0.0,
+            0.0, 1.0, 2.0
+        )
+
+        assertEquals(2.0, matrix[1, 2])
+    }
 }
\ No newline at end of file
diff --git a/kmath-koma/build.gradle.kts b/kmath-koma/build.gradle.kts
index 0e527e263..b95aaf3c8 100644
--- a/kmath-koma/build.gradle.kts
+++ b/kmath-koma/build.gradle.kts
@@ -8,8 +8,12 @@ repositories {
 
 kotlin {
     jvm {
-        compilations["main"].kotlinOptions.jvmTarget = "1.8"
-        compilations["test"].kotlinOptions.jvmTarget = "1.8"
+        compilations.all {
+            kotlinOptions {
+                jvmTarget = "1.8"
+                freeCompilerArgs += "-progressive"
+            }
+        }
     }
     js()
 
-- 
2.34.1


From a2ef50ab477125c4f2c477562fb0fda642cc54e8 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sun, 27 Jan 2019 09:45:32 +0300
Subject: [PATCH 6/6] Refactored Matrix features

---
 .../scientifik/kmath/linear/CMMatrix.kt       | 11 ++++--
 .../scientifik/kmath/linear/BufferMatrix.kt   |  3 ++
 .../kmath/linear/LUPDecomposition.kt          |  9 +++--
 .../kotlin/scientifik/kmath/linear/Matrix.kt  |  8 +++++
 .../scientifik/kmath/linear/MatrixFeatures.kt | 36 +++++++++++++++++--
 .../scientifik/kmath/linear/VirtualMatrix.kt  |  3 ++
 .../kmath/structures/ShortNDRing.kt           |  1 +
 .../scientifik.kmath.linear/KomaMatrix.kt     | 19 ++++++++--
 8 files changed, 78 insertions(+), 12 deletions(-)

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 aed4f4bb1..808b2768b 100644
--- a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
+++ b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
@@ -4,11 +4,16 @@ import org.apache.commons.math3.linear.*
 import org.apache.commons.math3.linear.RealMatrix
 import org.apache.commons.math3.linear.RealVector
 
-inline class CMMatrix(val origin: RealMatrix) : Matrix<Double> {
+class CMMatrix(val origin: RealMatrix, features: Set<MatrixFeature>? = null) : Matrix<Double> {
     override val rowNum: Int get() = origin.rowDimension
     override val colNum: Int get() = origin.columnDimension
 
-    override val features: Set<MatrixFeature> get() = emptySet()
+    override val features: Set<MatrixFeature> = features ?: sequence<MatrixFeature> {
+        if(origin is DiagonalMatrix) yield(DiagonalFeature)
+    }.toSet()
+
+    override fun suggestFeature(vararg features: MatrixFeature) =
+        CMMatrix(origin, this.features + features)
 
     override fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
 }
@@ -23,7 +28,7 @@ fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {
 
 fun RealMatrix.toMatrix() = CMMatrix(this)
 
-inline class CMVector(val origin: RealVector) : Point<Double> {
+class CMVector(val origin: RealVector) : Point<Double> {
     override val size: Int get() = origin.dimension
 
     override fun get(index: Int): Double = origin.getEntry(index)
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 5b65821db..3752f6db5 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
@@ -35,6 +35,9 @@ class BufferMatrix<T : Any>(
 
     override val shape: IntArray get() = intArrayOf(rowNum, colNum)
 
+    override fun suggestFeature(vararg features: MatrixFeature) =
+        BufferMatrix(rowNum, colNum, buffer, this.features + features)
+
     override fun get(index: IntArray): T = get(index[0], index[1])
 
     override fun get(i: Int, j: Int): T = buffer[i * colNum + j]
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 b2161baf2..85ddb8786 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
@@ -8,20 +8,19 @@ import scientifik.kmath.structures.MutableBufferFactory
 import scientifik.kmath.structures.NDStructure
 import scientifik.kmath.structures.get
 
-
 class LUPDecomposition<T : Comparable<T>>(
     private val elementContext: Ring<T>,
     internal val lu: NDStructure<T>,
     val pivot: IntArray,
     private val even: Boolean
-) : DeterminantFeature<T> {
+) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
 
     /**
      * Returns the matrix L of the decomposition.
      *
      * L is a lower-triangular matrix with [Ring.one] in diagonal
      */
-    val l: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
+    override val l: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1], setOf(LFeature)) { i, j ->
         when {
             j < i -> lu[i, j]
             j == i -> elementContext.one
@@ -35,7 +34,7 @@ class LUPDecomposition<T : Comparable<T>>(
      *
      * U is an upper-triangular matrix including the diagonal
      */
-    val u: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
+    override val u: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1], setOf(UFeature)) { i, j ->
         if (j >= i) lu[i, j] else elementContext.zero
     }
 
@@ -46,7 +45,7 @@ class LUPDecomposition<T : Comparable<T>>(
      * P is a sparse matrix with exactly one element set to [Ring.one] in
      * each row and each column, all other elements being set to [Ring.zero].
      */
-    val p: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
+    override val p: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
         if (j == pivot[i]) elementContext.one else elementContext.zero
     }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
index 68a81f077..533e77d61 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -116,6 +116,14 @@ interface Matrix<T : Any> : NDStructure<T> {
 
     val features: Set<MatrixFeature>
 
+    /**
+     * Suggest new feature for this matrix. The result is the new matrix that may or may not reuse existing data structure.
+     *
+     * The implementation does not guarantee to check that matrix actually have the feature, so one should be careful to
+     * add only those features that are valid.
+     */
+    fun suggestFeature(vararg features: MatrixFeature): Matrix<T>
+
     operator fun get(i: Int, j: Int): T
 
     override fun get(index: IntArray): T = get(index[0], index[1])
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt
index d9f4e58ca..6b45a14b1 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt
@@ -12,7 +12,7 @@ interface MatrixFeature
 object DiagonalFeature : MatrixFeature
 
 /**
- * Matix with this feature has all zero elements
+ * Matrix with this feature has all zero elements
  */
 object ZeroFeature : MatrixFeature
 
@@ -21,10 +21,42 @@ object ZeroFeature : MatrixFeature
  */
 object UnitFeature : MatrixFeature
 
+/**
+ * Inverted matrix feature
+ */
 interface InverseMatrixFeature<T : Any> : MatrixFeature {
     val inverse: Matrix<T>
 }
 
+/**
+ * A determinant container
+ */
 interface DeterminantFeature<T : Any> : MatrixFeature {
     val determinant: T
-}
\ No newline at end of file
+}
+
+@Suppress("FunctionName")
+fun <T: Any> DeterminantFeature(determinant: T) = object: DeterminantFeature<T>{
+    override val determinant: T = determinant
+}
+
+/**
+ * Lower triangular matrix
+ */
+object LFeature: MatrixFeature
+
+/**
+ * Upper triangular feature
+ */
+object UFeature: MatrixFeature
+
+/**
+ * TODO add documentation
+ */
+interface LUPDecompositionFeature<T : Any> : MatrixFeature {
+    val l: Matrix<T>
+    val u: Matrix<T>
+    val p: Matrix<T>
+}
+
+//TODO add sparse matrix feature
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt
index 98655ad48..1bab52902 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt
@@ -8,6 +8,9 @@ class VirtualMatrix<T : Any>(
 ) : Matrix<T> {
     override fun get(i: Int, j: Int): T = generator(i, j)
 
+    override fun suggestFeature(vararg features: MatrixFeature) =
+        VirtualMatrix(rowNum, colNum, this.features + features, generator)
+
     override fun equals(other: Any?): Boolean {
         if (this === other) return true
         if (other !is Matrix<*>) return false
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt
index 5c887f343..09e93483d 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt
@@ -15,6 +15,7 @@ class ShortNDRing(override val shape: IntArray) :
     override val zero by lazy { produce { ShortRing.zero } }
     override val one by lazy { produce { ShortRing.one } }
 
+    @Suppress("OVERRIDE_BY_INLINE")
     override inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Short): Buffer<Short> =
         ShortBuffer(ShortArray(size) { initializer(it) })
 
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 3b7894d18..343d5b8b9 100644
--- a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
+++ b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
@@ -48,10 +48,25 @@ class KomaMatrixContext<T : Any>(val factory: MatrixFactory<koma.matrix.Matrix<T
         KomaMatrix(a.toKoma().origin.inv())
 }
 
-inline class KomaMatrix<T : Any>(val origin: koma.matrix.Matrix<T>) : Matrix<T> {
+class KomaMatrix<T : Any>(val origin: koma.matrix.Matrix<T>, features: Set<MatrixFeature>? = null) :
+    Matrix<T> {
     override val rowNum: Int get() = origin.numRows()
     override val colNum: Int get() = origin.numCols()
-    override val features: Set<MatrixFeature> get() = emptySet()
+
+    override val features: Set<MatrixFeature> = features ?: setOf(
+        object : DeterminantFeature<T> {
+            override val determinant: T get() = origin.det()
+        },
+        object : LUPDecompositionFeature<T> {
+            private val lup by lazy { origin.LU() }
+            override val l: Matrix<T> get() = KomaMatrix(lup.second)
+            override val u: Matrix<T> get() = KomaMatrix(lup.third)
+            override val p: Matrix<T> get() = KomaMatrix(lup.first)
+        }
+    )
+
+    override fun suggestFeature(vararg features: MatrixFeature): Matrix<T> =
+        KomaMatrix(this.origin, this.features + features)
 
     override fun get(i: Int, j: Int): T = origin.getGeneric(i, j)
 }
-- 
2.34.1