### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-num-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 2 of 3
1027, -0.7791902451463131, -1.168773675878457, -1.8023918740364615,],
    [-0.7255813094089609, -0.6599481245008292, -0.6132935912751079, -0.6001656202344267, -0.9656995184226206,],
    [-0.2706435681806497, -0.37683693893675707, 0.1454012389605544, -0.2524856382131687, -0.1832861493887159,],
    [-0.25346494022471294, -0.034976909820802315, -0.40523345579661035, -0.3042534812505212, -0.7456694840765314,],
    [0.3232602101781864, 0.2187167384986574, 0.39770492818806313, 0.18484732428275386, 0.6840232171558648,],
    [-0.12627206077967576, 0.22810152457644187, -0.08358600426058016, 0.3855951837101906, 0.18596814875099849,],
    [-0.5648471145968404, -0.12230132226051842, -0.6255138866668446, -0.6991112162503683, -0.22180105879372478,],
    [-0.2514851779890921, 0.26868394934112794, 0.11129620314236724, 0.7323300527937038, -0.2037361807555147,],
    [0.5294793846838799, 0.12853319270334268, 0.255770526547394, 0.34313072822437585, -0.25510668175740325,],
    [-0.30702000923201295, 0.4330172700274453, 0.4160212523423857, -0.1256545168104251, 0.025424204736825898,];

  my $QTtau = QR-decomp($A);
  QR-QTmat($A, $QTtau, $B);
  ok ([&&] (gather for ^10 X ^5 -> ($i, $j) { take $B[$i;$j] } Z≅ (gather for ^10 X ^5 -> ($i, $j) { take @sol4[$i;$j] }))), 'QR-QTmat';

  my Math::Libgsl::Matrix $qr2  .= new(2, 2);
  my Math::Libgsl::Vector $rhs2 .= new(2);
  my Math::Libgsl::Vector $sol2 .= new(2);
  create-hilbert-matrix($qr2);
  $rhs2[$_] = $_ + 1e0 for ^2;
  $sol2[0] = -8e0;
  $sol2[1] = 18e0;

  my $d2 = QR-decomp($qr2);
  my ($q2, $r2) = QR-unpack($qr2, $d2);
  my $x2 = QR-QRsolve($q2, $r2, $rhs2);
  ok ([&&] ($x2[^2] Z≅ $sol2[^2])), 'QR-QRsolve';

  my Math::Libgsl::Matrix $qr3  .= new(3, 5);
  my Math::Libgsl::Vector $u3   .= new(3);
  my Math::Libgsl::Vector $v3   .= new(5);
  my Math::Libgsl::Vector $w3   .= new(3);
  my @sol5[3; 3] = [-0.8935042044928663, -0.12898229875765133, 0.43013230890129905],
                   [-0.2905430424289795, 0.8964013473979384, -0.3347377553836333],
                   [-0.3423959360881532, -0.42406154151296593, -0.8384132823136599];
  my @sol6[3; 5] =
    [-0.8601782732745459, 1.230031865411747, -0.4425078373682453, -0.06936085656765273, -1.4575361112775387],
    [0, -1.001388117214783, 0.09845588959891419, -0.08720390115362166, 0.7809101406822793],
    [0, 0, -0.10819274227603684, -0.0837858933928506, -0.15697099998011835];
  create-general-matrix($qr3);
  $u3[$_] = sin($_ + 1e0) for ^3;
  $v3[$_] = cos($_ + 2e0) + sin($_ * $_ + 3e0) for ^5;

  my $d3 = QR-decomp($qr3);
  my ($q3, $r3) = QR-unpack($qr3, $d3);

  for ^3 -> $j {
    my $sum = 0e0;
    for ^3 -> $i {
      $sum += $q3[$i; $j] * $u3[$i];
    }
    $w3[$j] = $sum;
  }

  QR-update($q3, $r3, $w3, $v3);

  ok ([&&] (gather for ^3 X ^3 -> ($i, $j) { take $q3[$i;$j] } Z≅ (gather for ^3 X ^3 -> ($i, $j) { take @sol5[$i;$j] })))
    &&
     ([&&] (gather for ^3 X ^5 -> ($i, $j) { take $r3[$i;$j] } Z≅ (gather for ^3 X ^5 -> ($i, $j) { take @sol6[$i;$j] }))),
    'QR-update';
}

