kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Golub-Kahan SVD algorithm for KMP tensors #499

Closed
grinisrit wants to merge 64 commits from dev into dev
pull from: dev
merge into: kscience:dev
Conversation 1 Commits 64 Files Changed 9 +17 -13
1 changed files with 17 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 67eb92357b - Show all commits

30
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
View File

@ -385,6 +385,10 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
var g = 0.0 var g = 0.0
var l = 0 var l = 0
val epsilon = 1e-10 val epsilon = 1e-10
val wStart = w.bufferStart
val wBuffer = w.mutableBuffer
for (i in 0 until n) { for (i in 0 until n) {
/* left-hand reduction */ /* left-hand reduction */
l = i + 1 l = i + 1
@ -427,7 +431,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
} }
} }
w.mutableBuffer.array()[w.bufferStart + i] = scale * g wBuffer[wStart + i] = scale * g
/* right-hand reduction */ /* right-hand reduction */
g = 0.0 g = 0.0
s = 0.0 s = 0.0
@ -468,7 +472,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
} }
} }
} }
anorm = max(anorm, (abs(w.mutableBuffer.array()[w.bufferStart + i]) + abs(rv1[i]))); anorm = max(anorm, (abs(wBuffer[wStart + i]) + abs(rv1[i])));
} }
for (i in n - 1 downTo 0) { for (i in n - 1 downTo 0) {
@ -497,7 +501,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
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.mutableBuffer.array()[w.bufferStart + i] g = wBuffer[wStart + i]
for (j in l until n) { for (j in l until n) {
this[i, j] = 0.0 this[i, j] = 0.0
} }
@ -541,7 +545,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
l = newl l = newl
break break
} }
if (abs(w.mutableBuffer.array()[w.bufferStart + nm]) + anorm == anorm) { if (abs(wBuffer[wStart + nm]) + anorm == anorm) {
l = newl l = newl
break break
} }
@ -556,9 +560,9 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
if (abs(f) + anorm == anorm) { if (abs(f) + anorm == anorm) {
break break
} }
g = w.mutableBuffer.array()[w.bufferStart + i] g = wBuffer[wStart + i]
h = pythag(f, g) h = pythag(f, g)
w.mutableBuffer.array()[w.bufferStart + i] = h wBuffer[wStart + i] = h
h = 1.0 / h h = 1.0 / h
c = g * h c = g * h
s = (-f) * h s = (-f) * h
@ -571,19 +575,19 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
} }
} }
z = w.mutableBuffer.array()[w.bufferStart + k] z = wBuffer[wStart + k]
if (l == k) { if (l == k) {
if (z < 0.0) { if (z < 0.0) {
w.mutableBuffer.array()[w.bufferStart + k] = -z wBuffer[wStart + k] = -z
for (j in 0 until n) for (j in 0 until n)
v[j, k] = -v[j, k] v[j, k] = -v[j, k]
} }
break break
} }
x = w.mutableBuffer.array()[w.bufferStart + l] x = wBuffer[wStart + l]
nm = k - 1 nm = k - 1
y = w.mutableBuffer.array()[w.bufferStart + nm] y = wBuffer[wStart + nm]
g = rv1[nm] g = rv1[nm]
h = rv1[k] h = rv1[k]
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
@ -596,7 +600,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
for (j in l until nm + 1) { for (j in l until nm + 1) {
i = j + 1 i = j + 1
g = rv1[i] g = rv1[i]
y = w.mutableBuffer.array()[w.bufferStart + i] y = wBuffer[wStart + i]
h = s * g h = s * g
g = c * g g = c * g
z = pythag(f, h) z = pythag(f, h)
@ -615,7 +619,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
v[jj, i] = z * c - x * s; v[jj, i] = z * c - x * s;
} }
z = pythag(f, h) z = pythag(f, h)
w.mutableBuffer.array()[w.bufferStart + j] = z wBuffer[wStart + j] = z
if (abs(z) > epsilon) { if (abs(z) > epsilon) {
z = 1.0 / z z = 1.0 / z
c = f * z c = f * z
@ -632,7 +636,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2
} }
rv1[l] = 0.0 rv1[l] = 0.0
rv1[k] = f rv1[k] = f
w.mutableBuffer.array()[w.bufferStart + k] = x wBuffer[wStart + k] = x
} }
} }