### ----------------------------------------------------
### -- Math::Matrix
### -- Licenses: Artistic-2.0
### -- Authors: Pierre Vigier, Herbert Breunung
### -- File: lib/Math/Matrix.rakumod
### ----------------------------------------------------
### -- Chunk 1 of 2
use Math::Matrix::Util;
use Math::Matrix::Type;

unit class Math::Matrix does Math::Matrix::Util;
use AttrX::Lazy;
use Math::Matrix::ArrayOfArray;

################################################################################
# attributes
################################################################################

has     @!rows    is required;           # primary content
has Int $!row-count is required;
has Int $!column-count is required;

has Bool $!is-zero is lazy;
has Bool $!is-identity is lazy;
has Bool $!is-square is lazy;
has Bool $!is-diagonal is lazy;
has Bool $!is-anti-diagonal is lazy;
has Bool $!is-diagonal-constant is lazy;
has Bool $!is-main-diagonal-constant is lazy;
has Bool $!is-catalecticant is lazy;
has Bool $!is-symmetric is lazy;
has Bool $!is-antisymmetric is lazy;
has Bool $!is-self-adjoint is lazy;
has Bool $!is-unitary is lazy;
has Bool $!is-orthogonal is lazy;
has Bool $!is-invertible is lazy;
has Bool $!is-positive-definite is lazy;
has Bool $!is-positive-semidefinite is lazy;

has Int     $!upper-bandwith is lazy;
has Int     $!lower-bandwith is lazy;
has Int     $!rank is lazy;
has Int     $!nullity is lazy;
has Rat     $!density is lazy;
has Numeric $!trace is lazy;
has Numeric $!determinant is lazy;
has Numeric $!condition is lazy;
has Numeric $!narrowest-element-type is lazy;
has Numeric $!widest-element-type is lazy;

has Str     $!gist;
has Rat     @!eigenvalues;

################################################################################
# private accessors - interface to util role
################################################################################

method !rows      { @!rows }
method !row-count  { $!row-count }
method !column-count{ $!column-count }
method !clone-cells  { Math::Matrix::ArrayOfArray::clone(@!rows) }

################################################################################
# public methods: constructors
################################################################################

multi method new( @m ) {
    Math::Matrix::ArrayOfArray::check-data( @m );
    self.bless( rows => @m );
}
multi method new (Str $m){
    Math::Matrix.new( $m.lines.map: { .words.map: {.Bool.Str eq $_ ?? .Bool !! .Numeric} } );
}

submethod BUILD( :@rows!, :$density, :$trace, :$determinant, :$rank, :$nullity,
                 :$is-zero, :$is-identity, :$is-symmetric, :$is-invertible) {
    @!rows         = Math::Matrix::ArrayOfArray::clone(@rows);
    $!row-count    = @rows.elems;
    $!column-count = @rows[0].elems;
    $!density      = $density if $density.defined;
    $!trace        = $trace if $trace.defined;
    $!determinant  = $determinant if $determinant.defined;
    $!rank         = $rank if $rank.defined;
    $!nullity      = $nullity if $nullity.defined;
    $!is-zero      = $is-zero if $is-zero.defined;
    $!is-identity  = $is-identity if $is-identity.defined;
    $!is-symmetric = $is-symmetric if $is-symmetric.defined;
    $!is-invertible= $is-invertible if $is-invertible.defined;
}

method clone { self.bless( rows => self!clone-cells() ) }

multi method new-zero(PosInt $size) {
    self.bless( rows => Math::Matrix::ArrayOfArray::new-zero($size, $size),
            determinant => 0, rank => 0, nullity => $size, density => 0.0, trace => 0,
            is-zero => True, is-identity => False, is-diagonal => True,
            is-square => True, is-symmetric => True, is-invertible => False );
}
multi method new-zero(Math::Matrix:U: PosInt $rows, PosInt $cols) {
    self.bless( rows => Math::Matrix::ArrayOfArray::new-zero($rows, $cols),
            determinant => 0, rank => 0, nullity => min($rows, $cols), density => 0.0, trace => 0,
            is-zero => True, is-identity => False, is-diagonal => ($cols == $rows), is-invertible => False );
}

