.gitignore

@@ -19,3 +19,4 @@ out/
 
 !/.idea/copyright/
 !/.idea/scopes/
+/kotlin-js-store/yarn.lock

benchmarks/build.gradle.kts

@@ -52,6 +52,8 @@ kotlin {
                 implementation(project(":kmath-viktor"))
                 implementation(project(":kmath-jafama"))
                 implementation(project(":kmath-multik"))
+                implementation(projects.kmath.kmathTensorflow)
+                implementation("org.tensorflow:tensorflow-core-platform:0.4.0")
                 implementation("org.nd4j:nd4j-native:1.0.0-M1")
                 //    uncomment if your system supports AVX2
                 //    val os = System.getProperty("os.name")
@@ -122,6 +124,11 @@ benchmark {
         include("JafamaBenchmark")
     }
 
+    configurations.register("tensorAlgebra") {
+        commonConfiguration()
+        include("TensorAlgebraBenchmark")
+    }
+
     configurations.register("viktor") {
         commonConfiguration()
         include("ViktorBenchmark")

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt

@@ -17,7 +17,9 @@ import space.kscience.kmath.multik.multikAlgebra
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.tensorflow.produceWithTF
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.tensorAlgebra
 import kotlin.random.Random
 
 @State(Scope.Benchmark)
@@ -34,9 +36,6 @@ internal class DotBenchmark {
             random.nextDouble()
         }
 
-        val tensor1 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12224)
-        val tensor2 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12225)
-
         val cmMatrix1 = CMLinearSpace { matrix1.toCM() }
         val cmMatrix2 = CMLinearSpace { matrix2.toCM() }
 
@@ -44,6 +43,16 @@ internal class DotBenchmark {
         val ejmlMatrix2 = EjmlLinearSpaceDDRM { matrix2.toEjml() }
     }
 
+
+    @Benchmark
+    fun tfDot(blackhole: Blackhole) {
+        blackhole.consume(
+            DoubleField.produceWithTF {
+                matrix1 dot matrix1
+            }
+        )
+    }
+
     @Benchmark
     fun cmDotWithConversion(blackhole: Blackhole) = CMLinearSpace {
         blackhole.consume(matrix1 dot matrix2)
@@ -64,13 +73,13 @@ internal class DotBenchmark {
         blackhole.consume(matrix1 dot matrix2)
     }
 
-//    @Benchmark
+    @Benchmark
-//    fun tensorDot(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+    fun tensorDot(blackhole: Blackhole) = with(DoubleField.tensorAlgebra) {
-//        blackhole.consume(matrix1 dot matrix2)
+        blackhole.consume(matrix1 dot matrix2)
-//    }
+    }
 
     @Benchmark
