### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-num-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 1 of 3
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::LinearAlgebra;
use Math::Libgsl::Vector;
use Math::Libgsl::Matrix;
use Math::Libgsl::Permutation;
use Math::Libgsl::Random;
use Math::Libgsl::BLAS;
use Math::Libgsl::Constants;
use NativeCall;

sub create-random-matrix(Math::Libgsl::Matrix $m, Math::Libgsl::Random $r --> Int)
{
  for ^$m.matrix.size1 X ^$m.matrix.size2 -> ($i, $j) {
    $m[$i;$j] = $r.get-uniform-pos;
  }
  return GSL_SUCCESS;
}

sub create-random-vector(Math::Libgsl::Vector $v, Math::Libgsl::Random $r --> Int)
{
  for ^$v.vector.size -> $i {
    $v[$i] = $r.get-uniform-pos;
  }
  return GSL_SUCCESS;
}

sub create-general-matrix(Math::Libgsl::Matrix $m)
{
  for ^$m.matrix.size1 X ^$m.matrix.size2 -> ($i, $j) {
    $m[$i; $j] = 1e0 / ($i + $j + 1e0);
  }
}

sub create-hilbert-matrix(Math::Libgsl::Matrix $m)
{
  for ^$m.matrix.size1 X ^$m.matrix.size2 -> ($i, $j) {
    $m[$i; $j] = 1e0/($i + $j + 1e0);
  }
}

sub create-rank-matrix(size_t $rank, Math::Libgsl::Matrix $m, Math::Libgsl::Random $r --> Int)
{
  my $M = $m.matrix.size1;
  my $N = $m.matrix.size2;
  my Math::Libgsl::Vector $u .= new: $M;
  my Math::Libgsl::Vector $v .= new: $N;

  for ^$rank -> $i {
    create-random-vector($u, $r);
    create-random-vector($v, $r);
    dger(1e0, $u, $v, $m);
  }

  return GSL_SUCCESS;
}

sub create-symm-matrix(Math::Libgsl::Matrix $m)
{
  my Math::Libgsl::Random $r .= new;
  for ^$m.matrix.size1 -> $i {
    for 0..$i -> $j {
      $m[$i; $j] = $r.get-uniform-pos;
    }
  }
  $m.transpose-tricpy($m, CblasLower, CblasNonUnit);
  return GSL_SUCCESS;
}

sub create-symm-band-matrix(Math::Libgsl::Matrix $m)
{
  my Math::Libgsl::Random $r .= new;
  for ^$m.size1 -> $i {
    my Math::Libgsl::Vector::View $vv1 .= new;
    my Math::Libgsl::Vector $v1 = $m.subdiagonal-view($vv1, $i);
    create-random-vector($v1, $r);
    if $i > 0 {
      my Math::Libgsl::Vector::View $vv2 .= new;
      my Math::Libgsl::Vector $v2 = $m.superdiagonal-view($vv2, $i);
      $v2.copy($v1);
    }
  }
}

sub create-posdef-band-matrix(Math::Libgsl::Matrix $m)
{
  create-symm-band-matrix($m);
  my $size = $m.matrix.size1;
  my $α = $size * 10e0;
  $m[$_; $_] = $m[$_; $_] + $α for ^$size;
}

sub symm-to-band(Math::Libgsl::Matrix $m, Math::Libgsl::Matrix $u)
{
  my $M = $m.matrix.size1;
  my Math::Libgsl::Vector::View $view1 .= new;
  my Math::Libgsl::Vector::View $view2 .= new;
  for ^$M {
    my Math::Libgsl::Vector $vv1 = $m.subdiagonal-view($view1, $_);
    my Math::Libgsl::Vector $vv2 = $u.subcol-view($view2, $_, 0, $M - $_);
    $vv2.copy($vv1);
  }
}