method new-identity( Int $size where * > 0 ) {
    self.bless( rows => Math::Matrix::ArrayOfArray::new-identity($size),
                determinant => 1, rank => $size, nullity => 0, density => 1/$size, trace => $size,
                is-zero => False, is-identity => True,
                is-square => True, is-diagonal => True, is-symmetric => True, is-invertible => True );
}

method new-diagonal( *@diag ){
    fail "Expect at least on number as parameter" if @diag == 0;
    fail "Expect an List of Number" unless @diag ~~ NumList;
    my Int $size = +@diag;
    my $det = [*](@diag.flat);
    self.bless( rows => Math::Matrix::ArrayOfArray::new-diagonal(@diag),
                determinant => $det, trace => [+] (@diag.flat),
                is-zero => False, is-square => True, is-diagonal => True, is-symmetric => True,
                is-invertible => $det != 0  );
}

method new-lower-triangular( @m ) {
    #don't want to trust outside of the class that a matrix is really triangular
    self.bless( rows => @m );
}

method new-upper-triangular( @m ) {
    #don't want to trust outside of the class that a matrix is really triangular
    self.bless( rows => @m );
}

method new-vector-product (@column_vector, @row_vector){
    fail "Expect two Arrays of Number" unless @column_vector ~~ NumArray and @row_vector ~~ NumArray;
    my @p;
    for ^+@column_vector X ^+@row_vector -> ($r, $c) {
        @p[$r][$c] = @column_vector[$r] * @row_vector[$c]
    }
    my $rank = ( all(@column_vector) ~~ 0 or all(@row_vector) ~~ 0) ?? 0 !! 1;
    self.bless( rows => @p, determinant => 0 , rank => $rank,
                is-invertible => (@column_vector.elems == @row_vector.elems and $rank >= @row_vector.elems) );
}

################################################################################
# end of constructor - start accessors
################################################################################

method element(Math::Matrix:D: Int:D $row, Int:D $column --> Numeric ) {
    self!check-index($row, $column);
    @!rows[$row][$column];
}

multi method AT-POS (Math::Matrix:D: Int:D $row){
    self!check-row-index($row);
    @!rows[$row];
}

method row(Math::Matrix:D: Int:D $row  --> List) {
    self!check-row-index($row);
    @!rows[$row];
}

method column(Math::Matrix:D: Int:D $column --> List) {
    self!check-column-index($column);
    (@!rows.keys.map:{ @!rows[$_;$column] }).list;
}

method diagonal(Math::Matrix:D: $start? = 0 --> List){
    fail "requested diagonal is outside of matrix boundaries" if $start >= $!row-count or $start <= -$!column-count;
    ($start > 0 ?? map { @!rows[$^i+$start;$^i] }, ^min($!column-count, $!row-count - $start)
                !! map { @!rows[$^i;$^i-$start] }, ^min($!row-count, $!column-count + $start) ).list;
}

method skew-diagonal(Math::Matrix:D: $start? = 0 --> List){
    fail "skew diagonal is only defined for square matrices" unless $.is-square;
    fail "requested skew diagonal is outside of matrix boundaries" if $start.abs >= $!row-count;
    ($start > 0 ?? map { @!rows[$!row-count -1 -$^i; $start+$^i] }, ^($!row-count - $start)
                !! map { @!rows[$!row-count -1 -$^i +$start;$^i] }, ^($!row-count + $start) ).list;
}