-    fun multikDot(blackhole: Blackhole) = with(Double.multikAlgebra) {
+    fun multikDot(blackhole: Blackhole) = with(DoubleField.multikAlgebra) {
         blackhole.consume(matrix1 dot matrix2)
     }
 
@@ -86,6 +95,6 @@ internal class DotBenchmark {
 
     @Benchmark
     fun doubleTensorDot(blackhole: Blackhole) = DoubleTensorAlgebra.invoke {
-        blackhole.consume(tensor1 dot tensor2)
+        blackhole.consume(matrix1 dot matrix2)
     }
 }

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt

@@ -0,0 +1,37 @@
+/*
+ * Copyright 2018-2021 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.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.linear.symmetric
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.tensors.core.tensorAlgebra
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class TensorAlgebraBenchmark {
+    companion object {
+        private val random = Random(12224)
+        private const val dim = 30
+
+        private val matrix = DoubleField.linearSpace.matrix(dim, dim).symmetric { _, _ -> random.nextDouble() }
+    }
+
+    @Benchmark
+    fun tensorSymEigSvd(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        blackhole.consume(matrix.symEigSvd(1e-10))
+    }
+
+    @Benchmark
+    fun tensorSymEigJacobi(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        blackhole.consume(matrix.symEigJacobi(50, 1e-10))
+    }
+}

build.gradle.kts

@@ -10,7 +10,7 @@ allprojects {
     }
 
     group = "space.kscience"
-    version = "0.3.0-dev-18"
+    version = "0.3.0-dev-19"
 }
 
 subprojects {

examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt

@@ -5,8 +5,8 @@
 
 package space.kscience.kmath.functions
 
-import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
+import space.kscience.kmath.interpolation.splineInterpolator
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
@@ -28,7 +28,7 @@ fun main() {
     val xs = 0.0..100.0 step 0.5
     val ys = xs.map(function)
 
-    val polynomial: PiecewisePolynomial<Double> = SplineInterpolator.double.interpolatePolynomials(xs, ys)
+    val polynomial: PiecewisePolynomial<Double> = DoubleField.splineInterpolator.interpolatePolynomials(xs, ys)
 
     val polyFunction = polynomial.asFunction(DoubleField, 0.0)
 

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt

@@ -28,6 +28,8 @@ public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T
 /**
  * Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]
  * Requires [UnivariateIntegrationNodes] or [IntegrationRange] and [IntegrandMaxCalls]
+ *
+ * TODO use context receiver for algebra
  */
 @UnstableKMathAPI
 public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(
@@ -98,6 +100,7 @@ public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
     }
 }
 
+@Suppress("unused")
 @UnstableKMathAPI
 public inline val DoubleField.splineIntegrator: UnivariateIntegrator<Double>
     get() = DoubleSplineIntegrator

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt

@@ -9,6 +9,7 @@ package space.kscience.kmath.interpolation
 
 import space.kscience.kmath.data.XYColumnarData
 import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.asFunction
 import space.kscience.kmath.functions.value
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
@@ -59,3 +60,33 @@ public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
     return interpolatePolynomials(pointSet)
 }
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    x: Buffer<T>,
+    y: Buffer<T>,
+): (T) -> T? = interpolatePolynomials(x, y).asFunction(algebra)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: Map<T, T>,
+): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: List<Pair<T, T>>,
+): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
+
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    x: Buffer<T>,
+    y: Buffer<T>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(x, y).asFunction(algebra, defaultValue)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: Map<T, T>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: List<Pair<T, T>>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt

@@ -22,6 +22,7 @@ internal fun <T : Comparable<T>> insureSorted(points: XYColumnarData<*, T, *>) {
  * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
  */
 public class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
+
     @OptIn(UnstableKMathAPI::class)
     override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
         require(points.size > 0) { "Point array should not be empty" }
@@ -37,3 +38,6 @@ public class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T
         }
     }
 }
+
+public val <T : Comparable<T>> Field<T>.linearInterpolator: LinearInterpolator<T>
+    get() = LinearInterpolator(this)

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt

@@ -63,8 +63,8 @@ public class SplineInterpolator<T : Comparable<T>>(
                 //Shift coefficients to represent absolute polynomial instead of one with an offset
                 val polynomial = Polynomial(
                     a - b * x0 + c * x02 - d * x03,
-                    b - 2*c*x0 + 3*d*x02,
+                    b - 2 * c * x0 + 3 * d * x02,
-                    c - 3*d*x0,
+                    c - 3 * d * x0,
                     d
                 )
                 cOld = c
@@ -72,8 +72,12 @@ public class SplineInterpolator<T : Comparable<T>>(
             }
         }
     }
-
-    public companion object {
-        public val double: SplineInterpolator<Double> = SplineInterpolator(DoubleField, ::DoubleBuffer)
-    }
 }
+
+
+public fun <T : Comparable<T>> Field<T>.splineInterpolator(
+    bufferFactory: MutableBufferFactory<T>,
+): SplineInterpolator<T> = SplineInterpolator(this, bufferFactory)
+
+public val DoubleField.splineInterpolator: SplineInterpolator<Double>
+    get() = SplineInterpolator(this, ::DoubleBuffer)

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/LinearInterpolatorTest.kt