subtest 'QR decomposition with column pivoting' => {
  my Math::Libgsl::Matrix $QR .= new: 5, 3;
  my Math::Libgsl::Vector $tau;
  my Math::Libgsl::Permutation $p;
  my Int $signum;
  my $*TOLERANCE = 1e-12;
  create-general-matrix($QR);
  my @sol1 = [-1.209797962930634, -0.4920144631391021, -0.688820248394743],
             [0.22626502892459002, -0.1404237167152902, -0.12984510634014051],
             [0.15084335261639334, 0.3706947936265647, 0.007470117121635775],
             [0.11313251446229501, 0.3598502556967054, 0.4980706174506804],
             [0.09050601156983601, 0.33463972482981114, 0.673319606696587];

  ($tau, $p, $signum) = QRPT-decomp($QR);

  ok ([&&] (gather for ^5 X ^3 -> ($i, $j) { take $QR[$i; $j] } Z≅ @sol1».[].flat))
    &&
     ([&&] ($tau[^3] Z≅ 1.8265842980736915, 1.450441410025935, 1.1754792908326592))
    &&
     ([&&] ($p.all Z== 0, 2, 1))
    &&
     $signum == -1,
    'QRPT-decomp';

  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  create-hilbert-matrix($hilb2);
  my Math::Libgsl::Matrix $QR1 .= new: 2, 2;
  my Math::Libgsl::Matrix $Q;
  my Math::Libgsl::Matrix $R;
  $QR1.copy($hilb2);
  my @solQ = [-0.8944271909999157, -0.44721359549995793],
             [-0.44721359549995793, 0.8944271909999159];
  my @solR = [-1.118033988749895, -0.5962847939999438],
             [0, 0.07453559924999298];
  ($Q, $R, $p, $signum) = QRPT-decomp2($QR1);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Q[$i; $j] } Z≅ @solQ».[].flat))
    &&
     ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $R[$i; $j] } Z≅ @solR».[].flat))
    &&
     ($p.all Z== 0, 1)
    &&
     $signum == 1,
    'QRPT-decomp2';

  my Math::Libgsl::Matrix $QR2 .= new: 2, 2;
  my Math::Libgsl::Vector $rhs .= new: 2;
  my @sol2 = -8, 18;
  $QR2.copy($hilb2);
  $rhs[$_] = $_ + 1e0 for ^2;

  ($tau, $p, $signum) = QRPT-decomp($QR2);
  my Math::Libgsl::Vector $x = QRPT-solve($QR2, $tau, $p, $rhs);
  ok ([&&] ($x[^2] Z≅ @sol2)), 'QRPT-solve';

  my Math::Libgsl::Matrix $QR3 .= new: 2, 2;
  my Math::Libgsl::Vector $u   .= new: 2;
  my Math::Libgsl::Vector $v   .= new: 2;
  my Math::Libgsl::Vector $w   .= new: 2;
  $QR3.copy($hilb2);
  $u[$_] = sin($_ + 1e0) for ^2;
  $v[$_] = cos($_ + 2e0) + sin($_ * $_ + 3e0) for ^2;
  my @solQ3 = [-0.9509858213330031, -0.3092344864719868],
              [-0.30923448647198676, 0.9509858213330031];
  my @solR3 = [-0.8081854498176496, 1.310435911733038],
              [0, -0.8935894339524059];

  ($tau, $p, $signum) = QRPT-decomp($QR3);
  ($Q, $R) = QR-unpack($QR3, $tau);

  for ^2 -> $j {
    my $sum = 0;
    for ^2 -> $i {
      $sum += $Q[$i; $j] * $u[$i];
    }
    $w[$j] = $sum;
  }

  QRPT-update($Q, $R, $p, $w, $v);

  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Q[$i; $j] } Z≅ @solQ3».[].flat))
    &&
     ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $R[$i; $j] } Z≅ @solR3».[].flat)),
     'QRPT-update';

  my Math::Libgsl::Matrix $QR4 .= new: 2, 2;
  my Math::Libgsl::Vector $res;
  $QR4.copy($hilb2);

  ($tau, $p, $signum) = QRPT-decomp($QR4);
  ($x, $res) = QRPT-lssolve($QR4, $tau, $p, $rhs);

  ok ([&&] ($x[^2] Z≅ @sol2».[].flat)), 'QRPT-lssolve';
  ok ([&&] ($res[^2] Z≅ (0, 0))), 'QRPT-lssolve residuals';

  my Math::Libgsl::Matrix $QR5 .= new: 2, 2;
  $QR5.copy($hilb2);

  ($tau, $p, $signum) = QRPT-decomp($QR5);
  my $rank = QRPT-rank($QR5, -1);
  ($x, $res) = QRPT-lssolve2($QR5, $tau, $p, $rhs, $rank);

  ok ([&&] ($x[^2] Z≅ @sol2 && $res[^2] Z≅ (0, 0))), 'QRPT-lssolve2 and residuals';
}

