kmath-core/api/kmath-core.api

@@ -2645,7 +2645,7 @@ public abstract interface class space/kscience/kmath/tensors/LinearOpsTensorAlge
 	public abstract fun inv (Lspace/kscience/kmath/nd/MutableStructureND;)Lspace/kscience/kmath/nd/MutableStructureND;
 	public abstract fun lu (Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Pair;
 	public abstract fun luPivot (Lspace/kscience/kmath/nd/MutableStructureND;Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Triple;
-	public abstract fun qr (Lspace/kscience/kmath/nd/MutableStructureND;)Lspace/kscience/kmath/nd/MutableStructureND;
+	public abstract fun qr (Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Pair;
 	public abstract fun svd (Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Triple;
 	public abstract fun symEig (Lspace/kscience/kmath/nd/MutableStructureND;Z)Lkotlin/Pair;
 }
@@ -2794,8 +2794,8 @@ public final class space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebr
 	public fun lu (Lspace/kscience/kmath/tensors/core/DoubleTensor;)Lkotlin/Pair;
 	public synthetic fun luPivot (Lspace/kscience/kmath/nd/MutableStructureND;Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Triple;
 	public fun luPivot (Lspace/kscience/kmath/tensors/core/DoubleTensor;Lspace/kscience/kmath/tensors/core/IntTensor;)Lkotlin/Triple;
-	public synthetic fun qr (Lspace/kscience/kmath/nd/MutableStructureND;)Lspace/kscience/kmath/nd/MutableStructureND;
+	public synthetic fun qr (Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Pair;
-	public fun qr (Lspace/kscience/kmath/tensors/core/DoubleTensor;)Lspace/kscience/kmath/tensors/core/DoubleTensor;
+	public fun qr (Lspace/kscience/kmath/tensors/core/DoubleTensor;)Lkotlin/Pair;
 	public synthetic fun svd (Lspace/kscience/kmath/nd/MutableStructureND;)Lkotlin/Triple;
 	public fun svd (Lspace/kscience/kmath/tensors/core/DoubleTensor;)Lkotlin/Triple;
 	public synthetic fun symEig (Lspace/kscience/kmath/nd/MutableStructureND;Z)Lkotlin/Pair;

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

@@ -84,7 +84,7 @@ public class DoubleLinearOpsTensorAlgebra :
         val seq = matrixSequence().zip((qTensor.matrixSequence().zip(rTensor.matrixSequence())))
         for ((matrix, qr) in seq) {
             val (q, r) = qr
-            qrHelper(matrix.as2D(), q.as2D(), r.as2D())
+            qrHelper(matrix.asTensor(), q.asTensor(), r.as2D())
         }
         return Pair(qTensor, rTensor)
     }

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

@@ -187,42 +187,32 @@ internal inline fun luMatrixInv(
     }
 }
 
-internal inline fun MutableStructure1D<Double>.dot(other: MutableStructure1D<Double>): Double {
+internal inline fun DoubleLinearOpsTensorAlgebra.qrHelper(
-    var res = 0.0
+    matrix: DoubleTensor,
-    for (i in 0 until size) {
+    q: DoubleTensor,
-        res += this[i] * other[i]
-    }
-    return res
-}
-
-internal inline fun MutableStructure1D<Double>.l2Norm(): Double {
-    var squareSum = 0.0
-    for (i in 0 until size) {
-        squareSum += this[i] * this[i]
-    }
-    return sqrt(squareSum)
-}
-
-internal inline fun qrHelper(
-    matrix: MutableStructure2D<Double>,
-    q: MutableStructure2D<Double>,
     r: MutableStructure2D<Double>
 ) {
-    //todo check square
+    checkSquareMatrix(matrix.shape)
-    val n = matrix.colNum
+    val n = matrix.shape[0]
+    val qM = q.as2D()
+    val matrixT = matrix.transpose(0,1)
+    val qT = q.transpose(0,1)
+
     for (j in 0 until n) {
-        val v = matrix.columns[j]
+        val v = matrixT[j]
+        val vv = v.as1D()
         if (j > 0) {
             for (i in 0 until j) {
-                r[i, j] = q.columns[i].dot(matrix.columns[j])
+                r[i, j] = qT[i].dot(matrixT[j]).value()
                 for (k in 0 until n) {
-                    v[k] = v[k] - r[i, j] * q.columns[i][k]
+                    val qTi = qT[i].as1D()
+                    vv[k] = vv[k] - r[i, j] * qTi[k]
                 }
             }
         }
-        r[j, j] = v.l2Norm()
+        r[j, j] = DoubleAnalyticTensorAlgebra { v.dot(v).sqrt().value() }
         for (i in 0 until n) {
-            q[i, j] = v[i] / r[j, j]
+            qM[i, j] = vv[i] / r[j, j]
         }
     }
 }