subtest 'LU decomposition' => {
  my $*TOLERANCE = 1e-10;
  my Math::Libgsl::Random $r      .= new;
  my Math::Libgsl::Matrix $m      .= new(10, 10);
  my Math::Libgsl::Matrix $A      .= new(10, 10);
  my Math::Libgsl::Vector $b      .= new(10);
  ok create-random-matrix($m, $r) == GSL_SUCCESS, 'can create a random matrix';
  ok create-random-vector($b, $r) == GSL_SUCCESS, 'can create a random vector';

  $A.copy($m);
  my ($signum, $p) = LU-decomp($A);
  ok $signum == -1, 'permutation sign';
  my @sol1 = 0.999741748906672e0, 0.16290987539105117e0, 0.28261780529282987e0, 0.9472010820172727e0,
             0.23165654274635017e0, 0.48497361433692276e0, 0.9574769565369934e0, 0.7443053431343287e0,
             0.5400436583440751e0, 0.7399529814720154e0, 0.499022899013741e0, 0.6425063204525122e0,
             0.0035177683747059374e0, -0.4303718665717976e0, 0.19518510304267706e0, -0.12037178526907077e0,
             -0.2357332556594393e0, 0.009633445365608229e0, 0.17063438514690155e0, 0.23054241973420175e0,
             0.18245795991271174e0, 0.7491098476129922e0, 0.686587142174087e0, 0.515560053175852e0,
             -0.027673826178128774e0, 0.4617911786127301e0, 0.6297272827275225e0, 0.5345827448408003e0,
             0.4718373602385013e0, 0.1708249594233042e0, 0.030343040333473376e0, 0.5114285392398699e0,
             0.2676634453699523e0, 0.9970836950352449e0, 0.48048655070164487e0, 0.8215462322357204e0,
             0.5885169099991214e0, 0.32801074453308043e0, 0.33068078288375435e0, -0.0037973478102040215e0,
             0.29660212857391954e0, 0.10753126521068211e0, -0.031002105323639315e0, 0.4307016533401094e0,
             0.42791590511605626e0, 0.16670586055124165e0, 0.2212917196271707e0, -0.24692378187564057e0,
             0.3929437579470931e0, -0.1295182202387101e0, 0.7601401051976446e0, 0.8323689084541823e0,
             0.14256080005331967e0, 0.37020662716221897e0, 0.016780335911895615e0, -0.47265430073583137e0,
             -0.36743050635638536e0, -0.37497975943257994e0, -0.5271513701209063e0, -0.5619476008133676e0,
             0.8310272499930922e0, 0.008441787366346746e0, 0.7416283853922676e0, -0.33711422046237366e0,
             0.8867308905978669e0, -0.013964733208345414e0, -0.9169834356683927e0, -0.6446160201656775e0,
             -0.4314128159761131e0, -0.07046165367728925e0, 0.1292761369854262e0, 0.7491641153779708e0,
             0.2455752112969821e0, 0.36332977914202347e0, -0.6094132213389469e0, 0.32729417590634935e0,
             -0.5309149553373499e0, -0.6415609049334897e0, 0.04143095385286538e0, 0.6502635955388084e0,
             0.0901308705469112e0, 0.5043457200542485e0, 0.09899258562919144e0, 0.5187021857147951e0,
             0.25440338506733157e0, 0.1650644818795948e0, -0.4873356330199755e0, -0.2184939875223114e0,
             0.4133272998756739e0, 0.4059714749698231e0, 0.22717650761798017e0, 0.6145249525674683e0,
             0.10831417283199905e0, 0.394994014273385e0, 0.7655289572064066e0, 0.7778936372387174e0,
             -0.2005664926106228e0, -0.053457146721091475e0, -0.17667793168009938e0, 0.8080448147049839e0;
  ok ([&&] (gather for ^10 X ^10 -> ($i, $j) { take $A[$i;$j] } Z≅ @sol1)), 'LU-decomp';

  my $x = LU-solve($A, $p, $b);
  my @sol2 = 0.7109223850796318e0, 0.6945023614807339e0, -0.6428828715661831e0, -1.367206776903287e0,
             0.4990441110612922e0, 1.648030263368056e0, 1.164745298345992e0, -0.3278943304709725e0,
             -1.2707664081909453e0, 0.2671149300169611e0;
  ok ([&&] ($x[^10] Z≅ @sol2)), 'LU-solve result';

  ok LU-refine($m, $A, $p, $b, $x) == GSL_SUCCESS, 'LU-refine returns success';
  my @sol3 = 0.7109223850796313e0, 0.6945023614807331e0, -0.6428828715661833e0, -1.3672067769032865e0,
             0.4990441110612925e0, 1.648030263368055e0, 1.164745298345991e0, -0.32789433047097205e0,
             -1.2707664081909444e0, 0.267114930016962e0;
  ok ([&&] ($x[^10] Z≅ @sol3)), 'LU-refine result';

  my Math::Libgsl::Matrix $inv = LU-invert($A, $p);
  my @sol4 = 0.4657186593446201e0, 0.11842664909095556e0, -0.5745520842328259e0, 0.997700202302949e0,
             -0.30150233298887197e0, -0.6502525913384113e0, -0.025583019469254797e0, -0.7327485270419323e0,
             -9.08379551210322e-06, 0.8485588014286635e0, -0.9310372141033609e0, 1.1122430561206031e0,
             0.23585639043945697e0, 0.6501744194139563e0, 1.0669209550911558e0, -1.498604584423743e0,
             -0.05688460319524452e0, -0.28761444402419967e0, -0.30182618992326277e0, 0.7519195174315875e0,
             -0.13975747784841086e0, 0.01744185264994647e0, -0.5444190245066698e0, -0.029226126200458367e0,
             -0.4471418078651977e0, 0.8705850423347158e0, 1.520506066650875e0, -0.6055026555588867e0,
             0.20731715359178657e0, -0.7683510643940016e0, -0.484714343280497e0, 1.1755080551007429e0,
             0.642806659841977e0, -1.6860302079181013e0, 0.6781996088727353e0, -0.4109539378113657e0,
             -1.3805455591370968e0, 1.704716841008282e0, 0.9409654059256847e0, -1.2149762263238228e0,
             0.34600001661870816e0, -0.6895421300276569e0, 0.014139715544870491e0, -0.37145596916822554e0,
             -1.538819567095952e0, 2.6527134833906554e0, 0.3478014142864313e0, -0.29587574200213185e0,
             -0.195436640672239e0, 0.0196669127817628e0, 0.27078108304098536e0, -1.1578991764481668e0,
             1.107664806488947e0, 1.3533661079208208e0, -0.09135869820058927e0, -0.7765895313980696e0,
             0.4557563412237237e0, -1.9892748952664914e0, -0.32595157232760663e0, 1.4620773430200138e0,
             0.48668664311866316e0, 0.21024191520950927e0, -1.0041429818187952e0, 2.302410718410038e0,
             0.7908229311116447e0, -0.34981371745271655e0, 1.3063769635169724e0, -1.309886770937435e0,
             -1.379187927769735e0, -0.5928607463246408e0, 0.8328205218033751e0, -0.5159432461318944e0,
             0.6043838183630114e0, -1.5468488442192012e0, -1.4163123335297516e0, 1.175845404727994e0,
             0.29945121999932e0, 0.36398619292948203e0, -0.33872276887925384e0, 0.22011732067763384e0,
             -0.8754841521716753e0, 0.4599457909478069e0, -0.4940608043808001e0, -0.3168080383371626e0,
             0.41672017513406745e0, -1.2155308366176127e0, -0.9121467311067085e0, 2.2046326806221974e0,
             1.0518442795965695e0, -0.09879718670995495e0, 0.48673022035417246e0, -1.0302669617632596e0,
             0.4310057899199748e0, -1.284509643905147e0, 0.11392958762402547e0, 1.2375551229359707e0,
             -0.4999433355027437e0, 0.21864867946043842e0, 0.41525430495842236e0, 0.1901304717182375e0;
  ok ([&&] (gather for ^10 X ^10 -> ($i, $j) { take $inv[$i;$j] } Z≅ @sol4)), 'LU-invert';

  is-approx LU-det($A, $signum), 0.017475253193805475, 'LU-det';
  is-approx LU-lndet($A), -4.0469695020513194, 'LU-lndet';
  ok LU-sgndet($A, $signum) == 1, 'LU-sgndet';
}

