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

1from __future__ import division 

2from __future__ import print_function 

3from builtins import range 

4"""Base class for iterative solvers.""" 

5 

6from block.block_base import block_base 

7import dolfin as df 

8import os 

9 

10class iterative(block_base): 

11 def __init__(self, A, precond=1.0, tolerance=1e-5, initial_guess=None, 

12 iter=None, maxiter=200, name=None, show=1, rprecond=None, 

13 nonconvergence_is_fatal=False, retain_guess=False, relativeconv=False, 

14 callback=None, **kwargs): 

15 

16 self.B = precond 

17 self.R = rprecond 

18 self.A = A 

19 self.initial_guess = initial_guess 

20 self.retain_guess = retain_guess 

21 self.nonconvergence_is_fatal = nonconvergence_is_fatal 

22 self.show = show 

23 self.callback = callback 

24 self.relativeconv = relativeconv 

25 self.kwargs = kwargs 

26 self.name = name if name else self.__class__.__name__ 

27 if iter is not None: 

28 tolerance = 0 

29 maxiter = iter 

30 self.tolerance = tolerance 

31 self.maxiter = maxiter 

32 self.cputime = 0. 

33 

34 def create_vec(self, dim=1): 

35 return self.A.create_vec(dim) 

36 

37 def matvec(self, b): 

38 from time import time 

39 from block.block_vec import block_vec 

40# from dolfin import log, info, Progress 

41 TRACE = 13 # dolfin.TRACE 

42 

43 T = time() 

44 

45 # If x and initial_guess are block_vecs, some of the blocks may be 

46 # scalars (although they are normally converted to vectors by bc 

47 # application). To be sure, call allocate() on them. 

48 

49 if isinstance(b, block_vec) and not b.allocated(): 49 ↛ 51line 49 didn't jump to line 51, because the condition on line 49 was never true

50 # Create a shallow copy to call allocate() on, to avoid changing the caller's copy of b 

51 b = block_vec(len(b), b.blocks) 

52 b.allocate(self.A, dim=0) 

53 

54 if self.initial_guess is not None: 

55 # Most (all?) solvers modify x, so make a copy to avoid changing 

56 # the caller's copy of x 

57 from block.block_util import copy 

58 x = copy(self.initial_guess) 

59 if isinstance(x, block_vec): 

60 x.allocate(self.A, dim=1) 

61 else: 

62 x = self.A.create_vec(dim=1) 

63 x.zero() 

64 

65 try: 

66# log(TRACE, self.__class__.__name__+' solve of '+str(self.A)) 

67# if self.B != 1.0: 

68# log(TRACE, 'Using preconditioner: '+str(self.B)) 

69# progress = Progress(self.name, self.maxiter) 

70 if self.tolerance < 0: 70 ↛ 71line 70 didn't jump to line 71, because the condition on line 70 was never true

71 tolerance = -self.tolerance 

72 relative = True 

73 else: 

74 tolerance = self.tolerance 

75 relative = False 

76 progress = 0 

77 x = self.method(self.B, self.AR, x, b, tolerance=tolerance, 

78 relativeconv=self.relativeconv, maxiter=self.maxiter, 

79 progress=progress, callback=self.callback, 

80 **self.kwargs) 

81 del progress # trigger final printout 

82 except Exception: 

83# from dolfin import warning 

84 print("Error solving " + self.name) 

85 raise 

86 x, self.residuals, self.alphas, self.betas = x 

87 

88 if self.tolerance == 0: 

89 msg = "done" 

90 elif self.converged: 90 ↛ 93line 90 didn't jump to line 93, because the condition on line 90 was never false

91 msg = "converged" 

92 else: 

93 msg = "NOT CONV." 

94 

95 cputime = time()-T 

96 

97 comm = df.MPI.comm_world 

98 # Get the average 

99 self.cputime = comm.allreduce(cputime)/comm.size 

100 if self.show == 1: 

101 print('%s %s [iter=%2d, time=%.2fs, res=%.1e]' \ 

102 % (self.name, msg, self.iterations, self.cputime, self.residuals[-1])) 

103 elif self.show >= 2: 