@@ -5,8 +5,6 @@
 
 package space.kscience.kmath.interpolation
 
-import space.kscience.kmath.functions.PiecewisePolynomial
-import space.kscience.kmath.functions.asFunction
 import space.kscience.kmath.operations.DoubleField
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -21,8 +19,8 @@ internal class LinearInterpolatorTest {
             3.0 to 4.0
         )
 
-        val polynomial: PiecewisePolynomial<Double> = LinearInterpolator(DoubleField).interpolatePolynomials(data)
+        //val polynomial: PiecewisePolynomial<Double> = DoubleField.linearInterpolator.interpolatePolynomials(data)
-        val function = polynomial.asFunction(DoubleField)
+        val function = DoubleField.linearInterpolator.interpolate(data)
         assertEquals(null, function(-1.0))
         assertEquals(0.5, function(0.5))
         assertEquals(2.0, function(1.5))

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt

@@ -5,8 +5,6 @@
 
 package space.kscience.kmath.interpolation
 
-import space.kscience.kmath.functions.PiecewisePolynomial
-import space.kscience.kmath.functions.asFunction
 import space.kscience.kmath.operations.DoubleField
 import kotlin.math.PI
 import kotlin.math.sin
@@ -21,9 +19,10 @@ internal class SplineInterpolatorTest {
             x to sin(x)
         }
 
-        val polynomial: PiecewisePolynomial<Double> = SplineInterpolator.double.interpolatePolynomials(data)
+        //val polynomial: PiecewisePolynomial<Double> = DoubleField.splineInterpolator.interpolatePolynomials(data)
+
+        val function = DoubleField.splineInterpolator.interpolate(data, Double.NaN)
 
-        val function = polynomial.asFunction(DoubleField, Double.NaN)
         assertEquals(Double.NaN, function(-1.0))
         assertEquals(sin(0.5), function(0.5), 0.1)
         assertEquals(sin(1.5), function(1.5), 0.1)

kmath-multik/src/test/kotlin/space/kscience/kmath/multik/MultikNDTest.kt

@@ -6,13 +6,38 @@
 package space.kscience.kmath.multik
 
 import org.junit.jupiter.api.Test
+import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.tensorAlgebra
+import kotlin.test.assertTrue
 
 internal class MultikNDTest {
     @Test
     fun basicAlgebra(): Unit = DoubleField.multikAlgebra{
         one(2,2) + 1.0
     }
+
+    @Test
+    fun dotResult(){
+        val dim = 100
+
+        val tensor1 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12224)
+        val tensor2 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12225)
+
+        val multikResult = with(DoubleField.multikAlgebra){
+            tensor1 dot tensor2
+        }
+
+        val defaultResult = with(DoubleField.tensorAlgebra){
+            tensor1 dot tensor2
+        }
+
+        assertTrue {
+            StructureND.contentEquals(multikResult, defaultResult)
+        }
+
+    }
 }

kmath-tensorflow/src/main/kotlin/space/kscience/kmath/tensorflow/TensorFlowAlgebra.kt