subtest 'QR decomposition' => {
  my Math::Libgsl::Matrix $qr    .= new(3, 5);
  create-general-matrix($qr);

  my $*TOLERANCE = 1e-12;
  my $tau = QR-decomp($qr);
  my @sol1 = -1.1666666666666667e0, -0.6428571428571428e0, -0.4499999999999999e0, -0.3476190476190476e0,
             -0.28367346938775506e0, 0.23076923076923073e0, -0.10171433030139929e0, -0.10533703247651757e0,
             -0.09697695053393682e0, -0.08758181082703652e0, 0.1538461538461538e0, 0.43298825542133207e0,
             0.0039013715732043597e0, 0.005852057359806512e0, 0.006688065554064614e0;
  ok ([&&] (gather for ^3 X ^5 -> ($i, $j) { take $qr[$i;$j] } Z≅ @sol1)), 'QR-decomp: Householder vectors';
  ok ([&&] ($tau[^3] Z≅ (1.8571428571428572e0, 1.6842405528389173e0, 0e0))), 'QR-decomp: Householder coefficients';
  my ($q, $r) = QR-unpack($qr, $tau);
  my @sol2 = -0.8571428571428572e0, 0.5016049165548453e0, 0.11704114719613047e0, -0.4285714285714285e0,
             -0.5684855720954914e0, -0.7022468831767833e0, -0.28571428571428564e0, -0.6520863915212988e0,
             0.7022468831767835e0;
  ok ([&&] (gather for ^3 X ^3 -> ($i, $j) { take $q[$i;$j] } Z≅ @sol2)), 'Q matrix';
  my @sol3 = -1.1666666666666667e0, -0.6428571428571428e0, -0.4499999999999999e0, -0.3476190476190476e0,
             -0.28367346938775506e0, 0e0, -0.10171433030139929e0, -0.10533703247651757e0, -0.09697695053393682e0,
             -0.08758181082703652e0, 0e0, 0e0, 0.0039013715732043597e0, 0.005852057359806512e0, 0.006688065554064614e0;
  ok ([&&] (gather for ^3 X ^5 -> ($i, $j) { take $r[$i;$j] } Z≅ @sol3)), 'R matrix';


  my Math::Libgsl::Matrix $hqr   .= new(2, 2);
  my Math::Libgsl::Vector $rhs   .= new(2);
  my Math::Libgsl::Vector $hsol  .= new(2);
  $rhs[$_] = $_ + 1e0 for ^2;
  $hsol[0] = -8e0; $hsol[1] = 18e0;
  create-general-matrix($hqr);

  my $hd = QR-decomp($hqr);
  my $hx = QR-solve($hqr, $hd, $rhs);
  ok ([&&] ($hx[^2] Z≅ $hsol[^2])), 'QR-solve result';

  my Math::Libgsl::Matrix $m53   .= new(5, 3);
  my Math::Libgsl::Vector $lsrhs .= new(5);
  my Math::Libgsl::Vector $lsr   .= new(5);
  my Math::Libgsl::Vector $lssol .= new(3);
  $lsrhs[$_] = $_ + 1e0 for ^5;
  $lssol[0] = 52.5992295702070e0;
  $lssol[1] = -337.7263113752073e0;
  $lssol[2] = 351.8823436427604e0;
  create-general-matrix($m53);

  my $lsd = QR-decomp($m53);
  my ($lsx, $lsres) = QR-lssolve($m53, $lsd, $lsrhs);
  ok ([&&] ($lsx[^3] Z≅ $lssol[^3])), 'QR-lssolve';

  my Math::Libgsl::Matrix $A     .= new(10, 10);
  my Math::Libgsl::Matrix $B     .= new(10, 5);
  my Math::Libgsl::Random $rnd   .= new;
  create-random-matrix($A, $rnd);
  create-random-matrix($B, $rnd);
  my @sol4[10; 5] =
    [-1.208620986625094, -1.054647224086