Cholesky decomposition: Difference between revisions
no edit summary
m (→{{header|F Sharp|F#}}: fix heading, as suggested on the Count examples/Full list/Tier 4 talk page) |
No edit summary |
||
Line 3,809:
End Function
</lang>
=={{header|Vlang}}==
{{trans|go}}
<lang vlang>import math
// Symmetric and Lower use a packed representation that stores only
// the Lower triangle.
struct Symmetric {
order int
ele []f64
}
struct Lower {
mut:
order int
ele []f64
}
// Symmetric.print prints a square matrix from the packed representation,
// printing the upper triange as a transpose of the Lower.
fn (s Symmetric) print() {
mut row, mut diag := 1, 0
for i, e in s.ele {
print("${e:10.5f} ")
if i == diag {
for j, col := diag+row, row; col < s.order; j += col {
print("${s.ele[j]:10.5f} ")
col++
}
println('')
row++
diag += row
}
}
}
// Lower.print prints a square matrix from the packed representation,
// printing the upper triangle as all zeros.
fn (l Lower) print() {
mut row, mut diag := 1, 0
for i, e in l.ele {
print("${e:10.5f} ")
if i == diag {
for _ in row..l.order {
print("${0.0:10.5f} ")
}
println('')
row++
diag += row
}
}
}
// cholesky_lower returns the cholesky decomposition of a Symmetric real
// matrix. The matrix must be positive definite but this is not checked.
fn (a Symmetric) cholesky_lower() Lower {
mut l := Lower{a.order, []f64{len: a.ele.len}}
mut row, mut col := 1, 1
mut dr := 0 // index of diagonal element at end of row
mut dc := 0 // index of diagonal element at top of column
for i, e in a.ele {
if i < dr {
d := (e - l.ele[i]) / l.ele[dc]
l.ele[i] = d
mut ci, mut cx := col, dc
for j := i + 1; j <= dr; j++ {
cx += ci
ci++
l.ele[j] += d * l.ele[cx]
}
col++
dc += col
} else {
l.ele[i] = math.sqrt(e - l.ele[i])
row++
dr += row
col = 1
dc = 0
}
}
return l
}
fn main() {
demo(Symmetric{3, [
f64(25),
15, 18,
-5, 0, 11]})
demo(Symmetric{4, [
f64(18),
22, 70,
54, 86, 174,
42, 62, 134, 106]})
}
fn demo(a Symmetric) {
println("A:")
a.print()
println("L:")
a.cholesky_lower().print()
}</lang>
{{out}}
<pre>
A:
25.00000 15.00000 -5.00000
15.00000 18.00000 0.00000
-5.00000 0.00000 11.00000
L:
5.00000 0.00000 0.00000
3.00000 3.00000 0.00000
-1.00000 1.00000 3.00000
A:
18.00000 22.00000 54.00000 42.00000
22.00000 70.00000 86.00000 62.00000
54.00000 86.00000 174.00000 134.00000
42.00000 62.00000 134.00000 106.00000
L:
4.24264 0.00000 0.00000 0.00000
5.18545 6.56591 0.00000 0.00000
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.3926
</pre>
=={{header|Wren}}==
| |||