/* linalg/test_cholesky.c

 *

 * Copyright (C) 2016 Patrick Alken

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_test.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_ieee_utils.h>

#include <gsl/gsl_permute_vector.h>

#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

#include <gsl/gsl_rng.h>

#include <gsl/gsl_permutation.h>

#include "test_common.c"

static int test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,

                                    const double expected_rcond, const double eps,

                                    const char * desc);

static int test_cholesky_decomp(gsl_rng * r);

int test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);

int test_cholesky_invert(gsl_rng * r);

static int test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,

                                     const double expected_rcond, const double eps,

                                     const char * desc);

static int test_pcholesky_decomp(gsl_rng * r);

int test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,

                             const gsl_vector * sol, const double eps,

                             const char * desc);

static int test_pcholesky_solve(gsl_rng * r);

int test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);

int test_pcholesky_invert(gsl_rng * r);

static int test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,

                                     const double expected_rcond, const double eps, const char * desc);

/* Hilbert matrix condition numbers, as calculated by LAPACK DPOSVX */

double hilb_rcond[] = { 1.000000000000e+00, 3.703703703704e-02, 1.336898395722e-03,

                        3.524229074890e-05, 1.059708198754e-06, 3.439939465186e-08,

                        1.015027593823e-09, 2.952221630602e-11, 9.093751565191e-13,

                        2.828277420229e-14, 8.110242564869e-16, 2.409320075800e-17 };

static int

test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,

                         const double expected_rcond, const double eps,

                         const char * desc)

{

  int s = 0;

  size_t i, j, N = m->size2;

  gsl_matrix * V  = gsl_matrix_alloc(N, N);

  gsl_matrix * A  = gsl_matrix_alloc(N, N);

  gsl_matrix * L  = gsl_matrix_calloc(N, N);

  gsl_matrix * LT = gsl_matrix_calloc(N, N);

  gsl_vector * S = gsl_vector_alloc(N);

  gsl_matrix_memcpy(V, m);

  if (scale)

    s += gsl_linalg_cholesky_decomp2(V, S);

  else

    s += gsl_linalg_cholesky_decomp1(V);

  /* compute L and LT */

  gsl_matrix_tricpy('L', 1, L, V);

  gsl_matrix_transpose_tricpy('L', 1, LT, L);

  if (scale)

    {

      /* L <- S^{-1} L, LT <- LT S^{-1} */

      for (i = 0; i < N; ++i)

        {

          double Si = gsl_vector_get(S, i);

          gsl_vector_view v = gsl_matrix_row(L, i);

          gsl_vector_view w = gsl_matrix_column(LT, i);

          gsl_vector_scale(&v.vector, 1.0 / Si);

          gsl_vector_scale(&w.vector, 1.0 / Si);

        }

    }

  /* compute A = L LT */

  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  for (i = 0; i < N; i++)

    {

      for (j = 0; j < N; j++)

        {

          double Aij = gsl_matrix_get(A, i, j);

          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(Aij, mij, eps,

                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, Aij, mij);

        }

    }

  if (expected_rcond > 0 && !scale)

    {

      gsl_vector *work = gsl_vector_alloc(3 * N);

      double rcond;

      gsl_linalg_cholesky_rcond(V, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,

                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",

                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);

    }

  gsl_matrix_free(V);

  gsl_matrix_free(A);

  gsl_matrix_free(L);

  gsl_matrix_free(LT);

  gsl_vector_free(S);

  return s;

}

static int