multi method submatrix(Math::Matrix:D: Int:D $row, Int:D $column --> Math::Matrix:D ){
    self!check-index($row, $column);
    my @rows = ^$!row-count;     @rows.splice($row,1);
    my @cols = ^$!column-count;  @cols.splice($column,1);
    self.submatrix( rows => @rows , columns => @cols);
}
multi method submatrix(Math::Matrix:D: :@rows    = (^$!row-count).list,
                                       :@columns = (^$!column-count).list --> Math::Matrix:D) {
    my @r = @rows.max    == Inf ?? (@rows.min    .. $!row-count-1).list    !! @rows.list;
    my @c = @columns.max == Inf ?? (@columns.min .. $!column-count-1).list !! @columns.list;
    fail "Need at least one row number" if @r == 0;
    fail "Need at least one column number" if @c == 0;
    self!check-indices(@r, @c);
    Math::Matrix.new([ @r.map( { [ @!rows[$^row][|@c] ] } ) ]);
}

################################################################################
# end of accessors - start with type conversion and handy shortcuts
################################################################################

method Bool(        Math::Matrix:D: --> Bool)   { ! self.is-zero }
method Numeric (    Math::Matrix:D: --> Numeric){   self.norm   }
method Str(         Math::Matrix:D: --> Str)    {   join("\n", @!rows.map: *.Str) }
method Range(       Math::Matrix:D: --> Range)  {   self.list.minmax }
method Array(       Math::Matrix:D: --> Array)  {   self!clone-cells }
method list(        Math::Matrix:D: --> List)   {   self.list-rows.flat.list }
method list-rows(   Math::Matrix:D: --> List)   {  (@!rows.map: {.flat}).list }
method list-columns(Math::Matrix:D: --> List)   {  ((^$!column-count).map: {self.column($_)}).list }
method Hash(        Math::Matrix:D: --> Hash)   {  ((^$!row-count).map: {$_ => @!rows[$_].kv.Hash}).Hash}

multi method gist(Math::Matrix:U: --> Str) { "({self.^name})" }
multi method gist(Math::Matrix:D: Int :$max-chars?, Int :$max-rows? --> Str) {
    if not $!gist.defined or ($max-chars.defined or $max-rows.defined) {
        my $max-width = (not $max-chars.defined or $max-chars < 5) ?? 80 !! $max-chars;
        my $max-heigth = (not $max-rows.defined or $max-rows < 2) ?? 20 !! $max-rows;
        my @fmt-content = @!rows.map: {    # all values in optimized complex format
            (.map: { $_ ~~ Bool   ?? %( re => $_,             im => '' ) !!
                    $_ ~~ Complex ?? %( re => $_.re.fmt("%g"),im => (($_.im >= 0 ??'+'!!'')~$_.im.fmt("%g")~'i') ) !!
                                     %( re => $_.fmt("%g"),   im => '' )
        }).Array};
        my @col-width;                     # width of the formatted element content in n column
        @fmt-content.map: {
            for .kv -> $ci, $val {
                @col-width[$ci]<re>.push: $val<re>.chars;
                @col-width[$ci]<im>.push: $val<im>.chars;
        }};
        my @max-width = @col-width.map: { %( re => $_<re>.max, im => $_<im>.max ) };
        my ($shown-cols, $width-index);
        for @max-width.kv -> $ci, $max {
            $width-index += 2 + $max<re> + $max<im>;
            if ($ci < @max-width.end and $width-index <= $max-width-3)
            or $width-index <= $max-width {$shown-cols++}
            else                          {last}
        }
        my $shown-rows = min @!rows.elems, $max-heigth;
        my $out;
        for @fmt-content.kv -> $ri, $row {
            if $ri == $shown-rows {$out ~= "  ...\n" ; last}
            for $row.kv -> $ci, $val {
                if $ci == $shown-cols { $out ~= ' ..' ; last}
                $out ~= (' ' x (@max-width[$ci]<re> - @col-width[$ci]<re>[$ri]) + 2) ~ $val<re> ~
                        $val<im> ~ (' ' x (@max-width[$ci]<im> - @col-width[$ci]<im>[$ri]));
            }
            $out ~= "\n";
        };
        $!gist = $out.chomp;
    }
    $!gist;
}

