Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/block_scheme.py: 44%

1from __future__ import division 

2 

3from builtins import range 

4"""These classes are typically not used directly, but returned by a call to 

5block_mat.scheme(). 

6""" 

7 

8from .block_base import block_base 

9 

10def block_jacobi(op): 

11 from .block_mat import block_mat 

12 m,n = op.blocks.shape 

13 assert m==n 

14 mat = block_mat(m,n) 

15 for i in range(m): 

16 mat[i,i] = op[i,i] 

17 return mat 

18 

19class block_gs(block_base): 

20 def __init__(self, op, reverse=False, truncated=False, symmetric=False, w=1.0): 

21 self.op = op 

22 self.range = list(range(len(op))) 

23 if reverse: 23 ↛ 24line 23 didn't jump to line 24, because the condition on line 23 was never true

24 self.range.reverse() 

25 if symmetric: 25 ↛ 27line 25 didn't jump to line 27, because the condition on line 25 was never false

26 self.matvec = self.matvec_symmetric 

27 elif truncated: 

28 self.matvec = self.matvec_truncated 

29 else: 

30 self.matvec = self.matvec_full 

31 self.weight = w 

32 

33 def sor_weighting(self, b, x): 

34 if self.weight != 1: 34 ↛ 35line 34 didn't jump to line 35, because the condition on line 34 was never true

35 for i in self.range: 

36 x[i] *= self.weight 

37 x[i] += (1-self.weight) * b[i] 

38 

39 def matvec_full(self, b): 

40 x = b.copy() 

41 for i in self.range: 

42 for j in self.range: 

43 if j==i: 

44 continue 

45 if self.op[i,j]: 

46 x[i] -= self.op[i,j]*x[j] 

47 x[i] = self.op[i,i]*x[i] 

48 self.sor_weighting(b, x) 

49 return x 

50 

51 def matvec_truncated(self, b): 

52 x = b.copy() 

53 for k,i in enumerate(self.range): 

54 for j in self.range[:k]: 

55 if self.op[i,j]: 

56 x[i] -= self.op[i,j]*x[j] 

57 x[i] = self.op[i,i]*x[i] 

58 self.sor_weighting(b, x) 

59 return x 

60 

61 def matvec_symmetric(self, b): 

62 x = b.copy() 

63 for k,i in enumerate(self.range): 

64 for j in self.range[:k]: 

65 if self.op[i,j]: 65 ↛ 64line 65 didn't jump to line 64, because the condition on line 65 was never false

66 x[i] -= self.op[i,j]*x[j] 

67 x[i] = self.op[i,i]*x[i] 

68 rev_range = list(reversed(self.range)) 

69 for k,i in enumerate(rev_range): 

70 for j in rev_range[:k]: 

71 if self.op[i,j]: 71 ↛ 70line 71 didn't jump to line 70, because the condition on line 71 was never false

72 x[i] -= self.op[i,i]*self.op[i,j]*x[j] 

73 self.sor_weighting(b, x) 

74 return x 

75 

76def blockscheme(op, scheme='jacobi', **kwargs): 

77 scheme = scheme.lower() 

78 if scheme == 'jacobi' or scheme == 'jac': 78 ↛ 79line 78 didn't jump to line 79, because the condition on line 78 was never true

79 return block_jacobi(op) 

80 

81 if scheme in ('gauss-seidel', 'gs', 'sor'): 81 ↛ 82line 81 didn't jump to line 82, because the condition on line 81 was never true

82 return block_gs(op, **kwargs) 

83 

84 if scheme in ('symmetric gauss-seidel', 'sgs', 'ssor'): 84 ↛ 87line 84 didn't jump to line 87, because the condition on line 84 was never false

85 return block_gs(op, symmetric=True, **kwargs) 

86 

87 if scheme in ('truncated gauss-seidel', 'tgs', 'tsor'): 

88 return block_gs(op, truncated=True, **kwargs) 

89 

90 raise TypeError('unknown scheme "%s"'%scheme)