LU decomposition: Difference between revisions
Content added Content deleted
m (→{{header|Phix}}: added syntax colouring, made p2js compatible) |
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: bit more idiomatic) |
||
| Line 3,963: | Line 3,963: | ||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
(formerly Perl 6) |
(formerly Perl 6) |
||
{{works with|Rakudo|2015-11-20}} |
|||
| ⚫ | |||
| ⚫ | |||
<lang perl6>for ( [1, 3, 5], # Test Matrices |
<lang perl6>for ( [1, 3, 5], # Test Matrices |
||
[2, 4, 7], |
[2, 4, 7], |
||
| Line 3,977: | Line 3,976: | ||
-> @test { |
-> @test { |
||
say-it 'A Matrix', @test; |
say-it 'A Matrix', @test; |
||
say-it( |
say-it( .[0], @(.[1]) ) for 'P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix' Z, lu @test; |
||
} |
} |
||
sub lu (@a) { |
sub lu (@a) { |
||
die unless @a.&is-square; |
die unless @a.&is-square; |
||
my $n = |
my $n = @a; |
||
my @P = pivotize @a; |
my @P = pivotize @a; |
||
my @Aʼ = mmult @P, @a; |
my @Aʼ = mmult @P, @a; |
||
my @L = matrix-ident $n; |
my @L = matrix-ident $n; |
||
my @U = matrix-zero $n; |
my @U = matrix-zero $n; |
||
for ^$n -> $i { |
for ^$n X ^$n -> ($i,$j) { |
||
| ⚫ | |||
for ^$n -> $j { |
|||
else { @L[$i;$j] = (@AP[$i;$j] - [+] map { @U[$_;$j] × @L[$i;$_] }, ^$j) / @U[$j;$j] } |
|||
| ⚫ | |||
} else { |
|||
@L[$i][$j] = (@Aʼ[$i][$j] - [+] map { @U[$_][$j] * @L[$i][$_] }, ^$j) / @U[$j][$j]; |
|||
} |
|||
} |
|||
} |
} |
||
@P, @Aʼ, @L, @U; |
|||
} |
} |
||
sub pivotize (@m) { |
sub pivotize (@m) { |
||
my $size = |
my $size = @m; |
||
my @id = matrix-ident $size; |
my @id = matrix-ident $size; |
||
for ^$size -> $i { |
for ^$size -> $i { |
||
my $max = @m[$i |
my $max = @m[$i;$i]; |
||
my $row = $i; |
my $row = $i; |
||
for $i ..^ $size -> $j { |
for $i ..^ $size -> $j { |
||
if @m[$j |
if @m[$j;$i] > $max { |
||
$max = @m[$j |
$max = @m[$j;$i]; |
||
$row = $j; |
$row = $j; |
||
} |
} |
||
} |
} |
||
if $row != $i |
@id[$row, $i] = @id[$i, $row] if $row != $i; |
||
@id[$row, $i] = @id[$i, $row] |
|||
} |
|||
} |
} |
||
@id |
@id |
||
} |
} |
||
sub is-square (@m) { so @m == all @m |
sub is-square (@m) { so @m == all @m } |
||
sub matrix-zero ($n, $m = $n) { map { [ flat 0 xx $n ] }, ^$m } |
sub matrix-zero ($n, $m = $n) { map { [ flat 0 xx $n ] }, ^$m } |
||
| Line 4,028: | Line 4,019: | ||
my @p; |
my @p; |
||
for ^@a X ^@b[0] -> ($r, $c) { |
for ^@a X ^@b[0] -> ($r, $c) { |
||
@p[$r |
@p[$r;$c] += @a[$r;$_] × @b[$_;$c] for ^@b; |
||
} |
} |
||
@p |
@p |
||
| Line 4,035: | Line 4,026: | ||
sub rat-int ($num) { |
sub rat-int ($num) { |
||
return $num unless $num ~~ Rat; |
return $num unless $num ~~ Rat; |
||
return $num.narrow if $num.narrow |
return $num.narrow if $num.narrow ~~ Int; |
||
$num.nude.join: '/'; |
$num.nude.join: '/'; |
||
} |
} |
||