diff options
author | Anton Khirnov <anton@khirnov.net> | 2020-01-30 13:18:00 +0100 |
---|---|---|
committer | Anton Khirnov <anton@khirnov.net> | 2020-01-30 13:18:00 +0100 |
commit | ef1c20a77ebcdfb7fe11f480479e5c86b9d18abf (patch) | |
tree | 0ddab034953df054a750fd3736fe17e268d96d7a | |
parent | e14f7d3c807d2a46c8c10553f11556d94bb8998c (diff) |
pssolve: fail when the matrix condition number gets smaller than DBL_EPSILON
-rw-r--r-- | pssolve.c | 8 |
1 files changed, 8 insertions, 0 deletions
@@ -17,6 +17,7 @@ */ #include <errno.h> +#include <float.h> #include <inttypes.h> #include <limits.h> #include <math.h> @@ -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; |