multi method raku(Math::Matrix:D: --> Str){ self.WHAT.raku ~ ".new(" ~ @!rows.raku ~ ")" }

################################################################################
# end of type conversion and handy shortcuts - start boolean matrix properties
################################################################################

method !build_is-square( Math::Matrix:D: --> Bool)   {$!column-count == $!row-count}
method !build_is-zero(    Math::Matrix:D: --> Bool)  {self.density() == 0}
method !build_is-identity( Math::Matrix:D: --> Bool) {$.is-diagonal and $.is-main-diagonal-constant and @!rows[0][0] == 1}
method !build_is-diagonal( Math::Matrix:D: --> Bool) {self.is-square and $.lower-bandwith == $.upper-bandwith == 0}
method !build_is-main-diagonal-constant(Math::Matrix:D: --> Bool){ [==](($.diagonal.flat).flat) }

method !build_is-diagonal-constant( Math::Matrix:D: --> Bool) {
    [&&](map { [==] $.diagonal($_).list }, -$!column-count+1 .. $!row-count-1);
}
method !build_is-catalecticant( Math::Matrix:D: --> Bool) {
    $.is-square and [&&](map { [==] $.skew-diagonal($_).list }, -$!column-count+1 .. $!row-count-1);
}

method is-triangular(Math::Matrix:D: Bool :$strict, Bool :$unit, Bool :$atomic, Bool :$upper,  Bool :$lower --> Bool) {
    return False unless $.is-square;
    return False if $unit and $strict;
    return False if $atomic and ($strict or ($unit.defined and not $unit));
    return False if $strict.defined and ($strict xor ($.is-main-diagonal-constant and @!rows[0][0] == 0));
    return False if $unit.defined and ($unit xor ($.is-main-diagonal-constant and @!rows[0][0] == 1));
    return False if $atomic.defined and ($atomic xor self!is-frobenius );
    return False if $upper.defined and ($upper xor $.lower-bandwith == 0);
    return False if $lower.defined and ($lower xor $.upper-bandwith == 0);
    $.lower-bandwith == 0 or $.upper-bandwith == 0;
}
method !is-frobenius(Math::Matrix:D: --> Bool){
    my $col;
    for ^$!row-count X ^$!column-count -> ($r, $c) {
        next if $r == $c;
        next if @!rows[$r][$c] == 0;
        if $col.defined { return False if $col != $c }
        else            { $col = $c  }
    }
    True;
}

method !build_is-anti-diagonal(Math::Matrix:D: --> Bool){
    for ^$!row-count X ^$!column-count -> ($r, $c) {
        next if $r == $!column-count - $c - 1 or @!rows[$r][$c] == 0;
        return False;
    }
    True;
}

method is-tridiagonal(Math::Matrix:D: --> Bool) {
    $.is-square and $.lower-bandwith < 2 and $.upper-bandwith < 2;
}

method is-diagonally-dominant( Math::Matrix:D: Bool :$strict = False,
                              Str :$along where {$^orient eq any <column row both>} = 'column' --> Bool) {
    return False unless self.is-square;
    my $greater = $strict ?? &[>] !! &[>=];
    my Bool $colwise;
    if $along ~~ any <column both> {
        $colwise = [and] map {my $c = $_; &$greater( @!rows[$c][$c] * 2,
                                                     [+](map {abs $_[$c]}, @!rows)) }, ^$!row-count;
    }
    return $colwise if $along eq 'column';
    my Bool $rowwise = [and] map { &$greater( @!rows[$^r][$^r] * 2,
                                              [+](map {abs $^c}, @!rows[$^r].flat)) }, ^$!row-count;
    return $rowwise if $along eq 'row';
    $colwise and $rowwise;
}

method !build_is-symmetric( Math::Matrix:D: --> Bool) {
    return False unless self.is-square;
    return True if $!row-count < 2;
    for ^($!row-count - 1) -> $r {
        for $r ^..^ $!row-count -> $c {
            return False unless @!rows[$r][$c] == @!rows[$c][$r];
        }
    }
    True;
}
method !build_is-antisymmetric(Math::Matrix:D: --> Bool) {
    return False unless self.is-square;
    return True if $!row-count < 2;
    for ^($!row-count - 1) -> $r {
        for $r ^..^ $!row-count -> $c {
            return False unless @!rows[$r][$c] == - @!rows[$c][$r];
        }
    }
    True;
}

method !build_is-self-adjoint(Math::Matrix:D: --> Bool) {
    self.is-square and self.T.conj ~~ self
}

method !build_is-unitary(Math::Matrix:D: --> Bool) {
    self.is-square and self.dot-product( self.T.conj ) ~~ Math::Matrix.new-identity( $!row-count );
}

method !build_is-orthogonal(Math::Matrix:D: --> Bool) {
    self.is-square and self.dot-product( self.T ) ~~ Math::Matrix.new-identity( $!row-count );
}

method !build_is-invertible(Math::Matrix:D: --> Bool) {
    self.is-square and self.determinant != 0;
}

# method is-definite(Math::Matrix:D: Bool :positive, Bool :negative, Bool :semi --> Bool){}

method !build_is-positive-definite (Math::Matrix:D: --> Bool) { # with Sylvester's criterion
    return False unless self.is-square;
    return False unless self.determinant > 0;
    my $sub = Math::Matrix.new( @!rows );
    for $!row-count - 1 ... 1 -> $r {
        $sub = $sub.submatrix(rows => 0..$r, columns => 0..$r);
        return False unless $sub.determinant > 0;
    }
    True;
}
method !build_is-positive-semidefinite (Math::Matrix:D: --> Bool) { # with Sylvester's criterion
    return False unless self.is-square;
    return False unless self.determinant >= 0;
    my $sub = Math::Matrix.new( @!rows );
    for $!row-count - 1 ... 1 -> $r {
        $sub = $sub.submatrix(rows => 0..$r, columns => 0..$r);
        return False unless $sub.determinant >= 0;
    }
    True;
}

################################################################################
# end of boolean matrix properties - start numeric matrix properties
################################################################################

method size(Math::Matrix:D: --> List)          {  $!row-count, $!column-count }

method !build_density(Math::Matrix:D: --> Rat) {
    my $valcount = 0;
    for ^$!row-count X ^$!column-count -> ($r, $c) { $valcount++ if @!rows[$r][$c] != 0 }
    $valcount / self.elems;
}

method !build_upper-bandwith(Math::Matrix:D: --> Int) {
    for $!column-count-1 ... 1  -> $i {
        return $i unless [&&](map * == 0, $.diagonal(-$i).list)
    }
    0;
}
method !build_lower-bandwith(Math::Matrix:D: --> Int) {
    for $!row-count-1 ... 1  -> $i {
        return $i unless [&&](map * == 0, $.diagonal($i).list)
    }
    0;
}
method bandwith(Math::Matrix:D: Str $which = '' --> Int) { max $.upper-bandwith, $.lower-bandwith }

method !build_trace(Math::Matrix:D: --> Numeric) {
    fail "trace is only defined for a square matrix" unless self.is-square;
    self.diagonal.sum;
}

method !build_rank(Math::Matrix:D: --> Int) {
    my $rank = 0;
    my @clone =  @!rows.clone();
    for ^$!column-count -> $c {            # make upper triangle via gauss elimination
        last if $rank == $!row-count;      # rank cant get bigger thean dim
        my $swap_row_nr = $rank;
        $swap_row_nr++ while $swap_row_nr < $!row-count and @clone[$swap_row_nr][$c] == 0;
        next if $swap_row_nr == $!row-count;
        (@clone[$rank], @clone[$swap_row_nr]) = (@clone[$swap_row_nr], @clone[$rank]);
        for $rank + 1 ..^ $!row-count -> $r {
            next if @clone[$r][$c] == 0;
            my $q = @clone[$rank][$c] / @clone[$r][$c];
            @clone[$r] = @clone[$rank] >>-<< $q <<*<< @clone[$r];
        }
        $rank++;
    }
    $rank;
}
method !build_nullity(Math::Matrix:D: --> Int) {
    min(self.size) - self.rank;
}