subtest 'LQ decomposition' => {
  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  my Math::Libgsl::Matrix $lq    .= new: 2, 2;
  my Math::Libgsl::Vector $rhs   .= new: 2;
  my $*TOLERANCE = 1e-12;
  create-hilbert-matrix($hilb2);
  $lq.copy-transpose($hilb2);
  $rhs[$_] = $_ + 1e0 for ^2;
  my @sol1 = -8, 18;

  my Math::Libgsl::Vector $d = LQ-decomp($lq);
  my Math::Libgsl::Vector $x = LQ-solve-T($lq, $d, $rhs);
  ok ([&&] ($x[^2] Z≅ @sol1)), 'LQ-decomp and LQ-solve-T';

  my ($q, $l) = LQ-unpack($lq, $d);
  $x = LQ-LQsolve($q, $l, $rhs);
  ok ([&&] ($x[^2] Z≅ @sol1)), 'LQ-unpack and LQ-LQsolve';

  my Math::Libgsl::Vector $residual;
  ($x, $residual) = LQ-lssolve-T($lq, $d, $rhs);
  ok ([&&] ($x[^2] Z≅ @sol1 && $residual[^2] Z≅ (0, 0))), 'LQ-lssolve-T and residuals';
}

subtest 'Complete Orthogonal Decomposition' => {
  my Math::Libgsl::Matrix $hilb2 .= new: 2, 2;
  my Math::Libgsl::Matrix $QRTZ  .= new: 2, 2;
  my Math::Libgsl::Matrix $Q;
  my Math::Libgsl::Matrix $R;
  my Math::Libgsl::Matrix $Z;
  my Math::Libgsl::Vector $tau-Q;
  my Math::Libgsl::Vector $tau-Z;
  my Math::Libgsl::Permutation $p;
  my Int $rank;
  my @solQRTZ1 = -1.118033988749895, -0.5962847939999438, 0.2360679774997897, 0.07453559924999298;
  my @soltauQ1 = 1.8944271909999157, 0;
  my @soltauZ1 = 0, 0;
  my @solQ1 = -0.8944271909999157, -0.44721359549995793, -0.44721359549995793, 0.8944271909999159;
  my @solR1 =-1.118033988749895, -0.5962847939999438, 0, 0.07453559924999298;
  my @solZ1 = 1, 0, 0, 1;
  create-hilbert-matrix($hilb2);
  my $*TOLERANCE = 1e-12;
  $QRTZ.copy($hilb2);

  ($tau-Q, $tau-Z, $p, $rank) = COD-decomp($QRTZ);

  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $QRTZ[$i; $j] } Z≅ @solQRTZ1))
    &&
     ([&&] ($tau-Q[^2] Z≅ @soltauQ1))
    &&
     ([&&] ($tau-Z[^2] Z≅ @soltauZ1))
    &&
      $rank == 2
    &&
      $p.all cmp (0, 1) eq 'Same',
    'COD-decomp';

  ($Q, $R, $Z) = COD-unpack($QRTZ, $tau-Q, $tau-Z, $rank);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Q[$i; $j] } Z≅ @solQ1))
    &&
     ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $R[$i; $j] } Z≅ @solR1))
    &&
     ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $Z[$i; $j] } Z≅ @solZ1)),
    'COD-unpack';

  my Math::Libgsl::Vector $x;
  my Math::Libgsl::Vector $residual;
  my Math::Libgsl::Vector $rhs .= new: 2;
  $rhs[$_] = $_ + 1e0 for ^2;

  ($x, $residual) = COD-lssolve($QRTZ, $tau-Q, $tau-Z, $p, $rank, $rhs);
  ok ([&&] ($x[^2] Z≅ (-8, 18))) && ([&&] ($residual[^2] Z≅ (0, 0))), 'COD-lssolve';

  my $lambda = 2.3;
  my Math::Libgsl::Matrix $A .= new: 100, 50;
  my Math::Libgsl::Random $r .= new;
  my Math::Libgsl::Vector $b .= new: 100;
  my @solx = -0.018462484585622357, 0.019575217751009073, 0.012640715473137413, -0.009169327686171697,
             -0.04536822106424766, 0.009123826669473341, -0.03799928794535317, 0.015231393529434173,
             -0.01562889258004679, -0.024739261667587557, -0.06108008108730741, -0.0029725782995958976,
             0.060291626542463984, -0.015578970867439236, 0.004320649763656334, -0.020416253782541305,
             -0.0031278793126650276, -0.06830124428136776, 0.0163151584806687, 0.05763413486690726,
             -0.03384175023225632, 0.0018458339310162344, 0.025280372595978887, 0.005602319166374679,
             0.10043963459510222, -0.03440777682701437, 0.06774975682369475, -0.005309708955107483,
             0.03477788960324378, -0.07427459094918588, -0.06477401548762148, -0.01715995695683791,
             -0.004913064749138054, 0.03784210527884247, 0.08678188070856876, -0.02190503466193547,
             0.053437937675444615, 0.030197104313360035, -0.0286552083076863, -0.05283492707176814,
             0.004981517482189463, 4.3112944918568944e-05, -0.009873396836500778, 0.04731192772452806,
             -0.007836638697682287, 0.0647771161650437, 0.023505975834774293, -0.027951346769336152,
             -0.020191048806093632, 0.0022644836663259713;
  my @solresidual = 0.07506552159992758, 0.3922469402268133, 0.13086496681398624, 0.171898235069902,
                    -0.465798865421583, -0.17738667573460443, 0.15558335374947, 0.07061611260292482,
                    -0.5060183552388352, -0.01709627035339239, 0.4127097735142271, 0.01468540698559401,
                    0.44060683205134393, 0.09758914841057965, -0.017985089969275073, 0.008399085064383666,
                    0.32241851030994995, 0.1256247366811053, 0.20227398986318143, 0.06165505285480713,
                    0.21831805784761152, -0.3154175853219755, -0.16321218499827234, -0.23245018586385802,
                    -0.23051098965385736, 0.34170338835721303, -0.03895174689122388, 0.43838425243233653,
                    -0.3965413684295542, 0.4123100221992462, 0.19865765410632205, 0.18345326705786552,
                    0.015448021520477171, 0.36257625430031765, -0.21713059669180795, 0.15031232718835835,
                    0.03781535235402038, -0.09934274141276186, -0.2622349867138063, 0.2537203014028584,
                    0.04612463056286892, -0.2799831927260623, 0.2523326073581078, 0.040078595169375245,
                    0.12537098382301437, 0.17660092906662037, 0.3674853031930402, 0.44687716271403377,
                    0.4641826891415557, -0.18294237412464584, -0.29701217728972185, 0.004545795944321354,
                    0.3107658018198267, -0.2504035820875857, -0.1892171518699156, -0.03903156019712084,
                    0.20474014761326123, -0.2490466922273761, -0.3819957010931251, 0.009866397059930743,
                    -0.2088024854211764, 0.29112099787475865, -0.46402669659424384, 0.37765553697981935,
                    0.41390859356011206, 0.036794514723160525, -0.18170789092485448, 0.23258826170761607,
                    -0.4118142497169783, 0.560014727781023, -0.06130241434366662, 0.03190945546976214,
                    0.41490019303245296, -0.174761217843828, 0.04458974469682233, -0.2918453899491667,
                    0.27334617031757047, -0.4076242326650097, 0.2986536971664846, 0.12720468080804648,
                    0.