test_cholesky_decomp(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_cholesky_decomp_eps(0, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp unscaled random");

      test_cholesky_decomp_eps(1, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp scaled random");

      if (N <= 12)

        {

          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)

            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_cholesky_decomp_eps(0, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp unscaled hilbert");

          test_cholesky_decomp_eps(1, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp scaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}

int

test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)

{

  int s = 0;

  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);

  gsl_matrix * c  = gsl_matrix_alloc(N, N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_cholesky_decomp1(v);

  s += gsl_linalg_cholesky_invert(v);

  /* c = m m^{-1} */

  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, v, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)

    {

      for (j = 0; j < N; ++j)

        {

          double cij = gsl_matrix_get(c, i, j);

          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, cij, expected);

        }

    }

  gsl_matrix_free(v);

  gsl_matrix_free(c);

  return s;

}

int

test_cholesky_invert(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_cholesky_invert_eps(m, N * GSL_DBL_EPSILON, "cholesky_invert unscaled random");

      if (N <= 4)

        {

          create_hilbert_matrix2(m);

          test_cholesky_invert_eps(m, N * 256.0 * GSL_DBL_EPSILON, "cholesky_invert unscaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}

static int

test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,

                          const double expected_rcond, const double eps, const char * desc)

{

  int s = 0;

  size_t i, j, N = m->size2;

  gsl_matrix * LDLT = gsl_matrix_alloc(N, N);

  gsl_matrix * V  = gsl_matrix_alloc(N, N);

  gsl_matrix * A  = gsl_matrix_alloc(N, N);

  gsl_matrix * L  = gsl_matrix_alloc(N, N);

  gsl_matrix * LT = gsl_matrix_alloc(N, N);

  gsl_vector * S = gsl_vector_alloc(N);

  gsl_vector * E = gsl_vector_alloc(N);

  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_vector_view D = gsl_matrix_diagonal(LDLT);

  gsl_matrix_memcpy(LDLT, m);

  s += gsl_linalg_mcholesky_decomp(LDLT, perm, E);

  /* check that the upper triangle of LDLT equals original matrix */

  for (i = 0; i < N; ++i)

    {

      for (j = i + 1; j < N; ++j)

        {

          double mij = gsl_matrix_get(m, i, j);

          double aij = gsl_matrix_get(LDLT, i, j);

          gsl_test_rel(aij, mij, 1.0e-12,

                       "%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, aij, mij);

        }

    }

  if (posdef)

    {

      /* ||E|| should be 0 */

      double norm = gsl_blas_dnrm2(E);

      s = norm != 0.0;

      gsl_test(s, "%s: (%zu,%zu): ||E|| = %.12e",

               desc, N, N, norm);

      /* check that D is decreasing */

      s = 0;

      for (i = 1; i < N; ++i)

        {

          double dprev = gsl_vector_get(&D.vector, i - 1);

          double di = gsl_vector_get(&D.vector, i);

          if (di > dprev)

            s = 1;

        }

      gsl_test(s, "%s: (%zu,%zu): D is not decreasing",

               desc, N, N);

    }

  /* compute L and LT */

  gsl_matrix_set_identity(L);

  gsl_matrix_set_identity(LT);

  gsl_matrix_tricpy('L', 0, L, LDLT);

  gsl_matrix_transpose_tricpy('L', 0, LT, L);

  /* compute (L sqrt(D)) and (sqrt(D) LT) */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_column(L, i);

      gsl_vector_view w = gsl_matrix_row(LT, i);

      double di = gsl_vector_get(&D.vector, i);

      gsl_vector_scale(&v.vector, sqrt(di));

      gsl_vector_scale(&w.vector, sqrt(di));

    }

  /* compute A = L D LT */

  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  /* compute V = P (S M S + E) P^T */

  gsl_matrix_memcpy(V, m);

  D = gsl_matrix_diagonal(V);

  /* compute S M S */

  if (scale)

    {

      gsl_linalg_cholesky_scale_apply(V, S);

      gsl_matrix_transpose_tricpy('L', 0, V, V);

    }

  /* compute S M S + E */

  gsl_vector_add(&D.vector, E);

  /* compute M P^T */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_row(V, i);

      gsl_permute_vector(perm, &v.vector);

    }

  /* compute P M P^T */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_column(V, i);

      gsl_permute_vector(perm, &v.vector);

    }

  for (i = 0; i < N; i++)

    {

      double Ei = gsl_vector_get(E, i);

      for (j = 0; j < N; j++)

        {

          double Aij = gsl_matrix_get(A, i, j); /* L D L^T */

          double Bij = gsl_matrix_get(V, i, j); /* P M P^T */

          double Cij;                           /* P M P^T + E */

          if (i == j)

            Cij = Bij + Ei*0;

          else

            Cij = Bij;

          gsl_test_rel(Aij, Cij, eps,

                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, Aij, Cij);

        }

    }

  if (expected_rcond > 0 && !scale)

    {

      gsl_vector *work = gsl_vector_alloc(3 * N);

      double rcond;

      gsl_linalg_mcholesky_rcond(LDLT, perm, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,

                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",

                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);

    }

  gsl_matrix_free(LDLT);

  gsl_matrix_free(V);

  gsl_matrix_free(A);

  gsl_matrix_free(L);

  gsl_matrix_free(LT);

  gsl_vector_free(S);

  gsl_vector_free(E);

  gsl_permutation_free(perm);

  return s;

}

static int

test_mcholesky_decomp(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_mcholesky_decomp_eps(1, 0, m, -1.0, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random posdef");

      create_symm_matrix(m, r);

      test_mcholesky_decomp_eps(0, 0, m, -1.0, 8192.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random symm");

      if (N <= 8)

        {

          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)

            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_mcholesky_decomp_eps(1, 0, m, expected_rcond, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}

int

test_mcholesky_solve_eps(const gsl_matrix * m, const gsl_vector * rhs,

                         const gsl_vector * sol, const double eps,

                         const char * desc)

{

  int s = 0;

  size_t i, N = m->size1;

  gsl_matrix * u  = gsl_matrix_alloc(N, N);

  gsl_vector * x = gsl_vector_calloc(N);

  gsl_vector * S = gsl_vector_alloc(N);

  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(u, m);

  s += gsl_linalg_mcholesky_decomp(u, perm, NULL);

  s += gsl_linalg_mcholesky_solve(u, perm, rhs, x);

  for (i = 0; i < N; i++)

    {

      double xi = gsl_vector_get(x, i);

      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps,

                   "%s: %3lu[%lu]: %22.18g   %22.18g\n",

                   desc, N, i, xi, yi);

    }

  gsl_vector_free(x);

  gsl_vector_free(S);

  gsl_matrix_free(u);

  gsl_permutation_free(perm);

  return s;

}

static int

test_mcholesky_solve(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      gsl_vector * rhs = gsl_vector_alloc(N);

      gsl_vector * sol = gsl_vector_alloc(N);

      create_posdef_matrix(m, r);

      create_random_vector(sol, r);

      gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

      test_mcholesky_solve_eps(m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "mcholesky_solve random");

      if (N <= 3)

        {

          create_hilbert_matrix2(m);

          gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

          test_mcholesky_solve_eps(m, rhs, sol, 1.0e3 * N * GSL_DBL_EPSILON, "mcholesky_solve hilbert");

        }

      gsl_matrix_free(m);

      gsl_vector_free(rhs);

      gsl_vector_free(sol);

    }

  return s;

}

int

test_mcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)

{

  int s = 0;

  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);

  gsl_matrix * c  = gsl_matrix_alloc(N, N);

  gsl_matrix * minv = gsl_matrix_alloc(N, N);

  gsl_vector * E = gsl_vector_alloc(N);

  gsl_permutation * p = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_mcholesky_decomp(v, p, E);

  s += gsl_linalg_mcholesky_invert(v, p, minv);

  /* c = m m^{-1} */

  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)

    {

      for (j = 0; j < N; ++j)

        {

          double cij = gsl_matrix_get(c, i, j);

          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, cij, expected);

        }

    }

  gsl_matrix_free(v);

  gsl_matrix_free(c);

  gsl_matrix_free(minv);

  gsl_vector_free(E);

  gsl_permutation_free(p);

  return s;

}

