Separate linear algebra utils into dedicated module

kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt

@@ -1,7 +1,5 @@
 package space.kscience.kmath.tensors.core
 
-import space.kscience.kmath.nd.MutableStructure1D
-import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.structures.*
 import space.kscience.kmath.tensors.TensorStructure
 
@@ -18,9 +16,6 @@ public open class BufferedTensor<T>(
     public val numel: Int
         get() = linearStructure.size
 
-    internal constructor(tensor: BufferedTensor<T>) :
-            this(tensor.shape, tensor.buffer, tensor.bufferStart)
-
     override fun get(index: IntArray): T = buffer[bufferStart + linearStructure.offset(index)]
 
     override fun set(index: IntArray, value: T) {
@@ -35,39 +30,6 @@ public open class BufferedTensor<T>(
 
     override fun hashCode(): Int = 0
 
-    internal fun vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
-        val n = shape.size
-        val vectorOffset = shape[n - 1]
-        val vectorShape = intArrayOf(shape.last())
-        for (offset in 0 until numel step vectorOffset) {
-            val vector = BufferedTensor(vectorShape, buffer, offset)
-            yield(vector)
-        }
-    }
-
-    internal fun matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
-        check(shape.size >= 2) { "todo" }
-        val n = shape.size
-        val matrixOffset = shape[n - 1] * shape[n - 2]
-        val matrixShape = intArrayOf(shape[n - 2], shape[n - 1])
-        for (offset in 0 until numel step matrixOffset) {
-            val matrix = BufferedTensor(matrixShape, buffer, offset)
-            yield(matrix)
-        }
-    }
-
-    internal inline fun forEachVector(vectorAction: (BufferedTensor<T>) -> Unit): Unit {
-        for (vector in vectorSequence()) {
-            vectorAction(vector)
-        }
-    }
-
-    internal inline fun forEachMatrix(matrixAction: (BufferedTensor<T>) -> Unit): Unit {
-        for (matrix in matrixSequence()) {
-            matrixAction(matrix)
-        }
-    }
-
 }
 
 public class IntTensor internal constructor(

kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt

@@ -1,11 +1,8 @@
 package space.kscience.kmath.tensors.core
 
-import space.kscience.kmath.nd.MutableStructure1D
-import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
 import space.kscience.kmath.nd.as1D
 import space.kscience.kmath.nd.as2D
-import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
-import kotlin.math.sqrt
 
 public class DoubleLinearOpsTensorAlgebra :
     LinearOpsTensorAlgebra<Double, DoubleTensor, IntTensor>,
@@ -15,48 +12,6 @@ public class DoubleLinearOpsTensorAlgebra :
 
     override fun DoubleTensor.det(): DoubleTensor = detLU()
 
-    private inline fun luHelper(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>, m: Int) {
-        for (row in 0 until m) pivots[row] = row
-
-        for (i in 0 until m) {
-            var maxVal = -1.0
-            var maxInd = i
-
-            for (k in i until m) {
-                val absA = kotlin.math.abs(lu[k, i])
-                if (absA > maxVal) {
-                    maxVal = absA
-                    maxInd = k
-                }
-            }
-
-            //todo check singularity
-
-            if (maxInd != i) {
-
-                val j = pivots[i]
-                pivots[i] = pivots[maxInd]
-                pivots[maxInd] = j
-
-                for (k in 0 until m) {
-                    val tmp = lu[i, k]
-                    lu[i, k] = lu[maxInd, k]
-                    lu[maxInd, k] = tmp
-                }
-
-                pivots[m] += 1
-
-            }
-
-            for (j in i + 1 until m) {
-                lu[j, i] /= lu[i, i]
-                for (k in i + 1 until m) {
-                    lu[j, k] -= lu[j, i] * lu[i, k]
-                }
-            }
-        }
-    }
-
     override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
 
         checkSquareMatrix(shape)
@@ -80,37 +35,6 @@ public class DoubleLinearOpsTensorAlgebra :
 
     }
 
-    private inline fun pivInit(
-        p: MutableStructure2D<Double>,
-        pivot: MutableStructure1D<Int>,
-        n: Int
-    ) {
-        for (i in 0 until n) {
-            p[i, pivot[i]] = 1.0
-        }
-    }
-
-    private inline fun luPivotHelper(
-        l: MutableStructure2D<Double>,
-        u: MutableStructure2D<Double>,
-        lu: MutableStructure2D<Double>,
-        n: Int
-    ) {
-        for (i in 0 until n) {
-            for (j in 0 until n) {
-                if (i == j) {
-                    l[i, j] = 1.0
-                }
-                if (j < i) {
-                    l[i, j] = lu[i, j]
-                }
-                if (j >= i) {
-                    u[i, j] = lu[i, j]
-                }
-            }
-        }
-    }
-
     override fun luPivot(
         luTensor: DoubleTensor,
         pivotsTensor: IntTensor
@@ -139,27 +63,6 @@ public class DoubleLinearOpsTensorAlgebra :
 
     }
 
