From 2d2c5aa6840902db6c5688d19e2960b124b5e0d3 Mon Sep 17 00:00:00 2001 From: Andrei Kislitsyn Date: Mon, 15 Mar 2021 15:18:21 +0300 Subject: [PATCH] matrixhelper --- .../kscience/kmath/tensors/BufferedTensor.kt | 53 ++++++++ .../tensors/RealLinearOpsTensorAlgebra.kt | 116 ++++++++++-------- 2 files changed, 120 insertions(+), 49 deletions(-) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt index 29605024d..4d89996a0 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/BufferedTensor.kt @@ -13,6 +13,7 @@ public open class BufferedTensor( TensorStrides(shape), buffer ) { + /* public operator fun get(i: Int, j: Int): T { check(this.dimension == 2) { "Not matrix" } @@ -24,8 +25,60 @@ public open class BufferedTensor( this[intArrayOf(i, j)] = value } + */ } +//todo make generator mb nextMatrixIndex? +public class InnerMatrix(private val tensor: BufferedTensor){ + private var offset: Int = 0 + private val n : Int = tensor.shape.size + //stride? + private val step = tensor.shape[n - 1] * tensor.shape[n - 2] + + public operator fun get(i: Int, j: Int): T { + val index = tensor.strides.index(offset) + index[n - 2] = i + index[n - 1] = j + return tensor[index] + } + + public operator fun set(i: Int, j: Int, value: T): Unit { + val index = tensor.strides.index(offset) + index[n - 2] = i + index[n - 1] = j + tensor[index] = value + } + + public fun makeStep(){ + offset += step + } +} + +public class InnerVector(private val tensor: BufferedTensor){ + private var offset: Int = 0 + private val n : Int = tensor.shape.size + //stride? + private val step = tensor.shape[n - 1] + + public operator fun get(i: Int): T { + val index = tensor.strides.index(offset) + index[n - 1] = i + return tensor[index] + } + + public operator fun set(i: Int, value: T): Unit { + val index = tensor.strides.index(offset) + index[n - 1] = i + tensor[index] = value + } + + public fun makeStep(){ + offset += step + } +} + + +//todo default buffer = arrayOf(0)??? public class IntTensor( shape: IntArray, buffer: IntArray diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealLinearOpsTensorAlgebra.kt index dcd740356..ab39fc6ac 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealLinearOpsTensorAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/RealLinearOpsTensorAlgebra.kt @@ -2,8 +2,7 @@ package space.kscience.kmath.tensors public class RealLinearOpsTensorAlgebra : LinearOpsTensorAlgebra, - RealTensorAlgebra() -{ + RealTensorAlgebra() { override fun eye(n: Int): RealTensor { val shape = intArrayOf(n, n) val buffer = DoubleArray(n * n) { 0.0 } @@ -26,55 +25,71 @@ public class RealLinearOpsTensorAlgebra : override fun RealTensor.lu(): Pair { // todo checks - val lu = this.copy() - val m = this.shape[0] - val pivot = IntArray(m) - - - // Initialize permutation array and parity - for (row in 0 until m) pivot[row] = row - var even = true - - for (i in 0 until m) { - var maxA = -1.0 - var iMax = i - - for (k in i until m) { - val absA = kotlin.math.abs(lu[k, i]) - if (absA > maxA) { - maxA = absA - iMax = k - } - } - - //todo check singularity - - if (iMax != i) { - - val j = pivot[i] - pivot[i] = pivot[iMax] - pivot[iMax] = j - even != even - - for (k in 0 until m) { - val tmp = lu[i, k] - lu[i, k] = lu[iMax, k] - lu[iMax, k] = tmp - } - - } - - 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] - } - } + val luTensor = this.copy() + val lu = InnerMatrix(luTensor) + //stride TODO!!! move to generator? + var matCnt = 1 + for (i in 0 until this.shape.size - 2) { + matCnt *= this.shape[i] } - return Pair(lu, IntTensor(intArrayOf(m), pivot)) + val n = this.shape.size + val m = this.shape[n - 1] + val pivotsShape = IntArray(n - 1) { i -> + this.shape[i] + } + val pivotsTensor = IntTensor( + pivotsShape, + IntArray(matCnt * m) { 0 } + ) + val pivot = InnerVector(pivotsTensor) + for (i in 0 until matCnt) { + for (row in 0 until m) pivot[row] = row + + for (i in 0 until m) { + var maxA = -1.0 + var iMax = i + + for (k in i until m) { + val absA = kotlin.math.abs(lu[k, i]) + if (absA > maxA) { + maxA = absA + iMax = k + } + } + + //todo check singularity + + if (iMax != i) { + + val j = pivot[i] + pivot[i] = pivot[iMax] + pivot[iMax] = j + + for (k in 0 until m) { + val tmp = lu[i, k] + lu[i, k] = lu[iMax, k] + lu[iMax, k] = tmp + } + + } + + 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] + } + } + } + lu.makeStep() + pivot.makeStep() + } + + return Pair(luTensor, pivotsTensor) } override fun luPivot(lu: RealTensor, pivots: IntTensor): Triple { + TODO() + /* // todo checks val n = lu.shape[0] val p = lu.zeroesLike() @@ -97,9 +112,12 @@ public class RealLinearOpsTensorAlgebra : } } } - return Triple(p, l, u) + + return Triple(p, l, u)*/ } + + override fun RealTensor.det(): RealTensor { TODO("Not yet implemented") } @@ -127,5 +145,5 @@ public class RealLinearOpsTensorAlgebra : } -public inline fun RealLinearOpsTensorAlgebra(block: RealTensorAlgebra.() -> R): R = +public inline fun RealLinearOpsTensorAlgebra(block: RealLinearOpsTensorAlgebra.() -> R): R = RealLinearOpsTensorAlgebra().block() \ No newline at end of file