From ef1c20a77ebcdfb7fe11f480479e5c86b9d18abf Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 30 Jan 2020 13:18:00 +0100 Subject: pssolve: fail when the matrix condition number gets smaller than DBL_EPSILON --- pssolve.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pssolve.c b/pssolve.c index 32a34b1..b320b63 100644 --- a/pssolve.c +++ b/pssolve.c @@ -17,6 +17,7 @@ */ #include +#include #include #include #include @@ -156,6 +157,9 @@ static int lu_invert(TDLogger *logger, const int N, double *mat, double *rhs, in "condition number %16.16g; forward error %16.16g backward error %16.16g\n", N, N, cond, ferr, berr); + if (cond < DBL_EPSILON) + ret = -EDOM; + free(mat_f); free(x); #endif @@ -235,9 +239,13 @@ int tdi_pssolve_solve(PSSolveContext *ctx, memcpy(coeffs, rhs, s->nb_coeffs * sizeof(*rhs)); ret = lu_invert(&ctx->logger, s->nb_coeffs, s->mat, coeffs, s->ipiv); + ctx->lu_solves_time += gettime() - start; ctx->lu_solves_count++; + if (ret < 0) + return ret; + ret = tdi_bicgstab_init(s->bicgstab, s->mat, coeffs); s->steps_since_inverse = 0; -- cgit v1.2.3