### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 4 of 4
($rhs);
  gsl_vector_free($x);
  gsl_vector_free($sol);
  return @res;
}

sub tridiag-cyc-solve
{
  my $M = 3;
  my gsl_vector $abovediag = gsl_vector_calloc($M);
  my gsl_vector $belowdiag = gsl_vector_calloc($M);
  my gsl_vector $diag      = gsl_vector_calloc($M);
  my gsl_vector $rhs       = gsl_vector_calloc($M);
  my gsl_vector $x         = gsl_vector_calloc($M);
  my gsl_vector $sol       = gsl_vector_calloc($M);
  my Bool @res;
  gsl_vector_set($sol, 0, (3/2).Num);
  gsl_vector_set($sol, 1, -.5e0);
  gsl_vector_set($sol, 2, .5e0);
  for ^$M -> $i {
    gsl_vector_set($diag,      $i, 1e0);
    gsl_vector_set($rhs,       $i, $i + 1e0);
    gsl_vector_set($abovediag, $i, 2e0);
    gsl_vector_set($belowdiag, $i, 1e0);
  }

  @res.push: gsl_linalg_solve_cyc_tridiag($diag, $abovediag, $belowdiag, $rhs, $x) == GSL_SUCCESS;
  @res.push: test-num64-vector($x, $sol);

  gsl_vector_free($abovediag);
  gsl_vector_free($belowdiag);
  gsl_vector_free($diag);
  gsl_vector_free($rhs);
  gsl_vector_free($x);
  gsl_vector_free($sol);
  return @res;
}

sub triang-invert
{
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $M = 3;
  my gsl_matrix $T   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $B   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $sol = gsl_matrix_calloc($M, $M);
  gsl_matrix_set_identity($sol);
  for ^$M -> $i {
    for 0..$i -> $j {
      my num64 $mij = gsl_rng_uniform($r);
      $mij = max($mij, .3e0) if $i == $j;
      gsl_matrix_set($T, $i, $j, $mij);
    }
  }
  gsl_matrix_memcpy($B, $T);
  gsl_linalg_tri_lower_invert($B);
  gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1e0, $T, $B);

  my $res = test-num64-matrix($B, $sol);

  gsl_rng_free($r);
  gsl_matrix_free($T);
  gsl_matrix_free($B);
  gsl_matrix_free($sol);
  return $res;
}

sub symm-to-band(gsl_matrix $m, gsl_matrix $u)
{
  my $M = $m.size1;
  my gsl_vector_view $view1 = alloc_gsl_vector_view();
  my gsl_vector_view $view2 = alloc_gsl_vector_view();
  for ^$M {
    my gsl_vector $vv1 = mgsl_matrix_subdiagonal($view1, $m, $_);
    my gsl_vector $vv2 = mgsl_matrix_subcolumn($view2, $u, $_, 0, $M - $_);
    gsl_vector_memcpy($vv2, $vv1);
  }
  free_gsl_vector_view($view1);
  free_gsl_vector_view($view2);
}

sub cholesky-band-decomp
{
  my $M = 3;
  my $*TOLERANCE = 1e-12;
  my gsl_rng $r       = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $v    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $l    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $lt   = gsl_matrix_calloc($M, $M);
  my gsl_vector $work = gsl_vector_alloc(3 * $M);
  my Bool @ret;

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

  gsl_linalg_cholesky_band_decomp($v);
  gsl_linalg_cholesky_band_unpack($v, $l);

  gsl_matrix_transpose_tricpy(CblasLower, CblasNonUnit, $lt, $l);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $l, $lt, 0e0, $a);

  @ret.push: test-num64-matrix($m, $a);

  my num64 ($rcond, $rcond_expected);
  gsl_matrix_memcpy($a, $m);
  gsl_linalg_cholesky_decomp1($a);
  gsl_linalg_cholesky_rcond($a, $rcond_expected, $work);
  gsl_linalg_cholesky_band_rcond($v, $rcond, $work);
  @ret.push: $rcond ≅ $rcond_expected;

  gsl_matrix_free($m);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($l);
  gsl_matrix_free($lt);
  gsl_vector_free($work);
  gsl_rng_free($r);
  return @ret;
}