int

test_mcholesky_invert(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 30;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_mcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "mcholesky_invert unscaled random");

      if (N <= 4)

        {

          create_hilbert_matrix2(m);

          test_mcholesky_invert_eps(m, 256.0 * N * GSL_DBL_EPSILON, "mcholesky_invert unscaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}

static int

test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,

                          const double expected_rcond, const double eps,

                          const char * desc)

{

  int s = 0;

  size_t i, j, N = m->size2;

  gsl_matrix * LDLT = gsl_matrix_alloc(N, N);

  gsl_matrix * V  = gsl_matrix_alloc(N, N);

  gsl_matrix * A  = gsl_matrix_alloc(N, N);

  gsl_matrix * L  = gsl_matrix_alloc(N, N);

  gsl_matrix * LT = gsl_matrix_alloc(N, N);

  gsl_vector * S = gsl_vector_alloc(N);

  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_vector_view D = gsl_matrix_diagonal(LDLT);

  gsl_matrix_memcpy(LDLT, m);

  if (scale)

    s += gsl_linalg_pcholesky_decomp2(LDLT, perm, S);

  else

    s += gsl_linalg_pcholesky_decomp(LDLT, perm);

  /* check that the upper triangle of LDLT equals original matrix */

  for (i = 0; i < N; ++i)

    {

      for (j = i + 1; j < N; ++j)

        {

          double mij = gsl_matrix_get(m, i, j);

          double aij = gsl_matrix_get(LDLT, i, j);

          gsl_test_rel(aij, mij, 1.0e-12,

                       "%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, aij, mij);

        }

    }

  /* check that D is decreasing */

  s = 0;

  for (i = 1; i < N; ++i)

    {

      double dprev = gsl_vector_get(&D.vector, i - 1);

      double di = gsl_vector_get(&D.vector, i);

      if (di > dprev)

        s = 1;

    }

  gsl_test(s, "%s: (%zu,%zu): D is not decreasing",

           desc, N, N);

  /* compute L and LT */

  gsl_matrix_set_identity(L);

  gsl_matrix_set_identity(LT);

  gsl_matrix_tricpy('L', 0, L, LDLT);

  gsl_matrix_transpose_tricpy('L', 0, LT, L);

  /* compute (L sqrt(D)) and (sqrt(D) LT) */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_column(L, i);

      gsl_vector_view w = gsl_matrix_row(LT, i);

      double di = gsl_vector_get(&D.vector, i);

      gsl_vector_scale(&v.vector, sqrt(di));

      gsl_vector_scale(&w.vector, sqrt(di));

    }

  /* compute A = L D LT */

  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);

  /* compute V = P S M S P^T */

  gsl_matrix_memcpy(V, m);

  /* compute S M S */

  if (scale)

    {

      gsl_linalg_cholesky_scale_apply(V, S);

      gsl_matrix_transpose_tricpy('L', 0, V, V);

    }

  /* compute M P^T */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_row(V, i);

      gsl_permute_vector(perm, &v.vector);

    }

  /* compute P M P^T */

  for (i = 0; i < N; ++i)

    {

      gsl_vector_view v = gsl_matrix_column(V, i);

      gsl_permute_vector(perm, &v.vector);

    }

  for (i = 0; i < N; i++)

    {

      for (j = 0; j < N; j++)

        {

          double Aij = gsl_matrix_get(A, i, j); /* L D L^T */

          double Bij = gsl_matrix_get(V, i, j); /* P M P^T */

          gsl_test_rel(Aij, Bij, eps,

                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, Aij, Bij);

        }

    }

  if (expected_rcond > 0 && !scale)

    {

      gsl_vector *work = gsl_vector_alloc(3 * N);

      double rcond;

      gsl_linalg_pcholesky_rcond(LDLT, perm, &rcond, work);

      gsl_test_rel(rcond, expected_rcond, 1.0e-6,

                   "%s rcond: (%3lu,%3lu): %22.18g   %22.18g\n",

                   desc, N, N, rcond, expected_rcond);

      gsl_vector_free(work);

    }

  gsl_matrix_free(LDLT);

  gsl_matrix_free(V);

  gsl_matrix_free(A);

  gsl_matrix_free(L);

  gsl_matrix_free(LT);

  gsl_vector_free(S);

  gsl_permutation_free(perm);

  return s;

}