-    private inline fun choleskyHelper(
-        a: MutableStructure2D<Double>,
-        l: MutableStructure2D<Double>,
-        n: Int
-    ) {
-        for (i in 0 until n) {
-            for (j in 0 until i) {
-                var h = a[i, j]
-                for (k in 0 until j) {
-                    h -= l[i, k] * l[j, k]
-                }
-                l[i, j] = h / l[j, j]
-            }
-            var h = a[i, i]
-            for (j in 0 until i) {
-                h -= l[i, j] * l[i, j]
-            }
-            l[i, i] = sqrt(h)
-        }
-    }
-
     override fun DoubleTensor.cholesky(): DoubleTensor {
         // todo checks
         checkSquareMatrix(shape)
@@ -185,14 +88,6 @@ public class DoubleLinearOpsTensorAlgebra :
         TODO("ANDREI")
     }
 
-    private fun luMatrixDet(luTensor: MutableStructure2D<Double>, pivotsTensor: MutableStructure1D<Int>): Double {
-        val lu = luTensor.as2D()
-        val pivots = pivotsTensor.as1D()
-        val m = lu.shape[0]
-        val sign = if ((pivots[m] - m) % 2 == 0) 1.0 else -1.0
-        return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
-    }
-
     public fun DoubleTensor.detLU(): DoubleTensor {
         val (luTensor, pivotsTensor) = lu()
         val n = shape.size
@@ -213,33 +108,6 @@ public class DoubleLinearOpsTensorAlgebra :
         return detTensor
     }
 