sub cholesky-band-solve
{
  my $M = 3;
  my gsl_rng $r      = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $u   = gsl_matrix_calloc($M, $M);
  my gsl_vector $rhs = gsl_vector_calloc($M);
  my gsl_vector $sol = gsl_vector_calloc($M);
  my gsl_vector $x   = gsl_vector_calloc($M);

  create-posdef-band-matrix($m);
  create-random-vector($rhs, $r);
  gsl_blas_dsymv(CblasLower, 1e0, $m, $sol, 0e0, $rhs);
  symm-to-band($m, $u);

  gsl_linalg_cholesky_band_decomp($u);
  gsl_linalg_cholesky_band_solve($u, $rhs, $x);

  my $ret = test-num64-vector($x, $sol);

  gsl_matrix_free($m);
  gsl_matrix_free($u);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_vector_free($x);
  gsl_rng_free($r);
  return $ret;
}

sub cholesky-band-invert
{
  my $M = 3;
  my gsl_rng $r       = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $v    = gsl_matrix_calloc($M, $M);
  my gsl_matrix $minv = gsl_matrix_calloc($M, $M);
  my gsl_matrix $c    = gsl_matrix_calloc($M, $M);

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

  gsl_linalg_cholesky_band_decomp($v);
  gsl_linalg_cholesky_band_invert($v, $minv);

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

  my $*TOLERANCE = 1e-12;
  my $ret = True;
  for ^$M X ^$M -> ($i, $j) {
    $ret &&= gsl_matrix_get($c, $i, $j) ≅ ($i == $j ?? 1e0 !! 0e0);
  }

  gsl_matrix_free($m);
  gsl_matrix_free($v);
  gsl_matrix_free($minv);
  gsl_matrix_free($c);
  gsl_rng_free($r);
  return $ret;
}

sub ldlt-band-decomp
{
  my $M = 3;
  my gsl_matrix $m  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $v  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $a  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $l  = gsl_matrix_calloc($M, $M);
  my gsl_matrix $lt = gsl_matrix_calloc($M, $M);
  my gsl_vector_view $d = alloc_gsl_vector_view();

  create-posdef-band-matrix($m);
  gsl_matrix_memcpy($v, $m);
  gsl_linalg_ldlt_decomp($v);
  gsl_matrix_tricpy(CblasLower, CblasUnit, $l, $v);
  my $diag = mgsl_matrix_diagonal($d, $l);

  gsl_vector_set_all($diag, 1e0);
  gsl_matrix_transpose_tricpy(CblasLower, CblasNonUnit, $lt, $l);
  $diag = mgsl_matrix_diagonal($d, $v);
  for ^$M -> $i {
    my gsl_vector_view $c = alloc_gsl_vector_view();
    my $col = mgsl_matrix_column($c, $l, $i);
    my $di = gsl_vector_get($diag, $i);
    gsl_vector_scale($col, $di);
  }
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1e0, $l, $lt, 0e0, $a);

  my $ret = test-num64-matrix($m, $a);

  gsl_matrix_free($m);
  gsl_matrix_free($v);
  gsl_matrix_free($a);
  gsl_matrix_free($l);
  gsl_matrix_free($lt);
  free_gsl_vector_view($d);
  return $ret;
}

sub ldlt-band-solve
{
  my $M = 3;
  my gsl_rng $r      = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m   = gsl_matrix_calloc($M, $M);
  my gsl_matrix $u   = gsl_matrix_calloc($M, $M);
  my gsl_vector $rhs = gsl_vector_calloc($M);
  my gsl_vector $sol = gsl_vector_calloc($M);
  my gsl_vector $x   = gsl_vector_calloc($M);

  create-posdef-band-matrix($m);
  create-random-vector($sol, $r);
  gsl_blas_dsymv(CblasLower, 1e0, $m, $sol, 0e0, $rhs);
  gsl_matrix_memcpy($u, $m);

  gsl_linalg_ldlt_decomp($u);
  gsl_linalg_ldlt_solve($u, $rhs, $x);

  my $ret = test-num64-vector($x, $sol);

  gsl_matrix_free($m);
  gsl_matrix_free($u);
  gsl_vector_free($rhs);
  gsl_vector_free($sol);
  gsl_vector_free($x);
  gsl_rng_free($r);
  return $ret;
}

