context("Cholesky solves")

set.seed(1, kind="Mersenne-Twister", normal.kind="Inversion")

test_that("Cholesky for ddiMatrix works", {
  n <- 100L
  M <- Diagonal(x=runif(n, 0.1, 1))
  y0 <- rnorm(n)
  y1 <- matrix(rnorm(3*n), n, 3)
  y2 <- rsparsematrix(n, 10, 0.1)
  ch <- build_chol(M)
  M <- Diagonal(x = 0.1 + 10*rnorm(n)^2)
  ch$update(M)
  expect_equal(ch$solve(y0), y0/M@x)
  expect_equal(ch$solve(y1), y1/M@x)
  expect_equal(ch$solve(y2), Diagonal(x=1/M@x) %*% y2)
  expect_equal(ch$solve(y0, system="Lt"), y0/sqrt(M@x))
  expect_equal(ch$solve(y0, system="Lt", systemP=TRUE), y0/sqrt(M@x))
  expect_equal(ch$solve(y0, system="L"), y0/sqrt(M@x))
  expect_equal(ch$solve(y1, system="Lt"), y1/sqrt(M@x))
  expect_equal(ch$solve(y2, system="Lt"), Diagonal(x=1/sqrt(M@x)) %*% y2)
})

test_that("Cholesky for matrix works", {
  n <- 12
  M <- crossprod(matrix(rnorm(n^2), n))
  diag(M) <- diag(M) + runif(n, 0.1, 1)
  ch <- build_chol(M)
  M <- crossprod(matrix(rnorm(n^2), n))
  diag(M) <- diag(M) + runif(n, 0.1, 1)
  ch$update(M)
  y0 <- rnorm(n)
  y1 <- matrix(rnorm(3*n), n, 3)
  expect_equal(ch$solve(y0), solve(M, y0))
  expect_equal(ch$solve(y1), solve(M, y1))
  expect_equal(ch$solve(y0, system="Lt"), solve(chol(M), y0))
  expect_equal(ch$solve(y0, system="L"), solve(t(chol(M)), y0))
  expect_equal(ch$solve(y1, system="Lt"), solve(chol(M), y1))
  x <- rnorm(n)
  expect_equal(ch$Ltimes(x), ch$cholM %m*v% x)
})

test_that("Cholesky for dsCMatrix works", {
  n <- 100
  M <- crossprod(rsparsematrix(n, n, density=0.01)) + Diagonal(n)
  ch <- build_chol(M)
  # in-place update only works for same sparsity pattern!!
  M <- M + Diagonal(n)
  ch$update(M)
  cholM <- Cholesky_dsC(M, perm=FALSE, super=NA)
  y0 <- rnorm(n)
  y1 <- matrix(rnorm(3*n), n, 3)
  y2 <- rsparsematrix(n, 10, 0.1)
  expect_equal(ch$solve(y0), as.vector(solve(cholM, y0)))
  expect_equal(ch$solve(y1), array(solve(cholM, y1)@x, dim(y1)))
  expect_equal(ch$solve(y2), solve(cholM, y2))
  expect_equal(ch$solve(y0, system="Lt"), as.vector(solve(cholM, y0, system="Lt")))
  expect_equal(ch$solve(y1, system="Lt"), array(solve(cholM, y1, system="Lt")@x, dim(y1)))
  expect_equal(ch$solve(y2, system="Lt"), solve(cholM, y2, system="Lt"))
  expect_equal(ch$solve(y0, system="L"), as.vector(solve(cholM, y0, system="L")))
  # NB disable following test; currently cannot control cholmod ordering
  #expect_equal(ch$solve(y0, system="L", systemP=TRUE), as.vector(solve(cholM, y0, system="L")))
  expect_equal(ch$solve(y1, system="L"), array(solve(cholM, y1, system="L")@x, dim(y1)))
  expect_equal(ch$solve(y2, system="L"), solve(cholM, y2, system="L"))
  # with permutation
  ch <- build_chol(M, control=chol_control(perm=TRUE))
  M <- M + Diagonal(x=runif(n))
  ch$update(M)
  cholM <- Cholesky_dsC(M, perm=TRUE, super=NA)
  expect_equal(ch$solve(y0), as.vector(solve(cholM, y0)))
  expect_equal(ch$solve(y1), array(solve(cholM, y1)@x, dim(y1)))
  expect_equal(ch$solve(y2), solve(cholM, y2))
  expect_equal(ch$solve(y0, system="Lt"), as.vector(solve(cholM, y0, system="Lt")))
  expect_equal(ch$solve(y1, system="Lt"), array(solve(cholM, y1, system="Lt")@x, dim(y1)))
  expect_equal(ch$solve(y2, system="Lt"), solve(cholM, y2, system="Lt"))
  expect_equal(ch$solve(y0, system="L"), as.vector(solve(cholM, y0, system="L")))
  expect_equal(ch$solve(y1, system="L"), array(solve(cholM, y1, system="L")@x, dim(y1)))
  expect_equal(ch$solve(y2, system="L"), solve(cholM, y2, system="L"))
})