@@ -199,8 +199,9 @@ public abstract class TensorFlowAlgebra<T, TT : TNumber, A : Ring<T>> internal c
 
     override fun StructureND<T>.dot(other: StructureND<T>): TensorFlowOutput<T, TT> = operate(other) { l, r ->
         ops.linalg.matMul(
-            if (l.asTensor().shape().numDimensions() == 1) ops.expandDims(l, ops.constant(0)) else l,
+            if (l.shape().numDimensions() == 1) ops.expandDims(l, ops.constant(0)) else l,
-            if (r.asTensor().shape().numDimensions() == 1) ops.expandDims(r, ops.constant(-1)) else r)
+            if (r.shape().numDimensions() == 1) ops.expandDims(r, ops.constant(-1)) else r
+        )
     }
 
     override fun diagonalEmbedding(

kmath-tensorflow/src/test/kotlin/space/kscience/kmath/tensorflow/DoubleTensorFlowOps.kt

@@ -4,6 +4,8 @@ import org.junit.jupiter.api.Test
 import space.kscience.kmath.nd.get
 import space.kscience.kmath.nd.structureND
 import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.sum
 import kotlin.test.assertEquals
 
 class DoubleTensorFlowOps {
@@ -18,6 +20,18 @@ class DoubleTensorFlowOps {
         assertEquals(3.0, res[0, 0])
     }
 
+    @Test
+    fun dot(){
+        val dim = 1000
+
+        val tensor1 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12224)
+        val tensor2 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12225)
+
+        DoubleField.produceWithTF {
+            tensor1 dot tensor2
+        }.sum()
+    }
+
     @Test
     fun extensionOps(){
         val res = DoubleField.produceWithTF {

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

@@ -9,10 +9,7 @@
 package space.kscience.kmath.tensors.core
 
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.*
-import space.kscience.kmath.nd.StructureND
-import space.kscience.kmath.nd.as1D
-import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.indices
@@ -885,7 +882,7 @@ public open class DoubleTensorAlgebra :
         return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
     }
 
-    override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEig(epsilon = 1e-15)
+    override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEigJacobi(maxIteration = 50, epsilon = 1e-15)
 
     /**
      * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices,
@@ -895,7 +892,7 @@ public open class DoubleTensorAlgebra :
      * and when the cosine approaches 1 in the SVD algorithm.
      * @return a pair `eigenvalues to eigenvectors`.
      */
-    public fun StructureND<Double>.symEig(epsilon: Double): Pair<DoubleTensor, DoubleTensor> {
+    public fun StructureND<Double>.symEigSvd(epsilon: Double): Pair<DoubleTensor, DoubleTensor> {
         checkSymmetric(tensor, epsilon)
 
         fun MutableStructure2D<Double>.cleanSym(n: Int) {
@@ -922,6 +919,151 @@ public open class DoubleTensorAlgebra :
         return eig to v
     }
 
+    public fun StructureND<Double>.symEigJacobi(maxIteration: Int, epsilon: Double): Pair<DoubleTensor, DoubleTensor> {
+        checkSymmetric(tensor, epsilon)
+
+        val size = this.dimension
+        val eigenvectors = zeros(this.shape)
+        val eigenvalues = zeros(this.shape.sliceArray(0 until size - 1))
+
+        var eigenvalueStart = 0
+        var eigenvectorStart = 0
+        for (matrix in tensor.matrixSequence()) {
+            val matrix2D = matrix.as2D()
+            val (d, v) = matrix2D.jacobiHelper(maxIteration, epsilon)
+
+            for (i in 0 until matrix2D.rowNum) {
+                for (j in 0 until matrix2D.colNum) {
+                    eigenvectors.mutableBuffer.array()[eigenvectorStart + i * matrix2D.rowNum + j] = v[i, j]
+                }
+            }
+
+            for (i in 0 until matrix2D.rowNum) {
+                eigenvalues.mutableBuffer.array()[eigenvalueStart + i] = d[i]
+            }
+
+            eigenvalueStart += this.shape.last()
+            eigenvectorStart += this.shape.last() * this.shape.last()
+        }
+
+        return eigenvalues to eigenvectors
+    }
+
+    private fun MutableStructure2D<Double>.jacobiHelper(
+        maxIteration: Int,
+        epsilon: Double
+    ): Pair<Structure1D<Double>, Structure2D<Double>> {
+        val n = this.shape[0]
+        val A_ = this.copy()
+        val V = eye(n)
+        val D = DoubleTensor(intArrayOf(n), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D()
+        val B = DoubleTensor(intArrayOf(n), (0 until this.rowNum).map { this[it, it] }.toDoubleArray()).as1D()
+        val Z = zeros(intArrayOf(n)).as1D()
+
+        // assume that buffered tensor is square matrix
+        operator fun BufferedTensor<Double>.get(i: Int, j: Int): Double {
+            return this.mutableBuffer.array()[bufferStart + i * this.shape[0] + j]
+        }
+
+        operator fun BufferedTensor<Double>.set(i: Int, j: Int, value: Double) {
+            this.mutableBuffer.array()[bufferStart + i * this.shape[0] + j] = value
+        }
+
+        fun maxOffDiagonal(matrix: BufferedTensor<Double>): Double {
+            var maxOffDiagonalElement = 0.0
+            for (i in 0 until n - 1) {
+                for (j in i + 1 until n) {
+                    maxOffDiagonalElement = max(maxOffDiagonalElement, abs(matrix[i, j]))
+                }
+            }
+            return maxOffDiagonalElement
+        }
+
+        fun rotate(a: BufferedTensor<Double>, s: Double, tau: Double, i: Int, j: Int, k: Int, l: Int) {
+            val g = a[i, j]
+            val h = a[k, l]
+            a[i, j] = g - s * (h + g * tau)
+            a[k, l] = h + s * (g - h * tau)
+        }
+
+        fun jacobiIteration(
+            a: BufferedTensor<Double>,
+            v: BufferedTensor<Double>,
+            d: MutableStructure1D<Double>,
+            z: MutableStructure1D<Double>,
+        ) {
+            for (ip in 0 until n - 1) {
+                for (iq in ip + 1 until n) {
+                    val g = 100.0 * abs(a[ip, iq])
+
+                    if (g <= epsilon * abs(d[ip]) && g <= epsilon * abs(d[iq])) {
+                        a[ip, iq] = 0.0
+                        continue
+                    }
+
+                    var h = d[iq] - d[ip]
+                    val t = when {
+                        g <= epsilon * abs(h) -> (a[ip, iq]) / h
+                        else -> {
+                            val theta = 0.5 * h / (a[ip, iq])
+                            val denominator = abs(theta) + sqrt(1.0 + theta * theta)
+                            if (theta < 0.0) -1.0 / denominator else 1.0 / denominator
+                        }
+                    }
+
+                    val c = 1.0 / sqrt(1 + t * t)
+                    val s = t * c
+                    val tau = s / (1.0 + c)
+                    h = t * a[ip, iq]
+                    z[ip] -= h
+                    z[iq] += h
+                    d[ip] -= h
+                    d[iq] += h
+                    a[ip, iq] = 0.0
+
+                    for (j in 0 until ip) {
+                        rotate(a, s, tau, j, ip, j, iq)
+                    }
+                    for (j in (ip + 1) until iq) {
+                        rotate(a, s, tau, ip, j, j, iq)
+                    }
+                    for (j in (iq + 1) until n) {
+                        rotate(a, s, tau, ip, j, iq, j)
+                    }
+                    for (j in 0 until n) {
+                        rotate(v, s, tau, j, ip, j, iq)
+                    }
+                }
+            }
+        }
+
+        fun updateDiagonal(
+            d: MutableStructure1D<Double>,
+            z: MutableStructure1D<Double>,
+            b: MutableStructure1D<Double>,
+        ) {
+            for (ip in 0 until d.size) {
+                b[ip] += z[ip]
+                d[ip] = b[ip]
+                z[ip] = 0.0
+            }
+        }
+
+        var sm = maxOffDiagonal(A_)
+        for (iteration in 0 until maxIteration) {
+            if (sm < epsilon) {
+                break
+            }
+
+            jacobiIteration(A_, V, D, Z)
+            updateDiagonal(D, Z, B)
+            sm = maxOffDiagonal(A_)
+        }
+
+        // TODO sort eigenvalues
+        return D to V.as2D()
+    }
+
     /**
      * Computes the determinant of a square matrix input, or of each square matrix in a batched input
      * using LU factorization algorithm.
@@ -997,5 +1139,6 @@ public open class DoubleTensorAlgebra :
 }
 
 public val Double.Companion.tensorAlgebra: DoubleTensorAlgebra.Companion get() = DoubleTensorAlgebra
+public val DoubleField.tensorAlgebra: DoubleTensorAlgebra.Companion get() = DoubleTensorAlgebra