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

1from __future__ import division 

2from .common import * 

3 

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

5 

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 # Is this correct?? Should we not transpose the preconditoner somewhere?? 

16 

17 # jobh 02/2011: Changed residual from sqrt(rho) to sqrt(inner(r,r)) due to negative rho 

18 

19 r = b - A*x 

20 Br = B*r 

21 ATBr = transpmult(A, Br) 

22 BATBr = B*ATBr 

23 

24 

25 rho = inner(BATBr,ATBr) 

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

27 raise RuntimeError('CGN: Preconditioner not positive-definite') 

28 rho1 = rho 

29 p = BATBr.copy() 

30 

31 iter = 0 

32 alphas = [] 

33 betas = [] 

34 residuals = [sqrt(rho)] 

35 

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

37 tolerance *= residuals[0] 

38 

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

40 Ap = A*p 

41 BAp = B*Ap; 

42 ATBAp = transpmult(A, BAp) 

43 alpha = rho/inner(p,ATBAp) 

44 x = x + alpha*p 

45 r = b-A*x 

46 Br = B*r 

47 ATBr = transpmult(A, Br) 

48 BATBr = B*ATBr 

49 

50 rho = inner(BATBr,ATBr) 

51 if rho < 0: 51 ↛ 52line 51 didn't jump to line 52, because the condition on line 51 was never true

52 raise RuntimeError('CGN: Preconditioner not positive-definite') 

53 beta = rho/rho1 

54 rho1 = rho 

55 p = BATBr+beta*p 

56 

57 residual = sqrt(rho) 

58 

59 # Call user provided callback with solution 

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

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

62 

63 iter += 1 

64 progress += 1 

65 alphas.append(alpha) 

66 betas.append(beta) 

67 residuals.append(residual) 

68 

69 return x, residuals, alphas, betas