From c4afb69ba9a2b598b0428645593c8ad400c457f4 Mon Sep 17 00:00:00 2001
From: margarita0303 <lashina.2000@gmail.com>
Date: Tue, 12 Jul 2022 17:19:19 +0300
Subject: [PATCH] added svdGolabKahan

---
 .../kmath/tensors/core/DoubleTensorAlgebra.kt | 28 +++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
index e9dc34748..174b1dde2 100644
--- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
@@ -882,6 +882,34 @@ public open class DoubleTensorAlgebra :
         return Triple(uTensor.transpose(), sTensor, vTensor.transpose())
     }
 
+    public fun StructureND<Double>.svdGolabKahan(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
+        val size = tensor.dimension
+        val commonShape = tensor.shape.sliceArray(0 until size - 2)
+        val (n, m) = tensor.shape.sliceArray(size - 2 until size)
+        val uTensor = zeros(commonShape + intArrayOf(min(n, m), n))
+        val sTensor = zeros(commonShape + intArrayOf(min(n, m)))
+        val vTensor = zeros(commonShape + intArrayOf(min(n, m), m))
+
+        val matrices = tensor.matrices
+        val uTensors = uTensor.matrices
+        val sTensorVectors = sTensor.vectors
+        val vTensors = vTensor.matrices
+
+        for (index in matrices.indices) {
+            val matrix = matrices[index]
+            val matrixSize = matrix.shape.reduce { acc, i -> acc * i }
+            val curMatrix = DoubleTensor(
+                matrix.shape,
+                matrix.mutableBuffer.array()
+                    .slice(matrix.bufferStart until matrix.bufferStart + matrixSize)
+                    .toDoubleArray()
+            )
+            curMatrix.as2D().svdGolabKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D())
+        }
+
+        return Triple(uTensor.transpose(), sTensor, vTensor)
+    }
+
     override fun StructureND<Double>.symEig(): Pair<DoubleTensor, DoubleTensor> = symEigJacobi(maxIteration = 50, epsilon = 1e-15)
 
     /**