|
Packit |
67cb25 |
/* linalg/test_cholesky.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2016 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_test.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_ieee_utils.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_permute_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_permutation.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "test_common.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps,
|
|
Packit |
67cb25 |
const char * desc);
|
|
Packit |
67cb25 |
static int test_cholesky_decomp(gsl_rng * r);
|
|
Packit |
67cb25 |
int test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);
|
|
Packit |
67cb25 |
int test_cholesky_invert(gsl_rng * r);
|
|
Packit |
67cb25 |
static int test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps,
|
|
Packit |
67cb25 |
const char * desc);
|
|
Packit |
67cb25 |
static int test_pcholesky_decomp(gsl_rng * r);
|
|
Packit |
67cb25 |
int test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,
|
|
Packit |
67cb25 |
const gsl_vector * sol, const double eps,
|
|
Packit |
67cb25 |
const char * desc);
|
|
Packit |
67cb25 |
static int test_pcholesky_solve(gsl_rng * r);
|
|
Packit |
67cb25 |
int test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc);
|
|
Packit |
67cb25 |
int test_pcholesky_invert(gsl_rng * r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps, const char * desc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Hilbert matrix condition numbers, as calculated by LAPACK DPOSVX */
|
|
Packit |
67cb25 |
double hilb_rcond[] = { 1.000000000000e+00, 3.703703703704e-02, 1.336898395722e-03,
|
|
Packit |
67cb25 |
3.524229074890e-05, 1.059708198754e-06, 3.439939465186e-08,
|
|
Packit |
67cb25 |
1.015027593823e-09, 2.952221630602e-11, 9.093751565191e-13,
|
|
Packit |
67cb25 |
2.828277420229e-14, 8.110242564869e-16, 2.409320075800e-17 };
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_cholesky_decomp_eps(const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * V = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * L = gsl_matrix_calloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * LT = gsl_matrix_calloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * S = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(V, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
s += gsl_linalg_cholesky_decomp2(V, S);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
s += gsl_linalg_cholesky_decomp1(V);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L and LT */
|
|
Packit |
67cb25 |
gsl_matrix_tricpy('L', 1, L, V);
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 1, LT, L);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* L <- S^{-1} L, LT <- LT S^{-1} */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Si = gsl_vector_get(S, i);
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_row(L, i);
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_matrix_column(LT, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_scale(&v.vector, 1.0 / Si);
|
|
Packit |
67cb25 |
gsl_vector_scale(&w.vector, 1.0 / Si);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute A = L LT */
|
|
Packit |
67cb25 |
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(A, i, j);
|
|
Packit |
67cb25 |
double mij = gsl_matrix_get(m, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(Aij, mij, eps,
|
|
Packit |
67cb25 |
"%s: (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, Aij, mij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (expected_rcond > 0 && !scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector *work = gsl_vector_alloc(3 * N);
|
|
Packit |
67cb25 |
double rcond;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_rcond(V, &rcond, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(rcond, expected_rcond, 1.0e-6,
|
|
Packit |
67cb25 |
"%s rcond: (%3lu,%3lu): %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, rcond, expected_rcond);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(work);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(V);
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(L);
|
|
Packit |
67cb25 |
gsl_matrix_free(LT);
|
|
Packit |
67cb25 |
gsl_vector_free(S);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_cholesky_decomp(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
test_cholesky_decomp_eps(0, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp unscaled random");
|
|
Packit |
67cb25 |
test_cholesky_decomp_eps(1, m, -1.0, 1.0e2 * N * GSL_DBL_EPSILON, "cholesky_decomp scaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 12)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double expected_rcond = -1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hilb_rcond[N - 1] > 1.0e-12)
|
|
Packit |
67cb25 |
expected_rcond = hilb_rcond[N - 1];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_cholesky_decomp_eps(0, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp unscaled hilbert");
|
|
Packit |
67cb25 |
test_cholesky_decomp_eps(1, m, expected_rcond, N * GSL_DBL_EPSILON, "cholesky_decomp scaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_cholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * v = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * c = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(v, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s += gsl_linalg_cholesky_decomp1(v);
|
|
Packit |
67cb25 |
s += gsl_linalg_cholesky_invert(v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c = m m^{-1} */
|
|
Packit |
67cb25 |
gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, v, 0.0, c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c should be the identity matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cij = gsl_matrix_get(c, i, j);
|
|
Packit |
67cb25 |
double expected = (i == j) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, cij, expected);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(v);
|
|
Packit |
67cb25 |
gsl_matrix_free(c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_cholesky_invert(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_cholesky_invert_eps(m, N * GSL_DBL_EPSILON, "cholesky_invert unscaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 4)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
test_cholesky_invert_eps(m, N * 256.0 * GSL_DBL_EPSILON, "cholesky_invert unscaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_mcholesky_decomp_eps(const int posdef, const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps, const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * LDLT = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * V = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * L = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * LT = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * S = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * E = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_permutation * perm = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view D = gsl_matrix_diagonal(LDLT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(LDLT, m);
|
|
Packit |
67cb25 |
s += gsl_linalg_mcholesky_decomp(LDLT, perm, E);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that the upper triangle of LDLT equals original matrix */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double mij = gsl_matrix_get(m, i, j);
|
|
Packit |
67cb25 |
double aij = gsl_matrix_get(LDLT, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(aij, mij, 1.0e-12,
|
|
Packit |
67cb25 |
"%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, aij, mij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (posdef)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* ||E|| should be 0 */
|
|
Packit |
67cb25 |
double norm = gsl_blas_dnrm2(E);
|
|
Packit |
67cb25 |
s = norm != 0.0;
|
|
Packit |
67cb25 |
gsl_test(s, "%s: (%zu,%zu): ||E|| = %.12e",
|
|
Packit |
67cb25 |
desc, N, N, norm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that D is decreasing */
|
|
Packit |
67cb25 |
s = 0;
|
|
Packit |
67cb25 |
for (i = 1; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dprev = gsl_vector_get(&D.vector, i - 1);
|
|
Packit |
67cb25 |
double di = gsl_vector_get(&D.vector, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (di > dprev)
|
|
Packit |
67cb25 |
s = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test(s, "%s: (%zu,%zu): D is not decreasing",
|
|
Packit |
67cb25 |
desc, N, N);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L and LT */
|
|
Packit |
67cb25 |
gsl_matrix_set_identity(L);
|
|
Packit |
67cb25 |
gsl_matrix_set_identity(LT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_tricpy('L', 0, L, LDLT);
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, LT, L);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute (L sqrt(D)) and (sqrt(D) LT) */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_column(L, i);
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_matrix_row(LT, i);
|
|
Packit |
67cb25 |
double di = gsl_vector_get(&D.vector, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_scale(&v.vector, sqrt(di));
|
|
Packit |
67cb25 |
gsl_vector_scale(&w.vector, sqrt(di));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute A = L D LT */
|
|
Packit |
67cb25 |
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute V = P (S M S + E) P^T */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(V, m);
|
|
Packit |
67cb25 |
D = gsl_matrix_diagonal(V);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute S M S */
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale_apply(V, S);
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, V, V);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute S M S + E */
|
|
Packit |
67cb25 |
gsl_vector_add(&D.vector, E);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute M P^T */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_row(V, i);
|
|
Packit |
67cb25 |
gsl_permute_vector(perm, &v.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute P M P^T */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_column(V, i);
|
|
Packit |
67cb25 |
gsl_permute_vector(perm, &v.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Ei = gsl_vector_get(E, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(A, i, j); /* L D L^T */
|
|
Packit |
67cb25 |
double Bij = gsl_matrix_get(V, i, j); /* P M P^T */
|
|
Packit |
67cb25 |
double Cij; /* P M P^T + E */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == j)
|
|
Packit |
67cb25 |
Cij = Bij + Ei*0;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
Cij = Bij;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(Aij, Cij, eps,
|
|
Packit |
67cb25 |
"%s: (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, Aij, Cij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (expected_rcond > 0 && !scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector *work = gsl_vector_alloc(3 * N);
|
|
Packit |
67cb25 |
double rcond;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_mcholesky_rcond(LDLT, perm, &rcond, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(rcond, expected_rcond, 1.0e-6,
|
|
Packit |
67cb25 |
"%s rcond: (%3lu,%3lu): %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, rcond, expected_rcond);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(work);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(LDLT);
|
|
Packit |
67cb25 |
gsl_matrix_free(V);
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(L);
|
|
Packit |
67cb25 |
gsl_matrix_free(LT);
|
|
Packit |
67cb25 |
gsl_vector_free(S);
|
|
Packit |
67cb25 |
gsl_vector_free(E);
|
|
Packit |
67cb25 |
gsl_permutation_free(perm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_mcholesky_decomp(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
test_mcholesky_decomp_eps(1, 0, m, -1.0, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random posdef");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_symm_matrix(m, r);
|
|
Packit |
67cb25 |
test_mcholesky_decomp_eps(0, 0, m, -1.0, 8192.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled random symm");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 8)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double expected_rcond = -1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hilb_rcond[N - 1] > 1.0e-12)
|
|
Packit |
67cb25 |
expected_rcond = hilb_rcond[N - 1];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_mcholesky_decomp_eps(1, 0, m, expected_rcond, 128.0 * N * GSL_DBL_EPSILON, "mcholesky_decomp unscaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_mcholesky_solve_eps(const gsl_matrix * m, const gsl_vector * rhs,
|
|
Packit |
67cb25 |
const gsl_vector * sol, const double eps,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, N = m->size1;
|
|
Packit |
67cb25 |
gsl_matrix * u = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_calloc(N);
|
|
Packit |
67cb25 |
gsl_vector * S = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_permutation * perm = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(u, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s += gsl_linalg_mcholesky_decomp(u, perm, NULL);
|
|
Packit |
67cb25 |
s += gsl_linalg_mcholesky_solve(u, perm, rhs, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double yi = gsl_vector_get(sol, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(xi, yi, eps,
|
|
Packit |
67cb25 |
"%s: %3lu[%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, i, xi, yi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(S);
|
|
Packit |
67cb25 |
gsl_matrix_free(u);
|
|
Packit |
67cb25 |
gsl_permutation_free(perm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_mcholesky_solve(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * rhs = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * sol = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
create_random_vector(sol, r);
|
|
Packit |
67cb25 |
gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_mcholesky_solve_eps(m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "mcholesky_solve random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
|
|
Packit |
67cb25 |
test_mcholesky_solve_eps(m, rhs, sol, 1.0e3 * N * GSL_DBL_EPSILON, "mcholesky_solve hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
gsl_vector_free(rhs);
|
|
Packit |
67cb25 |
gsl_vector_free(sol);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_mcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * v = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * c = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * minv = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * E = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_permutation * p = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(v, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s += gsl_linalg_mcholesky_decomp(v, p, E);
|
|
Packit |
67cb25 |
s += gsl_linalg_mcholesky_invert(v, p, minv);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c = m m^{-1} */
|
|
Packit |
67cb25 |
gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c should be the identity matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cij = gsl_matrix_get(c, i, j);
|
|
Packit |
67cb25 |
double expected = (i == j) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, cij, expected);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(v);
|
|
Packit |
67cb25 |
gsl_matrix_free(c);
|
|
Packit |
67cb25 |
gsl_matrix_free(minv);
|
|
Packit |
67cb25 |
gsl_vector_free(E);
|
|
Packit |
67cb25 |
gsl_permutation_free(p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_mcholesky_invert(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 30;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
test_mcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "mcholesky_invert unscaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 4)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
test_mcholesky_invert_eps(m, 256.0 * N * GSL_DBL_EPSILON, "mcholesky_invert unscaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_pcholesky_decomp_eps(const int scale, const gsl_matrix * m,
|
|
Packit |
67cb25 |
const double expected_rcond, const double eps,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * LDLT = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * V = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * L = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * LT = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * S = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_permutation * perm = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view D = gsl_matrix_diagonal(LDLT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(LDLT, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_decomp2(LDLT, perm, S);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_decomp(LDLT, perm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that the upper triangle of LDLT equals original matrix */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double mij = gsl_matrix_get(m, i, j);
|
|
Packit |
67cb25 |
double aij = gsl_matrix_get(LDLT, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(aij, mij, 1.0e-12,
|
|
Packit |
67cb25 |
"%s upper triangle: (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, aij, mij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that D is decreasing */
|
|
Packit |
67cb25 |
s = 0;
|
|
Packit |
67cb25 |
for (i = 1; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dprev = gsl_vector_get(&D.vector, i - 1);
|
|
Packit |
67cb25 |
double di = gsl_vector_get(&D.vector, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (di > dprev)
|
|
Packit |
67cb25 |
s = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test(s, "%s: (%zu,%zu): D is not decreasing",
|
|
Packit |
67cb25 |
desc, N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L and LT */
|
|
Packit |
67cb25 |
gsl_matrix_set_identity(L);
|
|
Packit |
67cb25 |
gsl_matrix_set_identity(LT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_tricpy('L', 0, L, LDLT);
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, LT, L);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute (L sqrt(D)) and (sqrt(D) LT) */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_column(L, i);
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_matrix_row(LT, i);
|
|
Packit |
67cb25 |
double di = gsl_vector_get(&D.vector, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_scale(&v.vector, sqrt(di));
|
|
Packit |
67cb25 |
gsl_vector_scale(&w.vector, sqrt(di));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute A = L D LT */
|
|
Packit |
67cb25 |
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, L, LT, 0.0, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute V = P S M S P^T */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(V, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute S M S */
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale_apply(V, S);
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, V, V);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute M P^T */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_row(V, i);
|
|
Packit |
67cb25 |
gsl_permute_vector(perm, &v.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute P M P^T */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_column(V, i);
|
|
Packit |
67cb25 |
gsl_permute_vector(perm, &v.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(A, i, j); /* L D L^T */
|
|
Packit |
67cb25 |
double Bij = gsl_matrix_get(V, i, j); /* P M P^T */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(Aij, Bij, eps,
|
|
Packit |
67cb25 |
"%s: (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, Aij, Bij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (expected_rcond > 0 && !scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector *work = gsl_vector_alloc(3 * N);
|
|
Packit |
67cb25 |
double rcond;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_pcholesky_rcond(LDLT, perm, &rcond, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(rcond, expected_rcond, 1.0e-6,
|
|
Packit |
67cb25 |
"%s rcond: (%3lu,%3lu): %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, rcond, expected_rcond);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(work);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(LDLT);
|
|
Packit |
67cb25 |
gsl_matrix_free(V);
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(L);
|
|
Packit |
67cb25 |
gsl_matrix_free(LT);
|
|
Packit |
67cb25 |
gsl_vector_free(S);
|
|
Packit |
67cb25 |
gsl_permutation_free(perm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_pcholesky_decomp(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
test_pcholesky_decomp_eps(0, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled random");
|
|
Packit |
67cb25 |
test_pcholesky_decomp_eps(1, m, -1.0, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 12)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double expected_rcond = -1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hilb_rcond[N - 1] > 1.0e-12)
|
|
Packit |
67cb25 |
expected_rcond = hilb_rcond[N - 1];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_pcholesky_decomp_eps(0, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp unscaled hilbert");
|
|
Packit |
67cb25 |
test_pcholesky_decomp_eps(1, m, expected_rcond, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_decomp scaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_pcholesky_solve_eps(const int scale, const gsl_matrix * m, const gsl_vector * rhs,
|
|
Packit |
67cb25 |
const gsl_vector * sol, const double eps,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, N = m->size1;
|
|
Packit |
67cb25 |
gsl_matrix * u = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_calloc(N);
|
|
Packit |
67cb25 |
gsl_vector * S = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_permutation * perm = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(u, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (scale)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_decomp2(u, perm, S);
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_solve2(u, perm, S, rhs, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_decomp(u, perm);
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_solve(u, perm, rhs, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double yi = gsl_vector_get(sol, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(xi, yi, eps,
|
|
Packit |
67cb25 |
"%s: %3lu[%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, i, xi, yi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(S);
|
|
Packit |
67cb25 |
gsl_matrix_free(u);
|
|
Packit |
67cb25 |
gsl_permutation_free(perm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
test_pcholesky_solve(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 50;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * rhs = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * sol = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
create_random_vector(sol, r);
|
|
Packit |
67cb25 |
gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_pcholesky_solve_eps(0, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled random");
|
|
Packit |
67cb25 |
test_pcholesky_solve_eps(1, m, rhs, sol, 64.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
|
|
Packit |
67cb25 |
test_pcholesky_solve_eps(0, m, rhs, sol, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_solve unscaled hilbert");
|
|
Packit |
67cb25 |
test_pcholesky_solve_eps(1, m, rhs, sol, 2048.0 * N * GSL_DBL_EPSILON, "pcholesky_solve scaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
gsl_vector_free(rhs);
|
|
Packit |
67cb25 |
gsl_vector_free(sol);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_pcholesky_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
size_t i, j, N = m->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * v = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * c = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * minv = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_permutation * p = gsl_permutation_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(v, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_decomp(v, p);
|
|
Packit |
67cb25 |
s += gsl_linalg_pcholesky_invert(v, p, minv);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c = m m^{-1} */
|
|
Packit |
67cb25 |
gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, m, minv, 0.0, c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c should be the identity matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cij = gsl_matrix_get(c, i, j);
|
|
Packit |
67cb25 |
double expected = (i == j) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
|
|
Packit |
67cb25 |
desc, N, N, i, j, cij, expected);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(v);
|
|
Packit |
67cb25 |
gsl_matrix_free(c);
|
|
Packit |
67cb25 |
gsl_matrix_free(minv);
|
|
Packit |
67cb25 |
gsl_permutation_free(p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_pcholesky_invert(gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = 0;
|
|
Packit |
67cb25 |
const size_t N_max = 30;
|
|
Packit |
67cb25 |
size_t N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (N = 1; N <= N_max; ++N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
create_posdef_matrix(m, r);
|
|
Packit |
67cb25 |
test_pcholesky_invert_eps(m, N * GSL_DBL_EPSILON, "pcholesky_invert unscaled random");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N <= 4)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_hilbert_matrix2(m);
|
|
Packit |
67cb25 |
test_pcholesky_invert_eps(m, 1024.0 * N * GSL_DBL_EPSILON, "pcholesky_invert unscaled hilbert");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|