LU decomposition: Difference between revisions
Content added Content deleted
mNo edit summary |
|||
Line 2,995: | Line 2,995: | ||
lu_backsub(lup, transpose([1, 1, -1, -1])); |
lu_backsub(lup, transpose([1, 1, -1, -1])); |
||
/* matrix([-204], [2100], [-4740], [2940]) */</lang> |
/* matrix([-204], [2100], [-4740], [2940]) */</lang> |
||
=={{header|Nim}}== |
|||
{{trans|Kotlin}} |
|||
{{libheader|strfmt}} |
|||
The matrices are represented by static arrays rather than sequences. This allows to use 1-based indexes which, in this case, is somewhat more natural than O-based indexes. Of course, as all is static, we have to rely heavily on generics. |
|||
For display, we use the third party module "strfmt" which allows to specify dynamically the format. |
|||
<lang Nim>import macros, strutils |
|||
import strfmt |
|||
type |
|||
Matrix[M, N: static int] = array[1..M, array[1..N, float]] |
|||
SquareMatrix[N: static int] = Matrix[N, N] |
|||
# Templates to allow to use more natural notation for indexing. |
|||
template `[]`(m: Matrix; i, j: int): float = m[i][j] |
|||
template `[]=`(m: Matrix; i, j: int; val: float) = m[i][j] = val |
|||
func `*`[M, N, P: static int](a: Matrix[M, N]; b: Matrix[N, P]): Matrix[M, P] = |
|||
## Matrix multiplication. |
|||
for i in 1..M: |
|||
for j in 1..P: |
|||
for k in 1..N: |
|||
result[i, j] += a[i, k] * b[k, j] |
|||
func pivotize[N: static int](m: SquareMatrix[N]): SquareMatrix[N] = |
|||
for i in 1..N: result[i, i] = 1 |
|||
for i in 1..N: |
|||
var max = m[i, i] |
|||
var row = i |
|||
for j in i..N: |
|||
if m[j, i] > max: |
|||
max = m[j, i] |
|||
row = j |
|||
if i != row: |
|||
swap result[i], result[row] |
|||
func lu[N: static int](m: SquareMatrix[N]): tuple[l, u, p: SquareMatrix[N]] = |
|||
result.p = m.pivotize() |
|||
let m2 = result.p * m |
|||
for j in 1..N: |
|||
result.l[j, j] = 1 |
|||
for i in 1..j: |
|||
var sum = 0.0 |
|||
for k in 1..<i: sum += result.u[k, j] * result.l[i, k] |
|||
result.u[i, j] = m2[i, j] - sum |
|||
for i in j..N: |
|||
var sum = 0.0 |
|||
for k in 1..<j: sum += result.u[k, j] * result.l[i, k] |
|||
result.l[i, j] = (m2[i, j] - sum) / result.u[j, j] |
|||
proc print(m: Matrix; title, f: string) = |
|||
echo '\n', title |
|||
for i in 1..m.N: |
|||
for j in 1..m.N: |
|||
stdout.write m[i, j].format(f), " " |
|||
stdout.write '\n' |
|||
when isMainModule: |
|||
const A1: SquareMatrix[3] = [[1.0, 3.0, 5.0], |
|||
[2.0, 4.0, 7.0], |
|||
[1.0, 1.0, 0.0]] |
|||
let (l1, u1, p1) = A1.lu() |
|||
echo "\nExample 2:" |
|||
A1.print("A:", "1.0f") |
|||
l1.print("L:", "8.5f") |
|||
u1.print("U:", "8.5f") |
|||
p1.print("P:", "1.0f") |
|||
const A2: SquareMatrix[4] = [[11.0, 9.0, 24.0, 2.0], |
|||
[ 1.0, 5.0, 2.0, 6.0], |
|||
[ 3.0, 17.0, 18.0, 1.0], |
|||
[ 2.0, 5.0, 7.0, 1.0]] |
|||
let (l2, u2, p2) = A2.lu() |
|||
echo "Example 1:" |
|||
A2.print("A:", "2.0f") |
|||
l2.print("L:", "8.5f") |
|||
u2.print("U:", "8.5f") |
|||
p2.print("P:", "1.0f")</lang> |
|||
{{out}} |
|||
<pre>Example 1: |
|||
A: |
|||
1 3 5 |
|||
2 4 7 |
|||
1 1 0 |
|||
L: |
|||
1.00000 0.00000 0.00000 |
|||
0.50000 1.00000 0.00000 |
|||
0.50000 -1.00000 1.00000 |
|||
U: |
|||
2.00000 4.00000 7.00000 |
|||
0.00000 1.00000 1.50000 |
|||
0.00000 0.00000 -2.00000 |
|||
P: |
|||
0 1 0 |
|||
1 0 0 |
|||
0 0 1 |
|||
Example 2: |
|||
A: |
|||
11 9 24 2 |
|||
1 5 2 6 |
|||
3 17 18 1 |
|||
2 5 7 1 |
|||
L: |
|||
1.00000 0.00000 0.00000 0.00000 |
|||
0.27273 1.00000 0.00000 0.00000 |
|||
0.09091 0.28750 1.00000 0.00000 |
|||
0.18182 0.23125 0.00360 1.00000 |
|||
U: |
|||
11.00000 9.00000 24.00000 2.00000 |
|||
0.00000 14.54545 11.45455 0.45455 |
|||
0.00000 0.00000 -3.47500 5.68750 |
|||
0.00000 0.00000 0.00000 0.51079 |
|||
P: |
|||
1 0 0 0 |
|||
0 0 1 0 |
|||
0 1 0 0 |
|||
0 0 0 1 </pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |