### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-num-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 3 of 3
3144962251158973, -0.2838600134919259, 0.015543879118609596, 0.030197381638951814,
                    -0.0999168041215961, -0.5444261330153957, 0.07097265642802611, -0.0210716214550671,
                    0.04054050575915193, -0.22010586048390524, -0.12700298259776704, -0.1342561981471622,
                    -0.23204045784812952, -0.32304635818876315, -0.3124666552755392, -0.13770312596167988,
                    -0.504456892309686, -0.0787879072472506, -0.04135346784253771, -0.22888456254352116;
  create-rank-matrix(26, $A, $r);
  create-random-vector($b, $r);

  ($tau-Q, $tau-Z, $p, $rank) = COD-decomp($A);
  ($x, $residual) = COD-lssolve2($A, $tau-Q, $tau-Z, $p, $rank, $b, $lambda);
  ok ([&&] ($x[^50] Z≅ @solx)) && ([&&] ($residual[^100] Z≅ @solresidual)), 'COD-lssolve2';
}

subtest 'Singular Value Decomposition' => {
  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  my Math::Libgsl::Matrix $A     .= new: 2, 2;
  my Math::Libgsl::Matrix $V;
  my Math::Libgsl::Vector $S;
  my Math::Libgsl::Vector $x;
  my @solV1 = -0.8816745987679439, -0.4718579255320244, -0.4718579255320244, 0.8816745987679439;
  my @solJ1 = 0.8816745987679439, -0.4718579255320244, 0.4718579255320244, 0.8816745987679439;
  my @solS1 = 1.2675918792439984, 0.06574145408933513;
  my $*TOLERANCE = 1e-12;
  create-hilbert-matrix($hilb2);
  $A.copy($hilb2);

  ($V, $S) = SV-decomp($A);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @solV1)) && ([&&] ($S[^2] Z≅ @solS1)), 'SV-decomp';

  $A.copy($hilb2);
  ($V, $S) = SV-decomp-jacobi($A);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @solJ1)) && ([&&] ($S[^2] Z≅ @solS1)), 'SV-decomp-jacobi';

  $A.copy($hilb2);
  ($V, $S) = SV-decomp-mod($A);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @solV1)) && ([&&] ($S[^2] Z≅ @solS1)), 'SV-decomp-mod';

  $A.copy($hilb2);
  my Math::Libgsl::Vector $b .= new: 2;
  $b[$_] = $_ + 1 for ^2;
  ($V, $S) = SV-decomp($A);
  $x = SV-solve($A, $V, $S, $b);
  ok ([&&] ($x[^2] Z≅ (-8, 18))), 'SV-solve';
}

subtest 'Cholesky Decomposition' => {
  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  my Math::Libgsl::Matrix $V     .= new: 2, 2;
  my Math::Libgsl::Vector $rhs   .= new: 2;
  my @sol1 = 1, 0.5, 0.5, 0.28867513459481287;
  my @sol2 = 4, -6, -6, 12;
  my $*TOLERANCE = 1e-12;
  create-hilbert-matrix($hilb2);
  $V.copy($hilb2);

  cholesky-decomp1($V);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @sol1)), 'cholesky-decomp1';
  is-approx cholesky-rcond($V), 0.03703703703703703, 'cholesky-rcond';
  cholesky-invert($V);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @sol2)), 'cholesky-invert';

  $V.copy($hilb2);
  $rhs[$_] = $_ + 1 for ^2;
  cholesky-decomp1($V);
  my $x = cholesky-solve($V, $rhs);
  ok ([&&] ($x[^2] Z≅ (-8, 18))), 'cholesky-solve';

  $V.copy($hilb2);
  my $s = cholesky-decomp2($V);
  ok ([&&] ($s[^2] Z≅ (1, 1.7320508075688774))), 'cholesky-decomp2';
  $x = cholesky-solve2($V, $s, $rhs);
  ok ([&&] ($x[^2] Z≅ (-8, 18))), 'cholesky-solve2';

  $V.copy($hilb2);
  my $d = cholesky-decomp-unit($V);
  ok ([&&] ($d[^2] Z≅ (1, 0.08333333333333333))), 'cholesky-decomp-unit';

  my Math::Libgsl::Matrix $m .= new: 2, 2;
  my @sol3 = 0.999741748906672, 0.16290987539105117, 0.16295195791234196, 0.25607132213460243;
  create-symm-matrix($m);
  my ($p, $E) = mcholesky-decomp($m);
  ok (([&&] (gather for ^2 X ^2 -> ($i, $j) { take $m[$i; $j] } Z≅ @sol3))
      &&
     ([&&] ($E[^2] Z≅ (0, 0)))
      &&
     ($p.all cmp (0, 1) eq 'Same')), 'mcholesky-decomp';

  create-symm-matrix($m);
  ($p, $) = mcholesky-decomp($m, :!perturbation);
  ok (([&&] (gather for ^2 X ^2 -> ($i, $j) { take $m[$i; $j] } Z≅ @sol3))
      &&
     ($p.all cmp (0,1) eq 'Same')), 'mcholesky-decomp, no perturbation';

  $rhs[0] = 1.0811966866021976; $rhs[1] = 0.5541741662503531;
  my @sol4 = 0.9997417489066718, 0.16290987539105165;
  $V.copy($hilb2);
  ($p, $) = mcholesky-decomp($V, :!perturbation);
  $x = mcholesky-solve($V, $p, $rhs);
  ok ([&&] ($x[^2] Z≅ @sol4)), 'mcholesky-solve';

  my @sol5 = 4, -6.000000000000001, -6.000000000000001, 12.000000000000002;
  my $Ainv = mcholesky-invert($V, $p);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Ainv[$i; $j] } Z≅ @sol5)), 'mcholesky-invert';

  my @sol6 = 1, 0.5, 0.5, 0.08333333333333331;
  $V.copy($hilb2);
  $p = pcholesky-decomp($V);
  ok (([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ @sol6)) && $p.all cmp (0,1) eq 'Same'), 'pcholesky-decomp';
  ok pcholesky-rcond($V, $p) ≅ 0.03703703703703703, 'pcholesky-rcond';

  $V.copy($hilb2);
  $rhs[0] = 1.0811966866021976; $rhs[1] = 0.5541741662503531;
  $x = pcholesky-solve($V, $p, $rhs);
  ok ([&&] ($x[^2] Z≅ (1.0608329521783162, 0.040727468847762904))), 'pcholesky-solve';

  my @sol7 = 1.75, -1.5, -1.5, 3;
  $Ainv = pcholesky-invert($V, $p);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Ainv[$i; $j] } Z≅ @sol7)), 'pcholesky-invert';
}

subtest 'LDLT Decomposition' => {
  if $gsl-version < v2.6 {
    my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
    create-hilbert-matrix($hilb2);
    throws-like { ldlt-decomp($hilb2) },
                X::Libgsl,
                message => /'version < v2.6'/,
                'LDLT functions will be available with v2.6+';
  } else {
    skip 'version ≥ 2.6';
  }
}

subtest 'Bidiagonalization' => {
  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  my Math::Libgsl::Matrix $A     .= new: 2, 2;
  my Math::Libgsl::Vector ($tau_U, $tau_V);
  my @sol1 = -1.118033988749895, -0.5962847939999438, 0.2360679774997897, 0.07453559924999298;
  create-hilbert-matrix($hilb2);
  $A.copy($hilb2);

  ($tau_U, $tau_V) = bidiag-decomp($A);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $A[$i; $j] } Z≅ @sol1))
      &&
     ([&&] ($tau_U[^2] Z≅ (1.8944271909999157, 0)))
      &&
     $tau_V[0] == 0, 'bidiag-decomp';

  my Math::Libgsl::Vector ($diag, $sdiag);
  my Math::Libgsl::Matrix ($U, $V);
  my @sol2 = -0.8944271909999157, -0.44721359549995793, -0.44721359549995793, 0.8944271909999159;

  ($U, $V, $diag, $sdiag) = bidiag-unpack($A, $tau_U, $tau_V);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $U[$i; $j] } Z≅ @sol2))
      &&
     ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $V[$i; $j] } Z≅ (1, 0, 0, 1)))
      &&
     ([&&] ($diag[^2] Z≅ (-1.118033988749895, 0.07453559924999298)))
      &&
     $sdiag[0] == -0.5962847939999438, 'bidiag-decomp';
}

subtest 'Householder Transformations' => {
  my Math::Libgsl::Matrix $A .= new: 2, 2;
  my Math::Libgsl::Vector $x .= new: 2;
  my @sol = -8, 18;
  create-hilbert-matrix($A);
  $x[$_] = $_ + 1 for ^2;

  HH-svx($A, $x);
  ok ([&&] ($x[^2] Z≅ @sol)), 'HH-svx';
}

subtest 'Tridiagonal Systems' => {
  my Math::Libgsl::Vector $diag1    .= new: 2;
  my Math::Libgsl::Vector $offdiag1 .= new: 1;
  my Math::Libgsl::Vector $rhs1     .= new: 2;
  my Math::Libgsl::Vector $x;
  my @sol1 = 0, 2;
  for ^2 -> $i {
    $diag1[$i] = 1;
    $rhs1[$i]  = 1 + $i;
  }
  $offdiag1[0] = .5;

  $x = tridiag-symm-solve($diag1, $offdiag1, $rhs1);
  ok ([&&] ($x[^2] Z≅ @sol1)), 'tridiag-symm-solve';

  my Math::Libgsl::Vector $diag2    .= new: 3;
  my Math::Libgsl::Vector $offdiag2 .= new: 3;
  my Math::Libgsl::Vector $rhs2     .= new: 3;
  my @sol2 = -2, 5, -2;
  $diag2[$_] = 1 for ^3;
  $offdiag2[$_] = 3 for ^3;
  $rhs2[0] = 7; $rhs2[1] = -7; $rhs2[2] = 7;

  $x = tridiag-symm-cyc-solve($diag2, $offdiag2, $rhs2);
  ok ([&&] ($x[^3] Z≅ @sol2)), 'tridiag-symm-cyc-solve';

  my Math::Libgsl::Vector $diag3      .= new: 3;
  my Math::Libgsl::Vector $abovediag3 .= new: 2;
  my Math::Libgsl::Vector $belowdiag3 .= new: 2;
  my Math::Libgsl::Vector $rhs3       .= new: 3;
  my @sol3 = .75, .75, 2.625;
  for ^3 -> $i {
    $diag3[$i] = 1;
    $rhs3[$i]  = 1 + $i;
  }
  for ^2 -> $i {
    $abovediag3[$i] = 1/3;
    $belowdiag3[$i] = 1/2;
  }

  $x = tridiag-solve($diag3, $abovediag3, $belowdiag3, $rhs3);
  ok ([&&] ($x[^3] Z≅ @sol3)), 'tridiag-solve';

  my Math::Libgsl::Vector $diag4      .= new: 3;
  my Math::Libgsl::Vector $abovediag4 .= new: 3;
  my Math::Libgsl::Vector $belowdiag4 .= new: 3;
  my Math::Libgsl::Vector $rhs4       .= new: 3;
  my @sol4 = .75, .75, 2.625;
  for ^3 -> $i {
    $diag4[$i] = 1;
    $rhs4[$i]  = 1 + $i;
    $abovediag4[$i] = 2;
    $belowdiag4[$i] = 1;
  }

  my @sol5 = 1.5, -.5, .5;
  $x = tridiag-cyc-solve($diag4, $abovediag4, $belowdiag4, $rhs4);
  ok ([&&] ($x[^3] Z≅ @sol5)), 'tridiag-cyc-solve';
}

subtest 'Triangular Systems' => {
  my Math::Libgsl::Matrix $A .= new: 3, 3;
  my Math::Libgsl::Random $r .= new;
  my @sol = 1.0002583178041833, 0, 0,
           -0.5431731930411399, 3.3333333333333335, 0,
           -1.694146057887794, -1.5922278099127871, 2.061967848224576;
  for ^3 -> $i {
    for 0..$i -> $j {
      my $mij = $r.get-uniform;
      $mij = max($mij, .3) if $i == $j;
      $A[$i; $j] = $mij;
    }
  }
  if $gsl-version < v2.6 {
    tri-lower-invert($A);
  }else{
    tri-invert(CblasLower, CblasNonUnit, $A);
  }
  ok ([&&] (gather for ^3 X ^3 -> ($i, $j) { take $A[$i; $j] } Z≅ @sol)), 'tri-lower-invert';
}

if $gsl-version > v2.5 {
  subtest 'Banded Cholesky decomposition' => {
    my $M = 3;
    my $*TOLERANCE = 1e-10;
    my Math::Libgsl::Matrix $m  .= new: $M, $M;
    my Math::Libgsl::Matrix $v  .= new: $M, $M;
    my Math::Libgsl::Matrix $a  .= new: $M, $M;
    my Math::Libgsl::Matrix $lt .= new: $M, $M;

    create-posdef-band-matrix($m);
    symm-to-band($m, $v);

    cholesky-band-decomp($v);
    my $l = cholesky-band-unpack($v);

    $lt.transpose-tricpy($l, CblasLower, CblasNonUnit);
    dgemm(CblasNoTrans, CblasNoTrans, 1, $l, $lt, 0, $a);

    ok ([&&] (gather for ^$M X ^$M -> ($i, $j) { take $m[$i;$j] } Z≅
        gather for ^$M X ^$M -> ($i, $j) { take $a[$i;$j] })), 'cholesky-band-decomp and cholesky-band-unpack';

    $a.copy($m);
    cholesky-decomp1($a);
    my $rcond-expected = cholesky-rcond($a);
    my $rcond = cholesky-band-rcond($v);

    ok $rcond-expected ≅ $rcond, 'cholesky-band-rcond';

    my Math::Libgsl::Random $r   .= new;
    my Math::Libgsl::Vector $rhs .= new: $M;
    my Math::Libgsl::Vector $sol .= new: $M;
    my Math::Libgsl::Matrix $c   .= new: $M, $M;

    create-random-vector($rhs, $r);
    dsymv(CblasLower, 1, $m, $sol, 0, $rhs);

    my $x = cholesky-band-solve($v, $rhs);

    ok ($sol[$_] ≅ $x[$_] for ^$M), 'cholesky-band-solve';

    my $minv = cholesky-band-invert($v);

    dsymm(CblasLeft, CblasUpper, 1, $m, $minv, 0, $c);

    my $ret = True;
    for ^$M X ^$M -> ($i, $j) {
      $ret &&= $c[$i; $j] ≅ ($i == $j ?? 1e0 !! 0e0);
    }
    ok $ret, 'cholesky_band_invert';
  }
  subtest 'Banded LDLT decomposition' => {
    my $M = 3;
    my Math::Libgsl::Matrix $m  .= new: $M, $M;
    my Math::Libgsl::Matrix $v  .= new: $M, $M;
    my Math::Libgsl::Matrix $a  .= new: $M, $M;
    my Math::Libgsl::Matrix $l  .= new: $M, $M;
    my Math::Libgsl::Matrix $lt .= new: $M, $M;
    my Math::Libgsl::Vector $d  .= new: $M;
    my Math::Libgsl::Vector::View $vd .= new;

    create-posdef-band-matrix($m);
    symm-to-band($m, $v);
    ldlt-band-decomp($v);
    ($l, $d) = ldlt-band-unpack($v);
    $lt.transpose-tricpy($l, CblasLower, CblasNonUnit);

    for ^$M -> $i {
      my Math::Libgsl::Vector::View $c .= new;
      my $col = $l.col-view($c, $i);
      my $di = $d[$i];
      $col.scale($di);
    }
    dgemm(CblasNoTrans, CblasNoTrans, 1, $l, $lt, 0, $a);

    ok ([&&] (gather for ^$M X ^$M -> ($i, $j) { take $m[$i;$j] } Z≅
        gather for ^$M X ^$M -> ($i, $j) { take $a[$i;$j] })), 'ldlt-band-decomp';

    my Math::Libgsl::Random $r   .= new;
    my Math::Libgsl::Vector $rhs .= new: $M;
    my Math::Libgsl::Vector $sol .= new: $M;

    create-random-vector($sol, $r);
    dsymv(CblasLower, 1, $m, $sol, 0, $rhs);

    my $x = ldlt-band-solve($v, $rhs);

    ok ($sol[$_] ≅ $x[$_] for ^$M), 'ldlt-band-solve';
  }
}

done-testing;


