Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/iterative/bicgstab.py: 73%

1from __future__ import division 

2from __future__ import print_function 

3from .common import * 

4 

5def precondBiCGStab(B, A, x, b, tolerance, maxiter, progress, relativeconv=False, callback=None): 

6 ##### 

7 # Adapted from code supplied by KAM (Simula PyCC; GPL license). This code 

8 # relicensed under LGPL v2.1 or later, in agreement with the original 

9 # authors: 

10 # Kent Andre Mardal <kent-and@simula.no> 

11 # Ola Skavhaug <skavhaug@simula.no> 

12 # Gunnar Staff <gunnaran@simula.no> 

13 ##### 

14 

15 r = b - A*x 

16 

17 p = r.copy() 

18 r0 = r.copy() 

19 rr0 = inner(r,r0) 

20 

21 iter = 0 

22 alphas = [] 

23 betas = [] 

24 residuals = [sqrt(rr0)] 

25 

26 if relativeconv: 26 ↛ 27line 26 didn't jump to line 27, because the condition on line 26 was never true

27 tolerance *= residuals[0] 

28 

29 while residuals[-1] > tolerance and iter < maxiter: 

30 Bp = B*p 

31 ABp = A*Bp 

32 alpha = rr0/inner(r0,ABp) 

33 s = r-alpha*ABp 

34 Bs = B*s 

35 ABs = A*Bs 

36 ABsABs = inner(ABs,ABs) 

37 sABs = inner(ABs,s) 

38 if ABsABs == 0.0 or sABs == 0.0: 38 ↛ 39line 38 didn't jump to line 39, because the condition on line 38 was never true

39 print ("BiCGStab breakdown (zero inner product)") 

40 return x, residuals, alphas, betas 

41 w = sABs/ABsABs 

42 x += alpha*Bp+w*Bs 

43 r = s - w*ABs 

44 rrn = inner(r,r0) 

45 

46 residual = sqrt(inner(r,r)) 

47 

48 if residual == 0.0: 48 ↛ 49line 48 didn't jump to line 49, because the condition on line 48 was never true

49 print ("BiCGStab breakdown (zero residual)") 

50 return x, residuals, alphas, betas 

51 

52 # Call user provided callback with solution 

53 if callable(callback): 53 ↛ 54line 53 didn't jump to line 54, because the condition on line 53 was never true

54 newres = callback(k=iter, x=x, r=residual) 

55 if newres is not None: 

56 residual = newres 

57 

58 beta = (rrn/rr0)*(alpha/w) 

59 if beta==0.0: 59 ↛ 60line 59 didn't jump to line 60, because the condition on line 59 was never true

60 print(("BiCGStab breakdown, beta=0, at iter=",iter," with residual=", residual)) 

61 return x, residuals, alphas, betas 

62 rr0 = rrn 

63 p = r+beta*(p-w*ABp) 

64 

65 iter += 1 

66 progress += 1 

67 alphas.append(alpha) 

68 betas.append(beta) 

69 residuals.append(residual) 

70 

71 return x, residuals, alphas, betas