aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2020-01-30 13:18:00 +0100
committerAnton Khirnov <anton@khirnov.net>2020-01-30 13:18:00 +0100
commitef1c20a77ebcdfb7fe11f480479e5c86b9d18abf (patch)
tree0ddab034953df054a750fd3736fe17e268d96d7a
parente14f7d3c807d2a46c8c10553f11556d94bb8998c (diff)
pssolve: fail when the matrix condition number gets smaller than DBL_EPSILON
-rw-r--r--pssolve.c8
1 files changed, 8 insertions, 0 deletions
diff --git a/pssolve.c b/pssolve.c
index 32a34b1..b320b63 100644
--- a/pssolve.c
+++ b/pssolve.c
@@ -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;