104 print('%s %s [iter=%2d, time=%.2fs, res=%.1e, true res=%.1e]' \ 

105 % (self.name, msg, self.iterations, self.cputime, self.residuals[-1], (self.A*x-b).norm('l2'))) 

106 if self.show == 3 and not "DOLFIN_NOPLOT" in os.environ: 106 ↛ 107line 106 didn't jump to line 107, because the condition on line 106 was never true

107 if comm.rank == 0: 

108 try: 

109 from matplotlib import pyplot 

110 pyplot.figure('%s convergence (show=3)'%self.name) 

111 pyplot.semilogy(self.residuals) 

112 pyplot.show(block=True) 

113 except: 

114 pass 

115 

116 if self.R is not None: 116 ↛ 117line 116 didn't jump to line 117, because the condition on line 116 was never true

117 x = self.R*x 

118 

119 if self.retain_guess: 119 ↛ 120line 119 didn't jump to line 120, because the condition on line 119 was never true

120 self.initial_guess = x 

121 

122 if True: 

123 # Tag vectors with iteration counts so that they can be picked up 

124 # by regression tests. 

125 x._regr_test_niter = self.iterations 

126 if isinstance(x, block_base): 

127 for block in x: 

128 block._regr_test_niter = self.iterations 

129 

130 if not self.converged and self.nonconvergence_is_fatal: 130 ↛ 131line 130 didn't jump to line 131, because the condition on line 130 was never true

131 raise RuntimeError('Not converged') 

132 

133 return x 

134 

135 def __call__(self, initial_guess=None, precond=None, tolerance=None, 

136 iter=None, maxiter=None, name=None, show=None, rprecond=None, 

137 nonconvergence_ok=None, callback=None): 

138 """Allow changing the parameters within an expression, e.g. x = Ainv(initial_guess=x) * b""" 

139 import copy 

140 ret = copy.copy(self) # shallow copy 

141 if precond is not None: ret.B = precond 

142 if rprecond is not None: ret.R = rprecond 

143 if initial_guess is not None: ret.initial_guess = initial_guess 

144 if nonconvergence_ok is not None: ret.nonconvergence_ok = nonconvergence_ok 

145 if show is not None: ret.show = show 

146 if name is not None: ret.name = name 

147 if tolerance is not None: ret.tolerance = tolerance 

148 if maxiter is not None: ret.maxiter = maxiter 

149 if callback is not None: ret.callback = callback 

150 if iter is not None: 

151 ret.tolerance = 0 

152 ret.maxiter = iter 

153 return ret 

154 

155 @property 

156 def AR(self): 

157 return self.A if self.R is None else self.A*self.R 

158 

159 @property 

160 def iterations(self): 

161 return len(self.residuals)-1 

162 @property 

163 def converged(self): 

164 eff_tolerance = self.tolerance 

165 if eff_tolerance == 0: 

166 return True 

167 if eff_tolerance < 0: 167 ↛ 168line 167 didn't jump to line 168, because the condition on line 167 was never true

168 eff_tolerance *= -self.residuals[0] 

169 if self.relativeconv: 

170 eff_tolerance *= self.residuals[0] 

171 return self.residuals[-1] <= eff_tolerance 

172 

173 def eigenvalue_estimates(self): 

174 ##### 

175 # Adapted from code supplied by KAM (Simula PyCC; GPL license), 

176 ##### 

177 

178 # eigenvalues estimates in terms of alphas and betas 

179 

180 import numpy 

181 

182 n = len(self.alphas) 

183 if n == 0: 183 ↛ 184line 183 didn't jump to line 184, because the condition on line 183 was never true

184 raise RuntimeError('Eigenvalues can not be estimated, no alphas/betas') 

185 M = numpy.zeros([n,n]) 

186 M[0,0] = 1/self.alphas[0] 

187 for k in range(1, n): 

188 M[k,k] = 1/self.alphas[k] + self.betas[k-1]/self.alphas[k-1] 

189 M[k,k-1] = numpy.sqrt(self.betas[k-1])/self.alphas[k-1] 

190 M[k-1,k] = M[k,k-1] 

191 e,v = numpy.linalg.eig(M) 

192 e.sort() 

193 return e 

194 

195 def __str__(self): 

196 return '<%d %s iterations on %s>'%(self.maxiter, self.name, self.A) 

197