diff options
author | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-07 18:24:48 +0000 |
---|---|---|
committer | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-07 18:24:48 +0000 |
commit | cca66b9ec4e494c1d919bff0f71a820d8afab1fa (patch) | |
tree | 146f39ded1c938019e1ed42d30923c2ac9e86789 /src/3rdparty/adaptagrams/libcola/tests/test_cg.cpp | |
parent | Initial commit. (diff) | |
download | inkscape-upstream.tar.xz inkscape-upstream.zip |
Adding upstream version 1.2.2.upstream/1.2.2upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to '')
-rw-r--r-- | src/3rdparty/adaptagrams/libcola/tests/test_cg.cpp | 164 |
1 files changed, 164 insertions, 0 deletions
diff --git a/src/3rdparty/adaptagrams/libcola/tests/test_cg.cpp b/src/3rdparty/adaptagrams/libcola/tests/test_cg.cpp new file mode 100644 index 0000000..13e7d02 --- /dev/null +++ b/src/3rdparty/adaptagrams/libcola/tests/test_cg.cpp @@ -0,0 +1,164 @@ +/* + * vim: ts=4 sw=4 et tw=0 wm=0 + * + * libcola - A library providing force-directed network layout using the + * stress-majorization method subject to separation constraints. + * + * Copyright (C) 2006-2008 Monash University + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library 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 + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library in the file LICENSE; if not, + * write to the Free Software Foundation, Inc., 59 Temple Place, + * Suite 330, Boston, MA 02111-1307 USA + * +*/ + +#include <stdio.h> +#include <cassert> +#include <libcola/conjugate_gradient.h> +#include <gsl/gsl_linalg.h> + +using std::valarray; + +static valarray<double> +outer_prod(valarray<double> x, valarray<double> y) { + valarray<double> result(x.size()*y.size()); + for(int j = 0; j < x.size(); j++) { + for(int i = 0; i < y.size(); i++) { + result[j*y.size() + i] = x[j]*y[i]; + } + } + return result; +} + +static double +uniform() { + return double(rand())/ RAND_MAX; +} + +static void +matrix_times_vector(valarray<double> const &matrix, /* m * n */ + valarray<double> const &vec, /* n */ + valarray<double> &result) /* m */ +{ + unsigned n = vec.size(); + unsigned m = result.size(); + assert(m*n == matrix.size()); + const double* mp = &matrix[0]; + for (unsigned i = 0; i < m; i++) { + double res = 0; + for (unsigned j = 0; j < n; j++) + res += *mp++ * vec[j]; + result[i] = res; + } +} + +static double +Linfty(valarray<double> const &vec) { + return std::max(vec.max(), -vec.min()); +} + +int +main (void) +{ + double tolerance = 1e-6; + for(int iters = 0; iters < 100; iters++) { + const unsigned N = unsigned(uniform()*40) + 1; + double A_data[N*N]; + printf("%ux%u matrix\n", N,N); + for(int r = 0; r < N; r++) { + for(int c = 0; c <= r; c++) { + A_data[r*N + c] = A_data[c*N + r] = fabs(uniform()); + } + } + + double * A_c[N]; + for(int i = 0; i < N; i++) + A_c[i] = &A_data[N*i]; + + double b_data[N]; + for(int i = 0; i < N; i++) + b_data[i] = uniform()*3; + std::valarray<double> b(b_data, N), xx(0.0, N); + std::valarray<double> A(A_data, N*N); + + conjugate_gradient(A, xx, b, N, tolerance, 2*N); + + gsl_matrix_view m + = gsl_matrix_view_array (A_data, N, N); + gsl_vector_view bgsl + = gsl_vector_view_array (b_data, N); + + gsl_vector *xgsl = gsl_vector_alloc (N); + + int s; + + gsl_permutation * p = gsl_permutation_alloc (N); + + gsl_linalg_LU_decomp (&m.matrix, p, &s); + + gsl_linalg_LU_solve (&m.matrix, p, &bgsl.vector, xgsl); + + std::valarray<double> gsl_xx(0.0, N); + for(unsigned i = 0; i < xx.size(); i++) { + gsl_xx[i] = gsl_vector_get(xgsl, i); + } + + double err = Linfty(xx - gsl_xx); + printf ("|xx-nxgsl|_infty = %g\n", err); + if(err > tolerance) { +#if 0 //dubious value + for(int r = 0; r < N; r++) { + for(int c = 0; c < N; c++) { + printf("%g ", A_data[r*N + c]); + } + printf("\n"); + } +#endif + printf("\nx njh-cg = \n"); + for(unsigned i = 0; i < xx.size(); i++) + printf("%g ", xx[i]); + printf("\nxgsl = "); + for(unsigned i = 0; i < xx.size(); i++) + printf("%g ", gsl_xx[i]); + printf("\n"); + valarray<double> result(0.0,N); + matrix_times_vector(A, xx, result); + result -= b; + printf("\nA xx -b = "); + for(unsigned i = 0; i < xx.size(); i++) + printf("%g ", result[i]); + printf("\n"); + matrix_times_vector(A, gsl_xx, result); + result -= b; + printf("\nA gsl_xx -b = "); + for(unsigned i = 0; i < xx.size(); i++) + printf("%g ", result[i]); + printf("\n"); + + printf("FAILED!!!!!!!!!!!!!!!!!!!!!!!!\n"); + exit(1); + } + } + return 0; +} +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 |