Mercurial > octave

comparison scripts/sparse/bicg.m @ 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 883564a561aa
children 07c2c42f457e
comparison
equal deleted inserted replaced
25199:8cff7eceee87 25200:b3ee0179d7b0
104 ## 104 ##
105 ## @item 3: The algorithm stagnated, i.e., the absolute value of the 105 ## @item 3: The algorithm stagnated, i.e., the absolute value of the
106 ## difference between the current iteration @var{x} and the previous is less 106 ## difference between the current iteration @var{x} and the previous is less
107 ## than @code{eps * norm (@var{x},2)}. 107 ## than @code{eps * norm (@var{x},2)}.
108 ## 108 ##
109 ## @item 4: The algorithm can't continue due to a division by zero. 109 ## @item 4: The algorithm could not continue because intermediate values
110 ## became too small or too large for reliable computation.
110 ## @end itemize 111 ## @end itemize
111 ## 112 ##
112 ## @item 113 ## @item
113 ## @var{relres} is the ratio of the final residual to its initial value, 114 ## @var{relres} is the ratio of the final residual to its initial value,
114 ## measured in the Euclidean norm. 115 ## measured in the Euclidean norm.
269 end_try_catch 270 end_try_catch
270 271
271 while ((flag != 2) && (iter < maxit) && (resvec(iter+1) >= norm_b * tol)) 272 while ((flag != 2) && (iter < maxit) && (resvec(iter+1) >= norm_b * tol))
272 v = Afun (p, "notransp", varargin{:}); 273 v = Afun (p, "notransp", varargin{:});
273 prod_qv = q' * v; 274 prod_qv = q' * v;
274 if (prod_qv == 0) 275 alpha = (s0' * prec_r0);
276 if (abs (prod_qv) <= eps * abs (alpha))
275 flag = 4; 277 flag = 4;
276 break 278 break
277 endif 279 endif
278 alpha = (s0' * prec_r0) / prod_qv; 280 alpha ./= prod_qv;
279 x += alpha * p; 281 x += alpha * p;
280 prod_rs = (s0' * prec_r0); # Product between r0 and s0 282 prod_rs = (s0' * prec_r0); # Product between r0 and s0
281 r0 -= alpha * v; 283 r0 -= alpha * v;
282 s0 -= conj (alpha) * Afun (q, "transp", varargin{:}); 284 s0 -= conj (alpha) * Afun (q, "transp", varargin{:});
283 prec_r0 = M1fun (r0, "notransp", varargin{:}); 285 prec_r0 = M1fun (r0, "notransp", varargin{:});
284 prec_s0 = s0; 286 prec_s0 = s0;
285 prec_r0 = M2fun (prec_r0, "notransp", varargin{:}); 287 prec_r0 = M2fun (prec_r0, "notransp", varargin{:});
288 beta = s0' * prec_r0;
289 if (abs (prod_rs) <= abs (beta))
290 flag = 4;
291 break;
292 endif
293 beta ./= prod_rs;
286 prec_s0 = M2fun (prec_s0, "transp", varargin{:}); 294 prec_s0 = M2fun (prec_s0, "transp", varargin{:});
287 prec_s0 = M1fun (prec_s0, "transp", varargin{:}); 295 prec_s0 = M1fun (prec_s0, "transp", varargin{:});
288 iter += 1; 296 iter += 1;
289 resvec(iter+1) = norm (r0); 297 resvec(iter+1) = norm (r0);
290 if (resvec(iter+1) <= resvec(iter_min+1)) 298 if (resvec(iter+1) <= resvec(iter_min+1))
293 endif 301 endif
294 if (norm (x - x_pr) <= norm (x) * eps) 302 if (norm (x - x_pr) <= norm (x) * eps)
295 flag = 3; 303 flag = 3;
296 break; 304 break;
297 endif 305 endif
298 if (prod_rs == 0)
299 flag = 4;
300 break;
301 endif
302 beta = (s0' * prec_r0) / prod_rs;
303 p = prec_r0 + beta*p; 306 p = prec_r0 + beta*p;
304 q = prec_s0 + conj (beta) * q; 307 q = prec_s0 + conj (beta) * q;
305 endwhile 308 endwhile
306 resvec = resvec(1:iter+1,1); 309 resvec = resvec(1:iter+1,1);
307 310