subtest 'LU decomposition' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m = gsl_matrix_calloc(10, 10);
  ok create-random-matrix($m, $r) == GSL_SUCCESS, 'can create a random matrix';
  my gsl_matrix_complex $cm = gsl_matrix_complex_alloc(10, 10);
  ok create-random-complex-matrix($cm, $r) == GSL_SUCCESS, 'can create a random complex matrix';
  my $*TOLERANCE = 1e-12;

  ok LU-decomp-num64($m), 'LU num64 decomposition';
  is-deeply LU-solve-num64($m, $r), Array[Bool].new(True xx 5), 'LU num64 solve';
  is-deeply LU-solve-complex64($cm, $r), Array[Bool].new(True xx 5), 'LU complex decomp and solve';
  is-deeply LU-invert-num64($m), Array[Bool].new(True xx 2), 'LU num64 invert';
  is-deeply LU-invert-complex64($cm), Array[Bool].new(True xx 2), 'LU complex invert';
  is-approx LU-det-num64($m), 0.017475253193805475, 'num64 determinant';
  my gsl_complex $cres = alloc_gsl_complex;
  LU-det-complex64($cm, $cres);
  ok $cres.dat[0] ≅ 1.2669771918585264 && $cres.dat[1] ≅ -0.4448915095767637, 'complex determinant';
  free_gsl_complex($cres);
  is-approx gsl_linalg_LU_lndet($m), -8.682022026520569, 'num64 ln determinant';
  is-approx gsl_linalg_complex_LU_lndet($cm), -1.3441985858994467, 'complex ln determinant';
  ok LU-sgndet-num64($m) == 1, 'num64 sign of determinant';
  is-approx LU-sgndet-complex64($cm), 0.9435212924956552-0.33131189324762783i, 'complex phase of determinant';

  gsl_matrix_free($m);
  gsl_matrix_complex_free($cm);
  gsl_rng_free($r);
}

subtest 'QR decomposition' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my gsl_matrix $m35 = create-general-matrix(3, 5);
  my gsl_matrix $m53 = create-general-matrix(5, 3);
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);
  my gsl_matrix $A = gsl_matrix_calloc(10, 10);
  my gsl_matrix $S = gsl_matrix_calloc(10, 10);
  my gsl_matrix $B = gsl_matrix_calloc(10, 5);
  create-random-matrix($A, $r);
  create-random-matrix($B, $r);
  create-random-matrix($S, $r);

  is-deeply QR-decomp($m35),      Array[Bool].new(True xx 4), 'QR decomposition';
  is-deeply QR-solve($hilb2),     Array[Bool].new(True xx 3), 'QR solve';
  is-deeply QR-lssolve($m53),     Array[Bool].new(True xx 3), 'QR least square solve';
  is-deeply QR-QTmat($A, $B),     Array[Bool].new(True xx 3), 'QR QTmat';
  is-deeply QR-QRsolve($hilb2),   Array[Bool].new(True xx 4), 'QR QRsolve';
  is-deeply QR-update($m35),      Array[Bool].new(True xx 4), 'QR update';

  gsl_matrix_free($S);
  gsl_matrix_free($A);
  gsl_matrix_free($B);
  gsl_matrix_free($hilb2);
  gsl_matrix_free($m53);
  gsl_matrix_free($m35);
}

subtest 'QR decomposition with column pivoting' => {
  my gsl_matrix $m53 = create-general-matrix(5, 3);
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);

  is-deeply QRPT-decomp($m53), Array[Bool].new(True xx 5), 'QRPT decomp';
  is-deeply QRPT-decomp2($hilb2), Array[Bool].new(True xx 3), 'QRPT decomp2';
  is-deeply QRPT-solve($hilb2), Array[Bool].new(True xx 3), 'QRPT solve';
  is-deeply QRPT-lssolve($hilb2), Array[Bool].new(True xx 4), 'QRPT lssolve';
  is-deeply QRPT-lssolve2($hilb2), Array[Bool].new(True xx 4), 'QRPT lssolve2';
  is-deeply QRPT-update($hilb2), Array[Bool].new(True xx 4), 'QRPT update';

  gsl_matrix_free($hilb2);
  gsl_matrix_free($m53);
}

subtest 'LQ decomposition' => {
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);

  is-deeply LQ-solve($hilb2), Array[Bool].new(True xx 3), 'LQ decomp';
  is-deeply LQ-LQsolve($hilb2), Array[Bool].new(True xx 4), 'LQ solve';
  is-deeply LQ-lssolve($hilb2), Array[Bool].new(True xx 4), 'LQ lssolve';
  is-deeply LQ-decomp($hilb2), Array[Bool].new(True xx 3), 'LQ lssolve';

  gsl_matrix_free($hilb2);
}