static int

test_pcholesky_decomp(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_pcholesky_decomp_eps(0, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled random");

      test_pcholesky_decomp_eps(1, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled random");

      if (N <= 12)

        {

          double expected_rcond = -1.0;

          if (hilb_rcond[N - 1] > 1.0e-12)

            expected_rcond = hilb_rcond[N - 1];

          create_hilbert_matrix2(m);

          test_pcholesky_decomp_eps(0, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled hilbert");

          test_pcholesky_decomp_eps(1, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}

int

test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,

                         const gsl_vector * sol, const double eps,

                         const char * desc)

{

  int s = 0;

  size_t i, N = m->size1;

  gsl_matrix * u  = gsl_matrix_alloc(N, N);

  gsl_vector * x = gsl_vector_calloc(N);

  gsl_vector * S = gsl_vector_alloc(N);

  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(u, m);

  if (scale)

    {

      s += gsl_linalg_pcholesky_decomp2(u, perm, S);

      s += gsl_linalg_pcholesky_solve2(u, perm, S, rhs, x);

    }

  else

    {

      s += gsl_linalg_pcholesky_decomp(u, perm);

      s += gsl_linalg_pcholesky_solve(u, perm, rhs, x);

    }

  for (i = 0; i < N; i++)

    {

      double xi = gsl_vector_get(x, i);

      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps,

                   "%s: %3lu[%lu]: %22.18g   %22.18g\n",

                   desc, N, i, xi, yi);

    }

  gsl_vector_free(x);

  gsl_vector_free(S);

  gsl_matrix_free(u);

  gsl_permutation_free(perm);

  return s;

}

static int

test_pcholesky_solve(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 50;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      gsl_vector * rhs = gsl_vector_alloc(N);

      gsl_vector * sol = gsl_vector_alloc(N);

      create_posdef_matrix(m, r);

      create_random_vector(sol, r);

      gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

      test_pcholesky_solve_eps(0, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled random");

      test_pcholesky_solve_eps(1, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled random");

      if (N <= 3)

        {

          create_hilbert_matrix2(m);

          gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);

          test_pcholesky_solve_eps(0, m, rhs, sol, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled hilbert");

          test_pcholesky_solve_eps(1, m, rhs, sol, 2048.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled hilbert");

        }

      gsl_matrix_free(m);

      gsl_vector_free(rhs);

      gsl_vector_free(sol);

    }

  return s;

}

int

test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)

{

  int s = 0;

  size_t i, j, N = m->size1;

  gsl_matrix * v  = gsl_matrix_alloc(N, N);

  gsl_matrix * c  = gsl_matrix_alloc(N, N);

  gsl_matrix * minv = gsl_matrix_alloc(N, N);

  gsl_permutation * p = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(v, m);

  s += gsl_linalg_pcholesky_decomp(v, p);

  s += gsl_linalg_pcholesky_invert(v, p, minv);

  /* c = m m^{-1} */

  gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)

    {

      for (j = 0; j < N; ++j)

        {

          double cij = gsl_matrix_get(c, i, j);

          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",

                       desc, N, N, i, j, cij, expected);

        }

    }

  gsl_matrix_free(v);

  gsl_matrix_free(c);

  gsl_matrix_free(minv);

  gsl_permutation_free(p);

  return s;

}

int

test_pcholesky_invert(gsl_rng * r)

{

  int s = 0;

  const size_t N_max = 30;

  size_t N;

  for (N = 1; N <= N_max; ++N)

    {

      gsl_matrix * m = gsl_matrix_alloc(N, N);

      create_posdef_matrix(m, r);

      test_pcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "pcholesky_invert unscaled random");

      if (N <= 4)

        {

          create_hilbert_matrix2(m);

          test_pcholesky_invert_eps(m, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_invert unscaled hilbert");

        }

      gsl_matrix_free(m);

    }

  return s;

}