teldufalsari/kmath
2
0
Fork 0
forked from kscience/kmath
Code Pull Requests 1 Activity

added comment about division by zero

Browse Source
This commit is contained in:
Margarita 2022-05-24 19:50:50 +03:00
parent 2fa39fff14
commit 86efe48217
1 changed files with 103 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

144
examples/src/main/kotlin/space/kscience/kmath/tensors/svdcmp.kt
View File

@ -164,26 +164,18 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
l = i l = i
} }
// тут все правильно считается // до этого момента все правильно считается
// дальше - нет
// println("w")
// w.print()
//
val eps = 0.000000001
// println("1.0 / w[2, 0] " + 1.0 / w[2, 0])
// println("w[2, 0] " + w[2, 0])
for (i in min(n, m) - 1 downTo 0) { for (i in min(n, m) - 1 downTo 0) {
l = i + 1 l = i + 1
g = w[i, 0] g = w[i, 0]
// println("w[i, 0] " + w[i, 0])
for (j in l until n) { for (j in l until n) {
this[i, j] = 0.0 this[i, j] = 0.0
} }
if (g != 0.0) { if (g != 0.0) {
// !!!!! вот тут деление на почти ноль
g = 1.0 / g g = 1.0 / g
// println("g " + g)
for (j in l until n) { for (j in l until n) {
s = 0.0 s = 0.0
for (k in l until m) { for (k in l until m) {
@ -204,15 +196,11 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
} }
} }
this[i, i] += 1.0 this[i, i] += 1.0
// println("matrix")
// this.print()
} }
println("matrix") println("matrix")
this.print() this.print()
// тут матрица должна выглядеть так:
// тут матрица должна выглядеть так:
// 0.134840 -0.762770 0.522117 // 0.134840 -0.762770 0.522117
// -0.269680 -0.476731 -0.245388 // -0.269680 -0.476731 -0.245388
// -0.404520 -0.190693 -0.527383 // -0.404520 -0.190693 -0.527383
@ -220,57 +208,131 @@ internal fun MutableStructure2D<Double>.svdcmp(v: MutableStructure2D<Double>) {
// -0.674200 0.381385 0.548193 // -0.674200 0.381385 0.548193
// var flag = 0 // var flag = 0
// var nm = 0 // var nm = 0
// var c = 0.0 // var c = 0.0
// var h = 0.0 // var h = 0.0
// var y = 0.0 // var y = 0.0
// var z = 0.0 // var z = 0.0
// var x = 0.0
// println("L = " + l)
// for (k in n - 1 downTo 0) { // for (k in n - 1 downTo 0) {
// for (its in 0 until 30) { // for (its in 1 until 30) {
// flag = 0 // flag = 1
// for (l in k downTo 0) { //
// nm = l - 1 // for (newl in k downTo 0) {
// if (abs(rv1[l]) < eps) { // nm = newl - 1
// if (abs(rv1[newl]) + anorm == anorm) {
// flag = 0 // flag = 0
//// println("break1") // l = newl
//// println("newl before break1 = " + newl)
// println("break1")
// break // break
// } // }
// if (abs(w[nm, 0]) < eps) { // if (abs(w[nm, 0] + anorm) == anorm) {
// l = newl
// println("break2") // println("break2")
// break // break
// } // }
// } // }
// //
// // l = 1 тут //// println("NEWL = " + l)
//
//// l = 0
// //
// if (flag != 0) { // if (flag != 0) {
// c = 0.0 // c = 0.0
// s = 0.0 // s = 1.0
// for (i in l until k) { // а точно ли такие границы? там немного отличается // for (i in l until k) {
// f=s*rv1[i] // f = s * rv1[i]
// rv1[i]=c*rv1[i] // rv1[i] = c * rv1[i]
// if (abs(f) < eps) { // if (abs(f) + anorm == anorm) {
// println("break3") // println("break3")
// break // break
// } // }
// g=w[i, 0] // h = pythag(f, g)
// h=pythag(f,g) // w[i, 0] = h
// w[i, 0]=h // h = 1.0 / h
// h=1.0/h // c = g * h
// c=g*h // s = (-f) * h
// s = -f*h // for (j in 0 until m) {
// for (j in 0 until m) { // точно ли такие границы? // y = this[j, nm]
// y=this[j, nm] // z = this[j, i]
// z=this[j, i] // this[j, nm] = y * c + z * s
// this[j, nm]=y*c+z*s // this[j, i] = z * c - y * 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
// } // }
// } // }
} }