method det(Math::Matrix:D: --> Numeric )        { self.determinant }  # the usual short name
method !build_determinant(Math::Matrix:D: --> Numeric) {
    fail "number of columns has to be same as number of rows" unless self.is-square;
    return 1            if $!row-count == 0;
    return @!rows[0][0] if $!row-count == 1;
    if $!row-count > 4 {
        #up to 4x4 naive method is fully usable
        return [*] $.diagonal.flat if $.is-triangular();
        try {
            my ($L, $U, $P) = $.decompositionLU();
            return $P.inverted.det * $L.det * $U.det;
        }
    }
    my $det = 0;
    for ( σ_permutations([^$!row-count]) ) {
        my $permutation = .key;
        my $product = .value;
        for $permutation.kv -> $i, $j { $product *= @!rows[$i][$j] };
        $det += $product;
    }
    $!determinant = $det;
}
method determinant-naive(Math::Matrix:D: --> Numeric) {
    fail "Number of columns has to be same as number of rows" unless self.is-square;
    return 1            if $!row-count == 0;
    return @!rows[0][0] if $!row-count == 1;
    my $det = 0;
    for ( σ_permutations([^$!row-count]) ) {
        my $permutation = .key;
        my $product = .value;
        for $permutation.kv -> $i, $j { $product *= @!rows[$i][$j] };
        $det += $product;
    }
    $det;
}
sub insert ($x, @xs) { ([flat @xs[0 ..^ $_], $x, @xs[$_ .. *]] for 0 .. @xs) }
sub order ($sg, @xs) { $sg > 0 ?? @xs !! @xs.reverse }

multi sub σ_permutations ([]) { [] => 1 }
multi sub σ_permutations ([$x, *@xs]) {
    σ_permutations(@xs).map({ |order($_.value, insert($x, $_.key)) }) Z=> |(1,-1) xx *
}

method minor(Math::Matrix:D: Int:D $row, Int:D $col --> Numeric) { $.submatrix($row, $col).determinant }

multi method norm(Math::Matrix:D: PosInt :$p = 2, PosInt :$q = $p --> Numeric) {
    my $norm = 0;
    for ^$!column-count -> $c {
        my $col_sum = 0;
        for ^$!row-count -> $r {  $col_sum += abs(@!rows[$r][$c]) ** $p }
        $norm += $col_sum ** ($q / $p);
    }
    $norm ** (1/$q);
}
multi method norm(Math::Matrix:D: PosInt $p   --> Numeric){ self.norm(:p<$p>,:q<$p>)}
multi method norm(Math::Matrix:D: 'frobenius' --> Numeric){ self.norm(:p<2>, :q<2>)}
multi method norm(Math::Matrix:D: 'euclidean' --> Numeric){ self.norm(:p<2>, :q<2>)}

multi method norm(Math::Matrix:D: 'max' --> Numeric)      { max            @!rows.map: {max .map: *.abs} }
multi method norm(Math::Matrix:D: 'row-sum' --> Numeric)  { max            @!rows.map: {[+] .map: *.abs} }
multi method norm(Math::Matrix:D: 'column-sum'--> Numeric){ max (^$!column-count).map: {[+] self.column($_).map: *.abs} }

method !build_condition(Math::Matrix:D:              --> Numeric) { $.norm() * $.inverted.norm       }

