### ----------------------------------------------------
### -- Math::Libgsl::LinearAlgebra
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/03-complex-linalg.rakutest
### ----------------------------------------------------
### -- Chunk 2 of 2
i,
    0.12790443345107724+0.30728043615672507i, -0.41382450080227173-0.2998641492064887i,
    0.4623993173835938+0.10501610829363456i, -0.1302183840605556-0.5391899849671723i,
    -0.26064320325805973-0.2336074980358772i, -0.44456074077415064-0.37464203257111206i,
    0.35063319526326-0.4731562134528027i, 0.10722005094747235+0.002689200302496455i,
    0.038970302790854845-0.24946202471045842i, 0.5718200105099123+0.27775047551860127i,
    0.9710901254280145-0.07050437150345107i, 0.8306187753097398+0.12161422917137432i;
  }

  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 ($sign, $p) = LU-cdecomp($A);
  ok ([&&] (gather for ^10 X ^10 -> ($i, $j) { take $A[$i;$j] }) Z≅ @res)
      &&
     ($sign == -1), 'LU-cdecomp';

  my $x = LU-csolve($A, $p, $b);
  ok ([&&] ($x[^10] Z≅ (-0.0998659594782898-0.36800875303814623i, 0.24869975073113937-0.09427842647983747i, -0.6927476660595517-0.5199882989184399i, -0.728372830967316+0.5084291163037656i, -0.7074957220025232-0.10164197708254931i, 0.7976172708839915-0.12091188140180395i, 0.7940740592055655+0.05203899306549519i, 0.2626638216803528+0.03408545618008574i, 0.39181503654427796-0.4760089055993606i, 0.5595743396086477+0.6508380692999388i))), 'LU-csolve';

  LU-crefine($m, $A, $p, $b, $x);
  ok ([&&] ($x[^10] Z≅ (-0.09986595947828966-0.36800875303814673i, 0.24869975073114026-0.0942784264798372i, -0.6927476660595523-0.5199882989184399i, -0.7283728309673158+0.5084291163037654i, -0.707495722002523-0.10164197708254923i, 0.7976172708839918-0.12091188140180327i, 0.7940740592055656+0.05203899306549554i, 0.2626638216803536+0.034085456180084586i, 0.3918150365442767-0.47600890559936093i, 0.5595743396086471+0.6508380692999393i))), 'LU-crefine result';
}

subtest 'Cholesky Decomposition' => {
  my $*TOLERANCE = 1e-12;
  my Complex @data = (
    66e0 + 0e0i, 0e0 + 64e0i, 126e0 + 63e0i, 124e0 - 62e0i, 61e0 - 61e0i, 60e0 + 60e0i, 0e0 - 59e0i,
    0e0 - 64e0i, 65e0 + 0e0i, 62e0 - 124e0i, -61e0 - 122e0i, -60e0 - 60e0i, 59e0 - 59e0i, -58e0 + 0e0i,
    126e0 - 63e0i, 62e0 + 124e0i, 308e0 + 0e0i, 180e0 - 240e0i, 59e0 - 177e0i, 174e0 + 58e0i, -57e0 - 114e0i,
    124e0 + 62e0i, -61e0 + 122e0i, 180e0 + 240e0i, 299e0 + 0e0i, 174e0 - 58e0i, 57e0 + 171e0i, 56e0 - 112e0i,
    61e0 + 61e0i, -60e0 + 60e0i, 59e0 + 177e0i, 174e0 + 58e0i, 119e0 + 0e0i, 0e0 + 112e0i, 55e0 - 55e0i,
    60e0 - 60e0i, 59e0 + 59e0i, 174e0 - 58e0i, 57e0 - 171e0i, 0e0 - 112e0i, 116e0 + 0e0i, -54e0 - 54e0i,
    0e0 + 59e0i, -58e0 + 0e0i, -57e0 + 114e0i, 56e0 + 112e0i, 55e0 + 55e0i, -54e0 + 54e0i, 60e0 + 0e0i
  );
  my Complex @sol = -0.524944196428570 + 0.209123883928571i,
                     1.052873883928572 + 0.712444196428571i,
                     0.117568824404762 + 0.443191964285714i,
                     0.412862723214286 - 0.356696428571429i,
                     0.815931919642858 - 0.265820312500000i,
                     0.777929687500000 + 0.119484747023810i,
                     1.058733258928571 - 0.132087053571429i;
  my Math::Libgsl::Matrix::Complex64 $m .= new: 7, 7;
  for ^7 X ^7 -> ($i, $j) {
    $m.set($i, $j, @data[$i * 7 + $j]);
  }
  my Math::Libgsl::Vector::Complex64 $sol .= new(7);
  $sol.set($_, @sol[$_]) for ^7;
  my Math::Libgsl::Vector::Complex64 $rhs .= new(7);
  $rhs.set($_, $_ + 1 + 0i) for ^7;

  ok cholesky-cdecomp($m) == GSL_SUCCESS, 'cholesky-cdecomp';
  my $x = cholesky-csolve($m, $rhs);
  ok ([&&] ($x[^7] Z≅ $sol[^7])), 'cholesky-csolve';

  my Math::Libgsl::Matrix::Complex64 $m1 .= new(2, 2);
  $m1[0;0] = 92.303+0i; $m1[0;1] = 10.858+1.798i; $m1[1;0] = 10.858-1.798i; $m1[1;1] = 89.027+0i;
  my Math::Libgsl::Matrix::Complex64 $inv .= new(2, 2);
  $inv[0;0] = 0.010995969489615481+0i; $inv[0;1] = -0.0013411014267384605-0.00022207592238678868i; $inv[1;0] = -0.0013411014267384605+0.00022207592238678868i; $inv[1;1] = 0.011400597254765159+0i;

  cholesky-cdecomp($m1);
  cholesky-cinvert($m1);
  ok ([&&] (gather for ^2 X ^2 -> ($i, $j) { take $m1[$i;$j] } Z≅ gather for ^2 X ^2 -> ($i, $j) { take $inv[$i;$j] })), 'cholesky-cinvert result';
}

done-testing;