-    private fun luMatrixInv(
-        lu: MutableStructure2D<Double>,
-        pivots: MutableStructure1D<Int>,
-        invMatrix: MutableStructure2D<Double>
-    ) {
-        val m = lu.shape[0]
-
-        for (j in 0 until m) {
-            for (i in 0 until m) {
-                if (pivots[i] == j) {
-                    invMatrix[i, j] = 1.0
-                }
-
-                for (k in 0 until i) {
-                    invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
-                }
-            }
-
-            for (i in m - 1 downTo 0) {
-                for (k in i + 1 until m) {
-                    invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
-                }
-                invMatrix[i, j] /= lu[i, i]
-            }
-        }
-    }
-
     public fun DoubleTensor.invLU(): DoubleTensor {
         val (luTensor, pivotsTensor) = lu()
         val invTensor = luTensor.zeroesLike()

kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt

@@ -1,8 +1,7 @@
 package space.kscience.kmath.tensors.core
 
-import space.kscience.kmath.nd.MutableStructure2D
-import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
+import space.kscience.kmath.nd.as2D
 import kotlin.math.abs
 
 public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, DoubleTensor> {
@@ -230,23 +229,6 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
         return this.view(other.shape)
     }
 
-    private inline fun dotHelper(
-        a: MutableStructure2D<Double>,
-        b: MutableStructure2D<Double>,
-        res: MutableStructure2D<Double>,
-        l: Int, m: Int, n: Int
-    ) {
-        for (i in 0 until l) {
-            for (j in 0 until n) {
-                var curr = 0.0
-                for (k in 0 until m) {
-                    curr += a[i, k] * b[k, j]
-                }
-                res[i, j] = curr
-            }
-        }
-    }
-
     override fun DoubleTensor.dot(other: DoubleTensor): DoubleTensor {
         if (this.shape.size == 1 && other.shape.size == 1) {
             return DoubleTensor(intArrayOf(1), doubleArrayOf(this.times(other).buffer.array().sum()))

kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt

@@ -0,0 +1,188 @@
+package space.kscience.kmath.tensors.core
+
+import space.kscience.kmath.nd.MutableStructure1D
+import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.as1D
+import space.kscience.kmath.nd.as2D
+import kotlin.math.sqrt
+
+
+internal inline fun <T> BufferedTensor<T>.vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
+    val n = shape.size
+    val vectorOffset = shape[n - 1]
+    val vectorShape = intArrayOf(shape.last())
+    for (offset in 0 until numel step vectorOffset) {
+        val vector = BufferedTensor(vectorShape, buffer, offset)
+        yield(vector)
+    }
+}
+
+internal inline fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
+    check(shape.size >= 2) { "todo" }
+    val n = shape.size
+    val matrixOffset = shape[n - 1] * shape[n - 2]
+    val matrixShape = intArrayOf(shape[n - 2], shape[n - 1])
+    for (offset in 0 until numel step matrixOffset) {
+        val matrix = BufferedTensor(matrixShape, buffer, offset)
+        yield(matrix)
+    }
+}
+
+internal inline fun <T> BufferedTensor<T>.forEachVector(vectorAction: (BufferedTensor<T>) -> Unit): Unit {
+    for (vector in vectorSequence()) {
+        vectorAction(vector)
+    }
+}
+
+internal inline fun <T> BufferedTensor<T>.forEachMatrix(matrixAction: (BufferedTensor<T>) -> Unit): Unit {
+    for (matrix in matrixSequence()) {
+        matrixAction(matrix)
+    }
+}
+
+
+internal inline fun dotHelper(
+    a: MutableStructure2D<Double>,
+    b: MutableStructure2D<Double>,
+    res: MutableStructure2D<Double>,
+    l: Int, m: Int, n: Int
+) {
+    for (i in 0 until l) {
+        for (j in 0 until n) {
+            var curr = 0.0
+            for (k in 0 until m) {
+                curr += a[i, k] * b[k, j]
+            }
+            res[i, j] = curr
+        }
+    }
+}
+
+internal inline fun luHelper(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>, m: Int) {
+    for (row in 0 until m) pivots[row] = row
+
+    for (i in 0 until m) {
+        var maxVal = -1.0
+        var maxInd = i
+
+        for (k in i until m) {
+            val absA = kotlin.math.abs(lu[k, i])
+            if (absA > maxVal) {
+                maxVal = absA
+                maxInd = k
+            }
+        }
+
+        //todo check singularity
+
+        if (maxInd != i) {
+
+            val j = pivots[i]
+            pivots[i] = pivots[maxInd]
+            pivots[maxInd] = j
+
+            for (k in 0 until m) {
+                val tmp = lu[i, k]
+                lu[i, k] = lu[maxInd, k]
+                lu[maxInd, k] = tmp
+            }
+
+            pivots[m] += 1
+
+        }
+
+        for (j in i + 1 until m) {
+            lu[j, i] /= lu[i, i]
+            for (k in i + 1 until m) {
+                lu[j, k] -= lu[j, i] * lu[i, k]
+            }
+        }
+    }
+}
+
+internal inline fun pivInit(
+    p: MutableStructure2D<Double>,
+    pivot: MutableStructure1D<Int>,
+    n: Int
+) {
+    for (i in 0 until n) {
+        p[i, pivot[i]] = 1.0
+    }
+}
+
+internal inline fun luPivotHelper(
+    l: MutableStructure2D<Double>,
+    u: MutableStructure2D<Double>,
+    lu: MutableStructure2D<Double>,
+    n: Int
+) {
+    for (i in 0 until n) {
+        for (j in 0 until n) {
+            if (i == j) {
+                l[i, j] = 1.0
+            }
+            if (j < i) {
+                l[i, j] = lu[i, j]
+            }
+            if (j >= i) {
+                u[i, j] = lu[i, j]
+            }
+        }
+    }
+}
+
+internal inline fun choleskyHelper(
+    a: MutableStructure2D<Double>,
+    l: MutableStructure2D<Double>,
+    n: Int
+) {
+    for (i in 0 until n) {
+        for (j in 0 until i) {
+            var h = a[i, j]
+            for (k in 0 until j) {
+                h -= l[i, k] * l[j, k]
+            }
+            l[i, j] = h / l[j, j]
+        }
+        var h = a[i, i]
+        for (j in 0 until i) {
+            h -= l[i, j] * l[i, j]
+        }
+        l[i, i] = sqrt(h)
+    }
+}
+
+internal inline fun luMatrixDet(luTensor: MutableStructure2D<Double>, pivotsTensor: MutableStructure1D<Int>): Double {
+    val lu = luTensor.as2D()
+    val pivots = pivotsTensor.as1D()
+    val m = lu.shape[0]
+    val sign = if ((pivots[m] - m) % 2 == 0) 1.0 else -1.0
+    return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
+}
+
+internal inline fun luMatrixInv(
+    lu: MutableStructure2D<Double>,
+    pivots: MutableStructure1D<Int>,
+    invMatrix: MutableStructure2D<Double>
+) {
+    val m = lu.shape[0]
+
+    for (j in 0 until m) {
+        for (i in 0 until m) {
+            if (pivots[i] == j) {
+                invMatrix[i, j] = 1.0
+            }
+
+            for (k in 0 until i) {
+                invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
+            }
+        }
+
+        for (i in m - 1 downTo 0) {
+            for (k in i + 1 until m) {
+                invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
+            }
+            invMatrix[i, j] /= lu[i, i]
+        }
+    }
+}