Mercurial > octave

changeset 25200:b3ee0179d7b0 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
bicg.m: Check for denominators much smaller than numerators bicg.m (bug #53589). * bicg.m: Set flag = 4 if the denominators are much smaller than the numerators, not just if they are equal to zero.
author Marco Caliari <marco.caliari@univr.it>
date Tue, 10 Apr 2018 08:39:24 +0200
parents 8cff7eceee87
children c80323fe4938
files scripts/sparse/bicg.m
diffstat 1 files changed, 11 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/sparse/bicg.m	Mon Apr 09 18:35:45 2018 -0500
+++ b/scripts/sparse/bicg.m	Tue Apr 10 08:39:24 2018 +0200
@@ -106,7 +106,8 @@
 ## difference between the current iteration @var{x} and the previous is less
 ## than @code{eps * norm (@var{x},2)}.
 ##
-## @item 4: The algorithm can't continue due to a division by zero.
+## @item 4: The algorithm could not continue because intermediate values
+## became too small or too large for reliable computation.
 ## @end itemize
 ##
 ## @item
@@ -271,11 +272,12 @@
   while ((flag != 2) && (iter < maxit) && (resvec(iter+1) >= norm_b * tol))
     v = Afun (p, "notransp", varargin{:});
     prod_qv = q' * v;
-    if (prod_qv == 0)
+    alpha = (s0' * prec_r0);
+    if (abs (prod_qv) <= eps * abs (alpha))
       flag = 4;
       break
     endif
-    alpha = (s0' * prec_r0) / prod_qv;
+    alpha ./= prod_qv;
     x += alpha * p;
     prod_rs = (s0' * prec_r0);  # Product between r0 and s0
     r0 -= alpha * v;
@@ -283,6 +285,12 @@
     prec_r0 = M1fun (r0, "notransp", varargin{:});
     prec_s0 = s0;
     prec_r0 = M2fun (prec_r0, "notransp", varargin{:});
+    beta = s0' * prec_r0;
+    if (abs (prod_rs) <= abs (beta))
+      flag = 4;
+      break;
+    endif
+    beta ./= prod_rs;
     prec_s0 = M2fun (prec_s0, "transp", varargin{:});
     prec_s0 = M1fun (prec_s0, "transp", varargin{:});
     iter += 1;
@@ -295,11 +303,6 @@
       flag = 3;
       break;
     endif
-    if (prod_rs == 0)
-      flag = 4;
-      break;
-    endif
-    beta = (s0' * prec_r0) / prod_rs;
     p = prec_r0 + beta*p;
     q = prec_s0 + conj (beta) * q;
   endwhile