subtest 'Complete Orthogonal Decomposition' => {
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);

  is-deeply COD-decomp($hilb2), Array[Bool].new(True xx 3), 'COD decomp';
  is-deeply COD-lssolve($hilb2), Array[Bool].new(True xx 4), 'COD lssolve';
  is-deeply COD-lssolve2, Array[Bool].new(True xx 4), 'COD lssolve2';

  gsl_matrix_free($hilb2);
}

subtest 'Singular Value Decomposition' => {
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);

  is-deeply SV-decomp($hilb2), Array[Bool].new(True xx 2), 'SV decomp';
  is-deeply SV-decomp-jacobi($hilb2), Array[Bool].new(True xx 2), 'SV decomp_jacobi';
  is-deeply SV-decomp-mod($hilb2), Array[Bool].new(True xx 2), 'SV decomp_mod';
  is-deeply SV-solve($hilb2), Array[Bool].new(True xx 3), 'SV solve';

  gsl_matrix_free($hilb2);
}

subtest 'Cholesky Decomposition' => {
  my gsl_matrix $hilb2 = create-hilbert-matrix(2);

  is-deeply Cholesky-decomp($hilb2), Array[Bool].new(True xx 3), 'Cholesky num64 decomp';
  is-deeply Cholesky-invert($hilb2), Array[Bool].new(True xx 3), 'Cholesky num64 invert';
  is-deeply Cholesky-solve($hilb2), Array[Bool].new(True xx 3), 'Cholesky num64 solve';
  is-deeply Cholesky-solve2($hilb2), Array[Bool].new(True xx 3), 'scaled Cholesky solve';
  is-deeply Cholesky-decomp-unit($hilb2), Array[Bool].new(True xx 2), 'Cholesky decomp unit';
  is-deeply Cholesky-complex-solve, Array[Bool].new(True xx 3), 'Cholesky complex64 solve';
  is-deeply Cholesky-complex-decomp, Array[Bool].new(True xx 19), 'Cholesky complex64 decomp';
  is-deeply Cholesky-complex-invert, Array[Bool].new(True xx 3), 'Cholesky complex64 invert';
  is-deeply mCholesky-decomp, Array[Bool].new(True xx 3), 'modified Cholesky decomp';
  is-deeply mCholesky-solve, Array[Bool].new(True xx 3), 'modified Cholesky solve';
  is-deeply mCholesky-invert, Array[Bool].new(True xx 3), 'modified Cholesky invert';
  is-deeply pCholesky-decomp($hilb2), Array[Bool].new(True xx 4), 'pivoted Cholesky decomp';
  is-deeply pCholesky-solve($hilb2), Array[Bool].new(True xx 3), 'pivoted Cholesky solve';
  is-deeply pCholesky-invert($hilb2), Array[Bool].new(True xx 3), 'pivoted Cholesky invert';

  gsl_matrix_free($hilb2);
}

my gsl_matrix $hilb2 = create-hilbert-matrix(2);
is-deeply bidiag-decomp($hilb2), Array[Bool].new(True xx 3), 'bidiag decomp';
is-deeply HH-solve($hilb2), Array[Bool].new(True xx 2), 'Householder solver';
gsl_matrix_free($hilb2);

subtest 'tridiagonal systems' => {
  is-deeply tridiag-symm-solve, Array[Bool].new(True xx 2), 'symmetric tridiagonal solver';
  is-deeply tridiag-symm-cyc-solve, Array[Bool].new(True xx 2), 'symmetric cyclic tridiagonal solver';
  is-deeply tridiag-solve, Array[Bool].new(True xx 2), 'tridiagonal solver';
  is-deeply tridiag-cyc-solve, Array[Bool].new(True xx 2), 'cyclic tridiagonal solver';
}

ok triang-invert, 'triangular system inversion';

if $gsl-version > v2.5 {
  subtest 'Banded Cholesky decomposition' => {
    is-deeply cholesky-band-decomp, Array[Bool].new(True xx 2), 'banded Cholesky decomp, unpack, rcond';
    ok cholesky-band-solve,  'banded Cholesky solve';
    ok cholesky-band-invert, 'banded Cholesky invert';
  }
  subtest 'Banded LDLT decomposition' => {
    ok ldlt-band-decomp, 'banded LDLT decomp';
    ok ldlt-band-solve, 'banded LDLT solve';
  }
}

done-testing;


