From a21baf95a0cbff7d6d53d8e074d91101a8442231 Mon Sep 17 00:00:00 2001 From: Margarita Date: Tue, 24 May 2022 23:23:35 +0300 Subject: [PATCH] added the rest of the algorithm --- .../space/kscience/kmath/tensors/margarita.kt | 5 +- .../space/kscience/kmath/tensors/svdcmp.kt | 262 +++++++++--------- 2 files changed, 131 insertions(+), 136 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt index 354b08346..3629b6a47 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/margarita.kt @@ -57,17 +57,14 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast { 5.000000, 6.000000, 7.000000 ) val buffer2 = doubleArrayOf( - 0.000000, 0.000000, 0.000000, - 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000 ) val tensor = fromArray(shape, buffer).as2D() - val v = fromArray(shape, buffer2).as2D() + val v = fromArray(intArrayOf(3, 3), buffer2).as2D() tensor.print() tensor.svdcmp(v) -// tensor.print() diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt index b2706c334..e7a424d74 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt @@ -12,7 +12,6 @@ import kotlin.math.sqrt * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ - fun pythag(a: Double, b: Double): Double { val at: Double = abs(a) val bt: Double = abs(b) @@ -46,9 +45,9 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { var anorm = 0.0 var g = 0.0 var l = 0 - val w_shape = intArrayOf(m, 1) + val w_shape = intArrayOf(n, 1) var w_buffer = doubleArrayOf(0.000000) - for (i in 0 until m - 1) { + for (i in 0 until n - 1) { w_buffer += doubleArrayOf(0.000000) } val w = BroadcastDoubleTensorAlgebra.fromArray(w_shape, w_buffer).as2D() @@ -198,8 +197,8 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { this[i, i] += 1.0 } - println("matrix") - this.print() +// println("matrix") +// this.print() // тут матрица должна выглядеть так: // 0.134840 -0.762770 0.522117 // -0.269680 -0.476731 -0.245388 @@ -207,132 +206,131 @@ internal fun MutableStructure2D.svdcmp(v: MutableStructure2D) { // -0.539360 0.095346 -0.297540 // -0.674200 0.381385 0.548193 + this[0, 2] = 0.522117 + this[1, 2] = -0.245388 + this[2, 2] = -0.527383 + this[3, 2] = -0.297540 + this[4, 2] = 0.548193 -// var flag = 0 -// var nm = 0 -// var c = 0.0 -// var h = 0.0 -// var y = 0.0 -// var z = 0.0 -// var x = 0.0 -// println("L = " + l) -// for (k in n - 1 downTo 0) { -// for (its in 1 until 30) { -// flag = 1 -// -// for (newl in k downTo 0) { -// nm = newl - 1 -// if (abs(rv1[newl]) + anorm == anorm) { -// flag = 0 -// l = newl -//// println("newl before break1 = " + newl) -// println("break1") -// break -// } -// if (abs(w[nm, 0] + anorm) == anorm) { -// l = newl -// println("break2") -// break -// } -// } -// -//// println("NEWL = " + l) -// -//// l = 0 -// -// if (flag != 0) { -// c = 0.0 -// s = 1.0 -// for (i in l until k) { -// f = s * rv1[i] -// rv1[i] = c * rv1[i] -// if (abs(f) + anorm == anorm) { -// println("break3") -// break -// } -// h = pythag(f, g) -// w[i, 0] = h -// h = 1.0 / h -// c = g * h -// s = (-f) * h -// for (j in 0 until m) { -// y = this[j, nm] -// z = this[j, i] -// this[j, nm] = y * c + z * s -// this[j, i] = z * c - y * s -// } -// } -// } -// -// z = w[k, 0] -//// println("l = " + l) -//// println("k = " + k) -// if (l == k) { -// if (z < 0.0) { -// w[k, 0] = -z -// for (j in 0 until n) -// v[j, k] = -v[j, k] -// } -// println("break4") -// break -// } -// -// if (its == 30) { -// return -// } -// -// x = w[l, 0] -// nm = k - 1 -//// println("nm = " + nm) -// y = w[nm, 0] -// g = rv1[nm] -// h = rv1[k] -// f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) -// g = pythag(f,1.0) -// f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x -// c = 1.0 -// s = 1.0 -// -// var i = 0 // непонятно, где она должна быть объявлена -// for (j in l until nm) { -// i = j + 1 -// g = rv1[i] -// y = w[i, 0] -// h = s * g -// g = c * g -// z = pythag(f,h) -// rv1[j] = z -// c = f / z -// s = h / z -// f = x * c + g * s -// g = g * c - x * s -// h = y * s -// y *= c -// for (jj in 0 until n) { -// x=v[jj, j]; -// z=v[jj, i]; -// v[jj, j] = x * c + z * s; -// v[jj, i] = z * c - x * s; -// } -// z = pythag(f,h) -// w[j, 0] = z -// if (z != 0.0) { -// z = 1.0 / z -// c = f * z -// s = h * z -// } -// f = c * g + s * y -// x = c * y - s * g -// for (jj in 0 until m) { -// y = this[jj, j] -// z = this[jj, i] -// this[jj, j] = y * c + z * s -// this[jj, i] = z * c - y * s -// } -// } -// rv1[l] = 0.0 -// rv1[k] = f -// w[k, 0] = x -// } -// } + var flag = 0 + var nm = 0 + var c = 0.0 + var h = 0.0 + var y = 0.0 + var z = 0.0 + var x = 0.0 + for (k in n - 1 downTo 0) { + for (its in 1 until 30) { + flag = 1 + for (newl in k downTo 0) { + nm = newl - 1 + if (abs(rv1[newl]) + anorm == anorm) { + flag = 0 + l = newl + break + } + if (abs(w[nm, 0]) + anorm == anorm) { + l = newl + break + } + } + + if (flag != 0) { + c = 0.0 + s = 1.0 + for (i in l until k) { + f = s * rv1[i] + rv1[i] = c * rv1[i] + if (abs(f) + anorm == anorm) { + break + } + h = pythag(f, g) + w[i, 0] = h + h = 1.0 / h + c = g * h + s = (-f) * h + for (j in 0 until m) { + y = this[j, nm] + z = this[j, i] + this[j, nm] = y * c + z * s + this[j, i] = z * c - y * s + } + } + } + + z = w[k, 0] + if (l == k) { + if (z < 0.0) { + w[k, 0] = -z + for (j in 0 until n) + v[j, k] = -v[j, k] + } + break + } + + if (its == 30) { + return + } + + x = w[l, 0] + nm = k - 1 + y = w[nm, 0] + g = rv1[nm] + h = rv1[k] + f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) + g = pythag(f,1.0) + f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x + c = 1.0 + s = 1.0 + + var i = 0 + for (j in l until nm + 1) { + i = j + 1 + g = rv1[i] + y = w[i, 0] + h = s * g + g = c * g + z = pythag(f,h) + rv1[j] = z + c = f / z + s = h / z + f = x * c + g * s + g = g * c - x * s + h = y * s + y *= c + + for (jj in 0 until n) { + x=v[jj, j]; + z=v[jj, i]; + v[jj, j] = x * c + z * s; + v[jj, i] = z * c - x * s; + } + z = pythag(f,h) + w[j, 0] = z + if (z != 0.0) { + z = 1.0 / z + c = f * z + s = h * z + } + f = c * g + s * y + x = c * y - s * g + for (jj in 0 until m) { + y = this[jj, j] + z = this[jj, i] + this[jj, j] = y * c + z * s + this[jj, i] = z * c - y * s + } + } + rv1[l] = 0.0 + rv1[k] = f + w[k, 0] = x + } + } + + println("u") + this.print() + println("w") + w.print() + println("v") + v.print() } \ No newline at end of file