method !build_narrowest-element-type(Math::Matrix:D: --> Numeric){
    return Bool if any( @!rows[*;*] ) ~~ Bool;
    return Int  if any( @!rows[*;*] ) ~~ Int;
    return Num  if any( @!rows[*;*] ) ~~ Num;
    return Rat  if any( @!rows[*;*] ) ~~ Rat;
    return FatRat if any( @!rows[*;*] ) ~~ FatRat;
    Complex;
}
method !build_widest-element-type(Math::Matrix:D: --> Numeric){
    return Complex if any( @!rows[*;*] ) ~~ Complex;
    return FatRat if any( @!rows[*;*] ) ~~ FatRat;
    return Rat   if any( @!rows[*;*] ) ~~ Rat;
    return Num  if any( @!rows[*;*] ) ~~ Num;
    return Int if any( @!rows[*;*] ) ~~ Int;
    Bool;
}

################################################################################
# end of numeric matrix properties - start create derivative matrices
################################################################################

method T(         Math::Matrix:D: --> Math::Matrix:D ) { self.transposed }
method transposed(Math::Matrix:D: --> Math::Matrix:D ) {
    my @transposed;
    for ^$!row-count X ^$!column-count -> ($r, $c) { @transposed[$c][$r] = @!rows[$r][$c] }
    Math::Matrix.new( @transposed ); # TODO use bless to preserve properties
}

method negated(Math::Matrix:D: --> Math::Matrix:D )       { self.map( - * ) }

method H(         Math::Matrix:D: --> Math::Matrix:D )    { self.conjugated.transposed }

method conj(      Math::Matrix:D: --> Math::Matrix:D )    { self.conjugated }
method conjugated(Math::Matrix:D: --> Math::Matrix:D )    { self.map( { $_.conj} ) }

method adjugated( Math::Matrix:D: --> Math::Matrix:D) {
    fail "Number of columns has to be same as number of rows" unless self.is-square;
    $!row-count == 1 ?? self.new([[1]])
                     !! self.map-index({ self.minor($^m, $^n) * self.cofactor-sign($^m, $^n) });
}

method inverted(Math::Matrix:D: --> Math::Matrix:D) {
    fail "Number of columns has to be same as number of rows" unless self.is-square;
    # fail "Matrix is not invertible, or singular because defect (determinant = 0)" if self.determinant == 0;
    my @clone = self!clone-cells();
    my @inverted = Math::Matrix::ArrayOfArray::new-identity( $!row-count );
    for ^$!row-count -> $c {
        my $swap_row_nr = $c;       # make sure that diagonal element != 0, later == 1
        $swap_row_nr++ while $swap_row_nr < $!row-count && @clone[$swap_row_nr][$c] == 0;
        fail "Matrix is not invertible, or singular because defect (determinant = 0)" if $swap_row_nr >= $!row-count;
        (@clone[$c], @clone[$swap_row_nr])       = (@clone[$swap_row_nr], @clone[$c]);
        (@inverted[$c], @inverted[$swap_row_nr]) = (@inverted[$swap_row_nr], @inverted[$c]);
        @inverted[$c] =  @inverted[$c] >>/>>  @clone[$c][$c];
        @clone[$c]    =  @clone[$c]    >>/>>  @clone[$c][$c];
        for $c + 1 ..^ $!row-count -> $r {
            @inverted[$r] = @inverted[$r]  >>-<<  @clone[$r][$c] <<*<< @inverted[$c];
            @clone[$r]    = @clone[$r]  >>-<<  @clone[$r][$c] <<*<< @clone[$c];
        }
    }
    for reverse(1 ..^ $!column-count) -> $c {
        for ^$c -> $r {
            @inverted[$r] = @inverted[$r]  >>-<<  @clone[$r][$c] <<*<< @inverted[$c];
            @clone[$r]    = @clone[$r]  >>-<<  @clone[$r][$c] <<*<< @clone[$c];
        }
    }
    Math::Matrix.new( @inverted );
}

