benchmarks/build.gradle

@@ -7,6 +7,7 @@ plugins {
 dependencies {
     compile project(":kmath-core")
     compile project(":kmath-coroutines")
+    compile project(":kmath-commons")
     //jmh project(':kmath-core')
 }
 

benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt

@@ -0,0 +1,48 @@
+package scientifik.kmath.linear
+
+import org.apache.commons.math3.linear.Array2DRowRealMatrix
+import kotlin.random.Random
+import kotlin.system.measureTimeMillis
+
+fun main() {
+    val random = Random(12224)
+    val dim = 100
+    //creating invertible matrix
+    val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+    val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
+    val matrix = l dot u
+
+    val n = 500 // iterations
+
+    val solver = LUSolver.real
+
+    repeat(50) {
+        val res = solver.inverse(matrix)
+    }
+
+    val inverseTime = measureTimeMillis {
+        repeat(n) {
+            val res = solver.inverse(matrix)
+        }
+    }
+
+    println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis")
+
+    //commons
+
+    val cmSolver = CMSolver
+
+    val commonsMatrix = Array2DRowRealMatrix(dim, dim)
+    matrix.elements().forEach { (index, value) -> commonsMatrix.setEntry(index[0], index[1], value) }
+
+    val commonsTime = measureTimeMillis {
+        val cm = matrix.toCM()
+        repeat(n) {
+            //overhead on coversion could be mitigated
+            val res = cmSolver.inverse(cm)
+        }
+    }
+
+
+    println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
+}

kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt

@@ -0,0 +1,83 @@
+package scientifik.kmath.linear
+
+import org.apache.commons.math3.linear.*
+import org.apache.commons.math3.linear.RealMatrix
+import org.apache.commons.math3.linear.RealVector
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.structures.DoubleBuffer
+
+inline class CMMatrix(val origin: RealMatrix) : Matrix<Double> {
+    override val rowNum: Int get() = origin.rowDimension
+    override val colNum: Int get() = origin.columnDimension
+
+    override val features: Set<MatrixFeature> get() = emptySet()
+
+    override fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
+}
+
+fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {
+    this
+} else {
+    //TODO add feature analysis
+    val array = Array(rowNum) { i -> DoubleArray(colNum) { j -> get(i, j) } }
+    CMMatrix(Array2DRowRealMatrix(array))
+}
+
+fun RealMatrix.toMatrix() = CMMatrix(this)
+
+inline class CMVector(val origin: RealVector) : Point<Double> {
+    override val size: Int get() = origin.dimension
+
+    override fun get(index: Int): Double = origin.getEntry(index)
+
+    override fun iterator(): Iterator<Double> = origin.toArray().iterator()
+}
+
+fun Point<Double>.toCM(): CMVector = if (this is CMVector) {
+    this
+} else {
+    val array = DoubleArray(size) { this[it] }
+    CMVector(ArrayRealVector(array))
+}
+
+fun RealVector.toPoint() = DoubleBuffer(this.toArray())
+
+object CMMatrixContext : MatrixContext<Double, RealField> {
+    override val elementContext: RealField get() = RealField
+
+    override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): Matrix<Double> {
+        val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
+        return CMMatrix(Array2DRowRealMatrix(array))
+    }
+
+    override fun point(size: Int, initializer: (Int) -> Double): Point<Double> {
+        val array = DoubleArray(size, initializer)
+        return CMVector(ArrayRealVector(array))
+    }
+
+    override fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = this.toCM().dot(other.toCM())
+}
+
+operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin))
+operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.subtract(other.origin))
+
+infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(this.origin.multiply(other.origin))
+
+object CMSolver : LinearSolver<Double, RealField> {
+
+    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
+        val decomposition = LUDecomposition(a.toCM().origin)
+        return decomposition.solver.solve(b.toCM().origin).toMatrix()
+    }
+
+    override fun solve(a: Matrix<Double>, b: Point<Double>): Point<Double> {
+        val decomposition = LUDecomposition(a.toCM().origin)
+        return decomposition.solver.solve(b.toCM().origin).toPoint()
+    }
+
+    override fun inverse(a: Matrix<Double>): Matrix<Double> {
+        val decomposition = LUDecomposition(a.toCM().origin)
+        return decomposition.solver.inverse.toMatrix()
+    }
+
+}

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt

@@ -2,13 +2,16 @@ package scientifik.kmath.linear
 
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.Ring
-import scientifik.kmath.structures.*
+import scientifik.kmath.structures.MutableBuffer
 import scientifik.kmath.structures.MutableBuffer.Companion.boxing
+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>,
-    private val lu: NDStructure<T>,
+    internal val lu: NDStructure<T>,
     val pivot: IntArray,
     private val even: Boolean
 ) : DeterminantFeature<T> {
@@ -62,19 +65,12 @@ class LUPDecomposition<T : Comparable<T>>(
 
 
 class LUSolver<T : Comparable<T>, F : Field<T>>(
-    override val context: MatrixContext<T, F>,
+    val context: MatrixContext<T, F>,
     val bufferFactory: MutableBufferFactory<T> = ::boxing,
     val singularityCheck: (T) -> Boolean
 ) : LinearSolver<T, F> {
 
 
-    /**
-     * In-place transformation for [MutableNDStructure], using given transformation for each element
-     */
-    private operator fun <T> MutableNDStructure<T>.set(i: Int, j: Int, value: T) {
-        this[intArrayOf(i, j)] = value
-    }
-
     private fun abs(value: T) =
         if (value > context.elementContext.zero) value else with(context.elementContext) { -value }
 
@@ -180,19 +176,19 @@ class LUSolver<T : Comparable<T>, F : Field<T>>(
                 for (col in 0 until a.rowNum) {
                     for (i in col + 1 until a.rowNum) {
                         for (j in 0 until b.colNum) {
-                            bp[i, j] -= bp[col, j] * l[i, col]
+                            bp[i, j] -= bp[col, j] * lu[i, col]
                         }
                     }
                 }
 
                 // Solve UX = Y
                 for (col in a.rowNum - 1 downTo 0) {
-                    for(j in 0 until b.colNum){
+                    for (j in 0 until b.colNum) {
-                        bp[col,j]/= u[col,col]
+                        bp[col, j] /= lu[col, col]
                     }
                     for (i in 0 until col) {
                         for (j in 0 until b.colNum) {
-                            bp[i, j] -= bp[col, j]  * u[i, col]
+                            bp[i, j] -= bp[col, j] * lu[i, col]
                         }
                     }
                 }
@@ -202,6 +198,8 @@ class LUSolver<T : Comparable<T>, F : Field<T>>(
         }
     }
 
+    override fun inverse(a: Matrix<T>): Matrix<T> = solve(a, context.one(a.rowNum, a.colNum))
+
     companion object {
         val real = LUSolver(MatrixContext.real, MutableBuffer.Companion::auto) { it < 1e-11 }
     }

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt

@@ -12,12 +12,9 @@ import scientifik.kmath.structures.asSequence
  * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
  */
 interface LinearSolver<T : Any, R : Ring<T>>  {
-    val context: MatrixContext<T,R>
-
-
     fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
     fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.toMatrix()).toVector()
-    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, context.one(a.rowNum, a.colNum))
+    fun inverse(a: Matrix<T>): Matrix<T>
 }
 
 /**

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt

@@ -200,3 +200,5 @@ fun <T : Any, R : Ring<T>> Matrix<T>.transpose(): Matrix<T> {
         setOf(TransposedFeature(this))
     ) { i, j -> get(j, i) }
 }
+
+infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }