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

1from __future__ import division 

2from __future__ import print_function 

3from .common import * 

4 

5def precondconjgrad(B, A, x, b, tolerance, maxiter, progress, relativeconv=False, robustresidual=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 z = B*r 

17 d = z 

18 rz = inner(r,z) 

19 if rz < 0: 19 ↛ 20line 19 didn't jump to line 20, because the condition on line 19 was never true

20 raise ValueError('Preconditioner not positive-definite') 

21 

22 iter = 0 

23 alphas = [] 

24 betas = [] 

25 residuals = [sqrt(rz)] 

26 

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

28 tolerance *= residuals[0] 

29 

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

31 z = A*d 

32 dz = inner(d,z) 

33 if dz == 0: 33 ↛ 34line 33 didn't jump to line 34, because the condition on line 33 was never true

34 print(f'ConjGrad breakdown') 

35 break 

36 alpha = rz/dz 

37 x += alpha*d 

38 if robustresidual: 38 ↛ 39line 38 didn't jump to line 39, because the condition on line 38 was never true

39 r = b - A*x 

40 else: 

41 r -= alpha*z 

42 z = B*r 

43 rz_prev = rz 

44 rz = inner(r,z) 

45 if rz < 0: 45 ↛ 46line 45 didn't jump to line 46, because the condition on line 45 was never true

46 print (f'Preconditioner not positive-definite') 

47 # Restore pre-breakdown state. Don't know if it helps any, but it's 

48 # consistent with returned quasi-residuals. 

49 x -= alpha*d 

50 break 

51 beta = rz/rz_prev 

52 d = z + beta*d 

53 

54 residual = sqrt(rz) 

55 

56 # Call user provided callback with solution 

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

58 callback(k=iter, x=x, r=residual) 

59 

60 iter += 1 

61 progress += 1 

62 alphas.append(alpha) 

63 betas.append(beta) 

64 residuals.append(residual) 

65 

66 return x, residuals, alphas, betas