From 10739e0d04b5b2a8a2d999723d02bf8691b3ede1 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 18 Feb 2024 12:27:46 +0300 Subject: [PATCH] Performance fixes --- .../kscience/kmath/benchmarks/DotBenchmark.kt | 13 +- .../benchmarks/MatrixInverseBenchmark.kt | 20 +- buildSrc/build.gradle.kts | 2 +- .../kmath/ejml/codegen/ejmlCodegen.kt | 443 ------------------ .../kscience/kmath/linear/dotPerformance.kt | 15 +- .../kscience/kmath/linear/lupPerformance.kt | 27 ++ .../kscience/kmath/series/seriesBuilder.kt | 6 +- gradle/wrapper/gradle-wrapper.properties | 2 +- .../kmath/linear/Float64LinearSpace.kt | 6 +- .../kscience/kmath/linear/LupDecomposition.kt | 12 +- .../kmath/operations/BufferAlgebra.kt | 10 +- .../space/kscience/kmath/structures/Buffer.kt | 2 + .../kmath/structures/BufferAccessor2D.kt | 13 +- .../kmath/structures/Float64Buffer.kt | 1 + .../kmath/structures/MutableBuffer.kt | 1 - .../linear/Float64ParallelLinearSpace.kt | 125 +++++ .../kmath/structures/parallelMutableBuffer.kt | 58 +++ .../kscience/kmath/stat/ValueAndErrorField.kt | 2 +- 18 files changed, 262 insertions(+), 496 deletions(-) delete mode 100644 buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt create mode 100644 examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt create mode 100644 kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt create mode 100644 kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt index 51be23192..a157a67b2 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt @@ -11,12 +11,11 @@ import kotlinx.benchmark.Scope import kotlinx.benchmark.State import space.kscience.kmath.commons.linear.CMLinearSpace import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM +import space.kscience.kmath.linear.Float64ParallelLinearSpace import space.kscience.kmath.linear.invoke import space.kscience.kmath.linear.linearSpace import space.kscience.kmath.operations.Float64Field -import space.kscience.kmath.operations.invoke import space.kscience.kmath.tensorflow.produceWithTF -import space.kscience.kmath.tensors.core.DoubleTensorAlgebra import space.kscience.kmath.tensors.core.tensorAlgebra import kotlin.random.Random @@ -72,12 +71,12 @@ internal class DotBenchmark { } @Benchmark - fun tensorDot(blackhole: Blackhole) = with(Float64Field.tensorAlgebra) { + fun multikDot(blackhole: Blackhole) = with(multikAlgebra) { blackhole.consume(matrix1 dot matrix2) } @Benchmark - fun multikDot(blackhole: Blackhole) = with(multikAlgebra) { + fun tensorDot(blackhole: Blackhole) = with(Float64Field.tensorAlgebra) { blackhole.consume(matrix1 dot matrix2) } @@ -87,12 +86,8 @@ internal class DotBenchmark { } @Benchmark - fun doubleDot(blackhole: Blackhole) = with(Float64Field.linearSpace) { + fun parallelDot(blackhole: Blackhole) = with(Float64ParallelLinearSpace) { blackhole.consume(matrix1 dot matrix2) } - @Benchmark - fun doubleTensorDot(blackhole: Blackhole) = DoubleTensorAlgebra.invoke { - blackhole.consume(matrix1 dot matrix2) - } } diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt index 01513ac55..1a9e09013 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt @@ -15,6 +15,7 @@ import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM import space.kscience.kmath.linear.invoke import space.kscience.kmath.linear.linearSpace import space.kscience.kmath.linear.lupSolver +import space.kscience.kmath.linear.parallel import space.kscience.kmath.operations.algebra import kotlin.random.Random @@ -38,16 +39,19 @@ internal class MatrixInverseBenchmark { } @Benchmark - fun cmLUPInversion(blackhole: Blackhole) { - CMLinearSpace { - blackhole.consume(lupSolver().inverse(matrix)) - } + fun kmathParallelLupInversion(blackhole: Blackhole) { + blackhole.consume(Double.algebra.linearSpace.parallel.lupSolver().inverse(matrix)) } @Benchmark - fun ejmlInverse(blackhole: Blackhole) { - EjmlLinearSpaceDDRM { - blackhole.consume(matrix.toEjml().inverted()) - } + fun cmLUPInversion(blackhole: Blackhole) = CMLinearSpace { + blackhole.consume(lupSolver().inverse(matrix)) } + + + @Benchmark + fun ejmlInverse(blackhole: Blackhole) = EjmlLinearSpaceDDRM { + blackhole.consume(matrix.toEjml().inverted()) + } + } diff --git a/buildSrc/build.gradle.kts b/buildSrc/build.gradle.kts index 08db49244..89bb32b15 100644 --- a/buildSrc/build.gradle.kts +++ b/buildSrc/build.gradle.kts @@ -17,7 +17,7 @@ val benchmarksVersion = spclibs.versions.kotlinx.benchmark.get() dependencies { api("space.kscience:gradle-tools:$toolsVersion") //plugins form benchmarks - api("org.jetbrains.kotlinx:kotlinx-benchmark-plugin:0.4.9") + api("org.jetbrains.kotlinx:kotlinx-benchmark-plugin:$benchmarksVersion") //api("org.jetbrains.kotlin:kotlin-allopen:$kotlinVersion") //to be used inside build-script only //implementation(spclibs.kotlinx.serialization.json) diff --git a/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt b/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt deleted file mode 100644 index fbe891508..000000000 --- a/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt +++ /dev/null @@ -1,443 +0,0 @@ -/* - * Copyright 2018-2024 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. - */ - -@file:Suppress("KDocUnresolvedReference") - -package space.kscience.kmath.ejml.codegen - -import org.intellij.lang.annotations.Language -import java.io.File - -private fun Appendable.appendEjmlVector(type: String, ejmlMatrixType: String) { - @Language("kotlin") val text = """/** - * [EjmlVector] specialization for [$type]. - */ -public class Ejml${type}Vector(override val origin: M) : EjmlVector<$type, M>(origin) { - init { - require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" } - } - - override val type: SafeType<${type}> get() = safeTypeOf() - - override operator fun get(index: Int): $type = origin[0, index] -}""" - appendLine(text) - appendLine() -} - -private fun Appendable.appendEjmlMatrix(type: String, ejmlMatrixType: String) { - val text = """/** - * [EjmlMatrix] specialization for [$type]. - */ -public class Ejml${type}Matrix(override val origin: M) : EjmlMatrix<$type, M>(origin) { - override val type: SafeType<${type}> get() = safeTypeOf() - - override operator fun get(i: Int, j: Int): $type = origin[i, j] -}""" - appendLine(text) - appendLine() -} - -private fun Appendable.appendEjmlLinearSpace( - type: String, - kmathAlgebra: String, - ejmlMatrixParentTypeMatrix: String, - ejmlMatrixType: String, - ejmlMatrixDenseType: String, - ops: String, - denseOps: String, - isDense: Boolean, -) { - @Language("kotlin") val text = """ - -/** - * [EjmlLinearSpace] implementation based on [CommonOps_$ops], [DecompositionFactory_${ops}] operations and - * [${ejmlMatrixType}] matrices. - */ -public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra}, $ejmlMatrixType>() { - /** - * The [${kmathAlgebra}] reference. - */ - override val elementAlgebra: $kmathAlgebra get() = $kmathAlgebra - - override val type: SafeType<${type}> get() = safeTypeOf() - - @Suppress("UNCHECKED_CAST") - override fun Matrix<${type}>.toEjml(): Ejml${type}Matrix<${ejmlMatrixType}> = when { - this is Ejml${type}Matrix<*> && origin is $ejmlMatrixType -> this as Ejml${type}Matrix<${ejmlMatrixType}> - else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } - } - - @Suppress("UNCHECKED_CAST") - override fun Point<${type}>.toEjml(): Ejml${type}Vector<${ejmlMatrixType}> = when { - this is Ejml${type}Vector<*> && origin is $ejmlMatrixType -> this as Ejml${type}Vector<${ejmlMatrixType}> - else -> Ejml${type}Vector(${ejmlMatrixType}(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } - }) - } - - override fun buildMatrix( - rows: Int, - columns: Int, - initializer: ${kmathAlgebra}.(i: Int, j: Int) -> ${type}, - ): Ejml${type}Matrix<${ejmlMatrixType}> = ${ejmlMatrixType}(rows, columns).also { - (0 until rows).forEach { row -> - (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } - } - }.wrapMatrix() - - override fun buildVector( - size: Int, - initializer: ${kmathAlgebra}.(Int) -> ${type}, - ): Ejml${type}Vector<${ejmlMatrixType}> = Ejml${type}Vector(${ejmlMatrixType}(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } - }) - - private fun T.wrapMatrix() = Ejml${type}Matrix(this) - private fun T.wrapVector() = Ejml${type}Vector(this) - - override fun Matrix<${type}>.unaryMinus(): Matrix<${type}> = this * elementAlgebra { -one } - - override fun Matrix<${type}>.dot(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.mult(toEjml().origin, other.toEjml().origin, out) - return out.wrapMatrix() - } - - override fun Matrix<${type}>.dot(vector: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.mult(toEjml().origin, vector.toEjml().origin, out) - return out.wrapVector() - } - - override operator fun Matrix<${type}>.minus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - - CommonOps_${ops}.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out,${ - if (isDense) "" else - """ - null, - null,""" - } - ) - - return out.wrapMatrix() - } - - override operator fun Matrix<${type}>.times(value: ${type}): Ejml${type}Matrix<${ejmlMatrixType}> { - val res = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.scale(value, toEjml().origin, res) - return res.wrapMatrix() - } - - override fun Point<${type}>.unaryMinus(): Ejml${type}Vector<${ejmlMatrixType}> { - val res = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.changeSign(toEjml().origin, res) - return res.wrapVector() - } - - override fun Matrix<${type}>.plus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - - CommonOps_${ops}.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out,${ - if (isDense) "" else - """ - null, - null,""" - } - ) - - return out.wrapMatrix() - } - - override fun Point<${type}>.plus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - - CommonOps_${ops}.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out,${ - if (isDense) "" else - """ - null, - null,""" - } - ) - - return out.wrapVector() - } - - override fun Point<${type}>.minus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> { - val out = ${ejmlMatrixType}(1, 1) - - CommonOps_${ops}.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out,${ - if (isDense) "" else - """ - null, - null,""" - } - ) - - return out.wrapVector() - } - - override fun ${type}.times(m: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> = m * this - - override fun Point<${type}>.times(value: ${type}): Ejml${type}Vector<${ejmlMatrixType}> { - val res = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.scale(value, toEjml().origin, res) - return res.wrapVector() - } - - override fun ${type}.times(v: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> = v * this - - @UnstableKMathAPI - override fun computeFeature(structure: Matrix<${type}>, type: KClass): F? { - structure.getFeature(type)?.let { return it } - val origin = structure.toEjml().origin - - return when (type) { - ${ - if (isDense) - """ InverseMatrixFeature::class -> object : InverseMatrixFeature<${type}> { - override val inverse: Matrix<${type}> by lazy { - val res = origin.copy() - CommonOps_${ops}.invert(res) - res.wrapMatrix() - } - } - - DeterminantFeature::class -> object : DeterminantFeature<${type}> { - override val determinant: $type by lazy { CommonOps_${ops}.det(origin) } - } - - SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<${type}> { - private val svd by lazy { - DecompositionFactory_${ops}.svd(origin.numRows, origin.numCols, true, true, false) - .apply { decompose(origin.copy()) } - } - - override val u: Matrix<${type}> by lazy { svd.getU(null, false).wrapMatrix() } - override val s: Matrix<${type}> by lazy { svd.getW(null).wrapMatrix() } - override val v: Matrix<${type}> by lazy { svd.getV(null, false).wrapMatrix() } - override val singularValues: Point<${type}> by lazy { ${type}Buffer(svd.singularValues) } - } - - QRDecompositionFeature::class -> object : QRDecompositionFeature<${type}> { - private val qr by lazy { - DecompositionFactory_${ops}.qr().apply { decompose(origin.copy()) } - } - - override val q: Matrix<${type}> by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> { - override val l: Matrix<${type}> by lazy { - val cholesky = - DecompositionFactory_${ops}.chol(structure.rowNum, true).apply { decompose(origin.copy()) } - - cholesky.getT(null).wrapMatrix().withFeature(LFeature) - } - } - - LupDecompositionFeature::class -> object : LupDecompositionFeature<${type}> { - private val lup by lazy { - DecompositionFactory_${ops}.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) } - } - - override val l: Matrix<${type}> by lazy { - lup.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix<${type}> by lazy { - lup.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val p: Matrix<${type}> by lazy { lup.getRowPivot(null).wrapMatrix() } - }""" else """ QRDecompositionFeature::class -> object : QRDecompositionFeature<$type> { - private val qr by lazy { - DecompositionFactory_${ops}.qr(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val q: Matrix<${type}> by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> { - override val l: Matrix<${type}> by lazy { - val cholesky = - DecompositionFactory_${ops}.cholesky().apply { decompose(origin.copy()) } - - (cholesky.getT(null) as ${ejmlMatrixParentTypeMatrix}).wrapMatrix().withFeature(LFeature) - } - } - - LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object : - LUDecompositionFeature<${type}>, DeterminantFeature<${type}>, InverseMatrixFeature<${type}> { - private val lu by lazy { - DecompositionFactory_${ops}.lu(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val l: Matrix<${type}> by lazy { - lu.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix<${type}> by lazy { - lu.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val inverse: Matrix<${type}> by lazy { - var a = origin - val inverse = ${ejmlMatrixDenseType}(1, 1) - val solver = LinearSolverFactory_${ops}.lu(FillReducing.NONE) - if (solver.modifiesA()) a = a.copy() - val i = CommonOps_${denseOps}.identity(a.numRows) - solver.solve(i, inverse) - inverse.wrapMatrix() - } - - override val determinant: $type by lazy { elementAlgebra.number(lu.computeDeterminant().real) } - }""" - } - - else -> null - }?.let{ - type.cast(it) - } - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p matrix. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix<${type}>, b: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> { - val res = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res) - return res.wrapMatrix() - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p vector. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix<${type}>, b: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> { - val res = ${ejmlMatrixType}(1, 1) - CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res) - return Ejml${type}Vector(res) - } -}""" - appendLine(text) - appendLine() -} - - -/** - * Generates routine EJML classes. - */ -fun ejmlCodegen(outputFile: String): Unit = File(outputFile).run { - parentFile.mkdirs() - - writer().use { - it.appendLine("/*") - it.appendLine(" * Copyright 2018-2024 KMath contributors.") - it.appendLine(" * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.") - it.appendLine(" */") - it.appendLine() - it.appendLine("/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */") - it.appendLine() - it.appendLine("package space.kscience.kmath.ejml") - it.appendLine() - it.appendLine("""import org.ejml.data.* -import org.ejml.dense.row.CommonOps_DDRM -import org.ejml.dense.row.CommonOps_FDRM -import org.ejml.dense.row.factory.DecompositionFactory_DDRM -import org.ejml.dense.row.factory.DecompositionFactory_FDRM -import org.ejml.sparse.FillReducing -import org.ejml.sparse.csc.CommonOps_DSCC -import org.ejml.sparse.csc.CommonOps_FSCC -import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC -import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC -import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC -import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC -import space.kscience.attributes.SafeType -import space.kscience.attributes.safeTypeOf -import space.kscience.kmath.linear.* -import space.kscience.kmath.linear.Matrix -import space.kscience.kmath.UnstableKMathAPI -import space.kscience.kmath.nd.StructureFeature -import space.kscience.kmath.structures.Float64 -import space.kscience.kmath.structures.Float32 -import space.kscience.kmath.operations.Float64Field -import space.kscience.kmath.operations.Float32Field -import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.operations.FloatField -import space.kscience.kmath.operations.invoke -import space.kscience.kmath.structures.Float64Buffer -import space.kscience.kmath.structures.Float32Buffer -import space.kscience.kmath.structures.DoubleBuffer -import space.kscience.kmath.structures.FloatBuffer -import kotlin.reflect.KClass -import kotlin.reflect.cast""") - it.appendLine() - it.appendEjmlVector("Double", "DMatrix") - it.appendEjmlVector("Float", "FMatrix") - it.appendEjmlMatrix("Double", "DMatrix") - it.appendEjmlMatrix("Float", "FMatrix") - it.appendEjmlLinearSpace("Double", "Float64Field", "DMatrix", "DMatrixRMaj", "DMatrixRMaj", "DDRM", "DDRM", true) - it.appendEjmlLinearSpace("Float", "Float32Field", "FMatrix", "FMatrixRMaj", "FMatrixRMaj", "FDRM", "FDRM", true) - - it.appendEjmlLinearSpace( - type = "Double", - kmathAlgebra = "Float64Field", - ejmlMatrixParentTypeMatrix = "DMatrix", - ejmlMatrixType = "DMatrixSparseCSC", - ejmlMatrixDenseType = "DMatrixRMaj", - ops = "DSCC", - denseOps = "DDRM", - isDense = false, - ) - - it.appendEjmlLinearSpace( - type = "Float", - kmathAlgebra = "Float32Field", - ejmlMatrixParentTypeMatrix = "FMatrix", - ejmlMatrixType = "FMatrixSparseCSC", - ejmlMatrixDenseType = "FMatrixRMaj", - ops = "FSCC", - denseOps = "FDRM", - isDense = false, - ) - } -} diff --git a/examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt b/examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt index 8874ca685..c59ffeae2 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt @@ -5,29 +5,24 @@ package space.kscience.kmath.linear -import space.kscience.kmath.operations.algebra import kotlin.random.Random -import kotlin.time.ExperimentalTime import kotlin.time.measureTime -@OptIn(ExperimentalTime::class) -fun main() { +fun main() = with(Float64ParallelLinearSpace) { val random = Random(12224) val dim = 1000 //creating invertible matrix - val matrix1 = Double.algebra.linearSpace.buildMatrix(dim, dim) { i, j -> + val matrix1 = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } - val matrix2 = Double.algebra.linearSpace.buildMatrix(dim, dim) { i, j -> + val matrix2 = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } val time = measureTime { - with(Double.algebra.linearSpace) { - repeat(10) { - matrix1 dot matrix2 - } + repeat(30) { + matrix1 dot matrix2 } } diff --git a/examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt b/examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt new file mode 100644 index 000000000..dcd380ea6 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt @@ -0,0 +1,27 @@ +/* + * Copyright 2018-2024 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.linear + +import kotlin.random.Random +import kotlin.time.measureTime + +fun main(): Unit = with(Float64LinearSpace) { + val random = Random(1224) + val dim = 500 + + //creating invertible matrix + val u = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = buildMatrix(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = l dot u + + val time = measureTime { + repeat(20) { + lupSolver().inverse(matrix) + } + } + + println(time) +} \ No newline at end of file diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt b/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt index f8a8f1d0b..9dfb0fdc9 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt @@ -7,7 +7,7 @@ package space.kscience.kmath.series import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.Float64Buffer import space.kscience.kmath.structures.asBuffer import space.kscience.kmath.structures.toDoubleArray import space.kscience.plotly.* @@ -21,7 +21,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) { val arrayOfRandoms = DoubleArray(20) { random.nextDouble() } - val series1: DoubleBuffer = arrayOfRandoms.asBuffer() + val series1: Float64Buffer = arrayOfRandoms.asBuffer() val series2: Series = series1.moveBy(3) val res = series2 - series1 @@ -42,7 +42,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) { Plotly.plot { series("series1", series1) series("series2", series2) - series("dif", res){ + series("dif", res) { mode = ScatterMode.lines line.color("magenta") } diff --git a/gradle/wrapper/gradle-wrapper.properties b/gradle/wrapper/gradle-wrapper.properties index e411586a5..17655d0ef 100644 --- a/gradle/wrapper/gradle-wrapper.properties +++ b/gradle/wrapper/gradle-wrapper.properties @@ -1,5 +1,5 @@ distributionBase=GRADLE_USER_HOME distributionPath=wrapper/dists -distributionUrl=https\://services.gradle.org/distributions/gradle-8.4-bin.zip +distributionUrl=https\://services.gradle.org/distributions/gradle-8.6-bin.zip zipStoreBase=GRADLE_USER_HOME zipStorePath=wrapper/dists diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt index 800f27546..f1290d727 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt @@ -54,11 +54,12 @@ public object Float64LinearSpace : LinearSpace { require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } val rows = this@dot.rows.map { it.linearize() } val columns = other.columns.map { it.linearize() } + val indices = 0 until this.rowNum return buildMatrix(rowNum, other.colNum) { i, j -> val r = rows[i] val c = columns[j] var res = 0.0 - for (l in r.indices) { + for (l in indices) { res += r[l] * c[l] } res @@ -69,10 +70,11 @@ public object Float64LinearSpace : LinearSpace { override fun Matrix.dot(vector: Point): Float64Buffer { require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" } val rows = this@dot.rows.map { it.linearize() } + val indices = 0 until this.rowNum return Float64Buffer(rowNum) { i -> val r = rows[i] var res = 0.0 - for (j in r.indices) { + for (j in indices) { res += r[j] * vector[j] } res diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt index 14828f2df..7738fdaa3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt @@ -11,6 +11,7 @@ import space.kscience.attributes.Attributes import space.kscience.attributes.PolymorphicAttribute import space.kscience.attributes.safeTypeOf import space.kscience.kmath.UnstableKMathAPI +import space.kscience.kmath.nd.* import space.kscience.kmath.operations.* import space.kscience.kmath.structures.* @@ -150,7 +151,14 @@ public fun > LinearSpace>.lup( for (row in col + 1 until m) lu[row, col] /= luDiag } - return GenericLupDecomposition(this@lup, lu.toStructure2D(), pivot.asBuffer(), even) + val shape = ShapeND(rowNum, colNum) + + val structure2D = BufferND( + RowStrides(ShapeND(rowNum, colNum)), + lu + ).as2D() + + return GenericLupDecomposition(this@lup, structure2D, pivot.asBuffer(), even) } } @@ -166,7 +174,7 @@ internal fun LinearSpace>.solve( ): Matrix { require(matrix.rowNum == lup.l.rowNum) { "Matrix dimension mismatch. Expected ${lup.l.rowNum}, but got ${matrix.colNum}" } - BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory).run { + with(BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory)) { elementAlgebra { // Apply permutations to b val bp = create { _, _ -> zero } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt index 327eb2881..6c260db9f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt @@ -77,13 +77,11 @@ private inline fun > BufferAlgebra.zipInline( return elementBufferFactory(l.size) { elementAlgebra.block(l[it], r[it]) } } -public fun BufferAlgebra.buffer(size: Int, initializer: (Int) -> T): MutableBuffer { - return elementBufferFactory(size, initializer) -} +public fun BufferAlgebra.buffer(size: Int, initializer: (Int) -> T): MutableBuffer = + elementBufferFactory(size, initializer) -public fun A.buffer(initializer: (Int) -> T): MutableBuffer where A : BufferAlgebra, A : WithSize { - return elementBufferFactory(size, initializer) -} +public fun A.buffer(initializer: (Int) -> T): MutableBuffer where A : BufferAlgebra, A : WithSize = + elementBufferFactory(size, initializer) public fun > BufferAlgebra.sin(arg: Buffer): Buffer = mapInline(arg) { sin(it) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index 71cab61e0..4a0ac9416 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -42,6 +42,8 @@ public inline fun BufferFactory(): BufferFactory = BufferFactory( */ public interface MutableBufferFactory : BufferFactory { override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer + + public companion object } /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt index e03ebcade..5c8541f22 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt @@ -6,12 +6,7 @@ package space.kscience.kmath.structures import space.kscience.attributes.SafeType -import space.kscience.kmath.nd.BufferND -import space.kscience.kmath.nd.ShapeND import space.kscience.kmath.nd.Structure2D -import space.kscience.kmath.nd.as2D -import kotlin.collections.component1 -import kotlin.collections.component2 /** * A context that allows to operate on a [MutableBuffer] as on 2d array @@ -32,10 +27,10 @@ internal class BufferAccessor2D( fun create(mat: Structure2D): MutableBuffer = create { i, j -> mat[i, j] } - //TODO optimize wrapper - fun MutableBuffer.toStructure2D(): Structure2D = BufferND( - type, ShapeND(rowNum, colNum) - ) { (i, j) -> get(i, j) }.as2D() +// //TODO optimize wrapper +// fun MutableBuffer.toStructure2D(): Structure2D = BufferND( +// type, ShapeND(rowNum, colNum) +// ) { (i, j) -> get(i, j) }.as2D() inner class Row(val buffer: MutableBuffer, val rowIndex: Int) : MutableBuffer { override val type: SafeType get() = buffer.type diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt index a64abfc3d..ca9e8ef9f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt @@ -39,6 +39,7 @@ public value class Float64Buffer(public val array: DoubleArray) : PrimitiveBuffe } } +@Deprecated("Use Float64Buffer instead", ReplaceWith("Float64Buffer")) public typealias DoubleBuffer = Float64Buffer /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt index 44d2ecb5c..0acfd13a1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt @@ -77,7 +77,6 @@ public inline fun MutableBuffer( size: Int, initializer: (Int) -> T, ): MutableBuffer = when (type.kType) { - typeOf() -> TODO() typeOf() -> Int8Buffer(size) { initializer(it) as Int8 } as MutableBuffer typeOf() -> MutableBuffer.short(size) { initializer(it) as Int16 } as MutableBuffer typeOf() -> MutableBuffer.int(size) { initializer(it) as Int32 } as MutableBuffer diff --git a/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt new file mode 100644 index 000000000..9e05fded1 --- /dev/null +++ b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt @@ -0,0 +1,125 @@ +/* + * Copyright 2018-2024 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.linear + +import space.kscience.kmath.PerformancePitfall +import space.kscience.kmath.nd.* +import space.kscience.kmath.operations.Float64BufferOps +import space.kscience.kmath.operations.Float64Field +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.Float64Buffer +import space.kscience.kmath.structures.asBuffer +import java.util.stream.IntStream + + +public object Float64ParallelLinearSpace : LinearSpace { + + override val elementAlgebra: Float64Field = Float64Field + + + override fun buildMatrix( + rows: Int, + columns: Int, + initializer: Float64Field.(i: Int, j: Int) -> Double, + ): Matrix { + val shape = ShapeND(rows, columns) + val indexer = BufferAlgebraND.defaultIndexerBuilder(shape) + + val buffer = IntStream.range(0, indexer.linearSize).parallel().mapToDouble { offset -> + val (i, j) = indexer.index(offset) + elementAlgebra.initializer(i, j) + }.toArray().asBuffer() + + return MutableBufferND( + indexer, + buffer + ).as2D() + } + + override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Float64Buffer = + IntStream.range(0, size).parallel().mapToDouble{ Float64Field.initializer(it) }.toArray().asBuffer() + + override fun Matrix.unaryMinus(): Matrix = Floa64FieldOpsND { + asND().map { -it }.as2D() + } + + override fun Matrix.plus(other: Matrix): Matrix = Floa64FieldOpsND { + require(shape.contentEquals(other.shape)) { "Shape mismatch on Matrix::plus. Expected $shape but found ${other.shape}" } + asND().plus(other.asND()).as2D() + } + + override fun Matrix.minus(other: Matrix): Matrix = Floa64FieldOpsND { + require(shape.contentEquals(other.shape)) { "Shape mismatch on Matrix::minus. Expected $shape but found ${other.shape}" } + asND().minus(other.asND()).as2D() + } + + // Create a continuous in-memory representation of this vector for better memory layout handling + private fun Buffer.linearize() = if (this is Float64Buffer) { + this.array + } else { + DoubleArray(size) { get(it) } + } + + @OptIn(PerformancePitfall::class) + override fun Matrix.dot(other: Matrix): Matrix { + require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } + val rows = this@dot.rows.map { it.linearize() } + val columns = other.columns.map { it.linearize() } + val indices = 0 until this.rowNum + return buildMatrix(rowNum, other.colNum) { i, j -> + val r = rows[i] + val c = columns[j] + var res = 0.0 + for (l in indices) { + res += r[l] * c[l] + } + res + } + } + + @OptIn(PerformancePitfall::class) + override fun Matrix.dot(vector: Point): Float64Buffer { + require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" } + val rows = this@dot.rows.map { it.linearize() } + val indices = 0 until this.rowNum + return Float64Buffer(rowNum) { i -> + val r = rows[i] + var res = 0.0 + for (j in indices) { + res += r[j] * vector[j] + } + res + } + + } + + override fun Matrix.times(value: Double): Matrix = Floa64FieldOpsND { + asND().map { it * value }.as2D() + } + + public override fun Point.plus(other: Point): Float64Buffer = Float64BufferOps.run { + this@plus + other + } + + public override fun Point.minus(other: Point): Float64Buffer = Float64BufferOps.run { + this@minus - other + } + + public override fun Point.times(value: Double): Float64Buffer = Float64BufferOps.run { + scale(this@times, value) + } + + public operator fun Point.div(value: Double): Float64Buffer = Float64BufferOps.run { + scale(this@div, 1.0 / value) + } + + public override fun Double.times(v: Point): Float64Buffer = v * this + + +} + +public val Float64LinearSpace.parallel: Float64ParallelLinearSpace get() = Float64ParallelLinearSpace diff --git a/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt new file mode 100644 index 000000000..c57dba41e --- /dev/null +++ b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt @@ -0,0 +1,58 @@ +/* + * Copyright 2018-2024 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.structures + +import space.kscience.attributes.SafeType +import java.util.stream.Collectors +import java.util.stream.IntStream +import kotlin.reflect.typeOf + +/** + * A Java stream-based parallel version of [MutableBuffer]. + * There is no parallelization for [Int8], [Int16] and [Float32] types. + * They are processed sequentially. + */ +@Suppress("UNCHECKED_CAST") +public fun MutableBuffer.Companion.parallel( + type: SafeType, + size: Int, + initializer: (Int) -> T, +): MutableBuffer = when (type.kType) { + typeOf() -> Int8Buffer(size) { initializer(it) as Int8 } as MutableBuffer + typeOf() -> Int16Buffer(size) { initializer(it) as Int16 } as MutableBuffer + typeOf() -> IntStream.range(0, size).parallel().map { initializer(it) as Int32 }.toArray() + .asBuffer() as MutableBuffer + + typeOf() -> IntStream.range(0, size).parallel().mapToLong { initializer(it) as Int64 }.toArray() + .asBuffer() as MutableBuffer + + typeOf() -> Float32Buffer(size) { initializer(it) as Float } as MutableBuffer + typeOf() -> IntStream.range(0, size).parallel().mapToDouble { initializer(it) as Float64 }.toArray() + .asBuffer() as MutableBuffer + //TODO add unsigned types + else -> IntStream.range(0, size).parallel().mapToObj { initializer(it) }.collect(Collectors.toList()).asMutableBuffer(type) +} + +public class ParallelBufferFactory(override val type: SafeType) : MutableBufferFactory { + override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer { + TODO("Not yet implemented") + } + +} + +/** + * A Java stream-based parallel alternative to [MutableBufferFactory]. + * See [MutableBuffer.Companion.parallel] for details. + */ +public fun MutableBufferFactory.Companion.parallel( + type: SafeType, +): MutableBufferFactory = object : MutableBufferFactory { + + override val type: SafeType = type + + override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer = MutableBuffer.parallel(type, size, builder) + +} \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/ValueAndErrorField.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/ValueAndErrorField.kt index a79abff9f..96acf559a 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/ValueAndErrorField.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/ValueAndErrorField.kt @@ -59,7 +59,7 @@ public object ValueAndErrorField : Field { ValueAndError(a.value * value, a.dispersion * value.pow(2)) - private class ValueAndErrorBuffer(val values: DoubleBuffer, val ds: DoubleBuffer) : MutableBuffer { + private class ValueAndErrorBuffer(val values: Float64Buffer, val ds: Float64Buffer) : MutableBuffer { init { require(values.size == ds.size) }