method rref(                    Math::Matrix:D: --> Math::Matrix:D) { self.reduced-row-echelon-form }
method reduced-row-echelon-form(Math::Matrix:D: --> Math::Matrix:D) {
    my @ref = self!clone-cells();
    my $lead = 0;
    MAIN: for ^$!row-count -> $r {
        last MAIN if $lead >= $!column-count;
        my $i = $r;
        while @ref[$i][$lead] == 0 {
            $i++;
            if $!row-count == $i {
                $i = $r;
                $lead++;
                last MAIN if $lead == $!column-count;
            }
        }
        @ref[$i, $r] = @ref[$r, $i];
        my $lead_value = @ref[$r][$lead];
        @ref[$r] »/=» $lead_value;
        for ^$!row-count -> $n {
            next if $n == $r;
            @ref[$n] »-=» @ref[$r] »*» @ref[$n][$lead];
        }
        $lead++;
    }
    return Math::Matrix.new( @ref );
}

################################################################################
# end of derivative matrices - start decompositions
################################################################################

# LU factorization with optional partial pivoting and optional diagonal matrix
method LU-decomposition(Math::Matrix:D: Bool :$pivot = False, Bool :$diagonal = False, Bool :$Crout = False) {
    fail "Not an square matrix" unless self.is-square;
    fail "Has to be invertible when not using pivoting" if not $pivot and not self.is-invertible;
    fail "Crout algorithm makes only sense when not decomposing into a diagonal matrix" if $Crout and $diagonal;

    my $size = $!row-count;
    my @L = Math::Matrix::ArrayOfArray::new-identity( $size );
    my Array @U = self!clone-cells( );
    my @P = Math::Matrix::ArrayOfArray::new-identity( $size );

    for 0 .. $size-2 -> $c {
        if $pivot {
            my $maxrow = $c;
            for $c+1 ..^$size -> $r { $maxrow = $c if @U[$maxrow][$c] < @U[$r][$c] }
            (@U[$maxrow], @U[$c]) = (@U[$c], @U[$maxrow]);
            (@P[$maxrow], @P[$c]) = (@P[$c], @P[$maxrow]);
        }

        for $c+1 ..^$size -> $r {
            next if @U[$r][$c] == 0;
            my $q = @L[$r][$c] = @U[$r][$c] / @U[$c][$c];
	    @U[$r;*] = @U[$r] >>-<< $q <<*<< @U[$c];
        }
    }
    #@U = self!AoA-clone( @U );

    if $diagonal {
        my @D;
        for 0 ..^ $size -> $c {
            push @D, @U[$c][$c];
            @U[$c][$c] = 1;
        }
        $pivot ?? (Math::Matrix.new-lower-triangular(@L), Math::Matrix.new-diagonal(@D),
                   Math::Matrix.new-upper-triangular(@U), Math::Matrix.new(@P))
               !! (Math::Matrix.new-lower-triangular(@L), Math::Matrix.new-diagonal(@D),
                   Math::Matrix.new-upper-triangular(@U));
    } elsif $Crout {
        for 0 ..^ $size -> $i {
            next if @U[$i][$i] == 1;
            fail "Matrix can not be inverted" if @U[$i][$i] == 0;
            my $f = @U[$i][$i];
            for $i ..^ $size -> $j {
                @L[$j][$i] *= $f;
                @U[$i][$j] /= $f;
            }
        }
    }
    $pivot ?? (Math::Matrix.new-lower-triangular(@L), Math::Matrix.new-upper-triangular(@U), Math::Matrix.new(@P))
           !! (Math::Matrix.new-lower-triangular(@L), Math::Matrix.new-upper-triangular(@U));
}

method Cholesky-decomposition(Math::Matrix:D: Bool :$diagonal = False) {
    fail "Matrix is not symmetric"         unless self.is-symmetric;
    fail "Matrix is not positive definite" unless self.is-positive-definite;
    my @L = self!clone-cells();
    for 0 ..^$!row-count -> $k {
        @L[$k][$k] -= @L[$k][$_]**2 for 0 .. $k-1;
        @L[$k][$k]  = @L[$k][$k].sqrt;
        for $k+1 ..^ $!row-count -> $i {
            @L[$i][$k] -= @L[$