Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/iterative/minres2.py: 8%
64 statements
« prev ^ index » next coverage.py v7.2.1, created at 2023-03-20 13:03 +0000
2"""
3@license: Copyright (C) 2006
4Author: Gunnar Staff
6Simula Research Laboratory AS
8This file is part of PyCC.
10PyCC free software; you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation; either version 2 of the License, or
13(at your option) any later version.
15PyCC is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with PyFDM; if not, write to the Free Software
22Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25MinRes: A generic implementation of the Minimum Residual method.
26"""
29MinResError = "Error in MinRes"
31import numpy as _n #/numarray
32#from Utils import _n,inner
33from math import sqrt,fabs
35def inner(u,v):
36 """Compute innerproduct of u and v.
37 It is not computed here, we just check type of operands and
38 call a suitable innerproduct method"""
40 if isinstance(u,_n.ndarray):
41 # assume both are numarrays:
42 return _n.dot(u,v)
43 else:
44 # assume that left operand implements inner
45 return u.inner(v)
48def minres(B, A, x, b, tolerance, maxiter, progress, relativeconv=False, rit_=1, callback=None):
49 """
50 precondMinRes(B,A,x,b): Solve Ax = b with the preconditioned minimum
51 residual method.
53 @param B: Preconditioner supporting the __mul__ operator for operating on
54 x. B must be symmetric positive definite
56 @param A: Matrix or discrete operator. Must support A*x to operate on x.
57 A must be symmetric, but may be indefinite.
59 @param x: Anything that A can operate on, as well as support inner product
60 and multiplication/addition/subtraction with scalar and something of the
61 same type
63 @param b: Right-hand side, probably the same type as x.
65 @param tolerance: Convergence criterion
66 @type tolerance: float
68 @param relativeconv: Control relative vs. absolute convergence criterion.
69 That is, if relativeconv is True, ||r||_2/||r_init||_2 is used as
70 convergence criterion, else just ||r||_2 is used. Actually ||r||_2^2 is
71 used, since that save a few cycles.
73 @type relativeconv: bool
75 @return: the solution x.
77 DESCRIPTION:
79 The method implements a left preconditioned Minimum residual method
80 for symmetric indefinite linear systems. The preconditioning
81 operator has to be symmetric and positive definite.
83 The minimum residual method solves systems of the form BAx=Bb,
84 where A is symmetric indefinite and B is symmetric positive
85 definite. The linear system is symmetric with respect to the inner
86 product $ \((\cdot,\cdot)_B = (B^{-1}\cdot,\cdot)\). The iterate
87 x_k is determined by minimization of the norm of the residual $ \[
88 \|B(b - A y)\|_B \] over the Krylov space $ \(span\{Bb, BABb,
89 \ldots, (BA)^{k-1}Bb\}\). $ Here the norm is defined by the inner
90 product \((\cdot,\cdot)_B\).
92 The default convergence monitor is $ \[ \rho_k = \|B(b - A x_k)\|_B
93 = (B(b - A x_k), b - A x_k).\] $ The residual \(b - A x_k\) is not
94 computed during the iteration, hence a direct computation of this
95 quantity reqire an additional matrix vector product. In the
96 algorithm it is computed recursively. Unfortunately this
97 computations accumulates error and it may be necessary to compute
98 the exact residual every update frequency iteration.
102 """
104 callback_converged = False
106 residuals = []
108 rit = rit_
109 d = A*x
110 r = b - d
111 s = B*r
112 rho = inner(r,s)
115 po = s.copy()
116 qo = A*po
118 p = r.copy()
119 p *= 0.0
120 q = p.copy()
121 u = p.copy()
123 if relativeconv:
124 tolerance *= sqrt(rho)
126 residuals.append(sqrt(rho))
127 iter = 0
128 #print "tolerance ", tolerance
129 # Alloc w
130 #while sqrt(inner(r,r)) > tolerance:# and iter<maxiter:
131 while (sqrt(rho) > tolerance and not callback_converged) and iter < maxiter:
132 #print "iter, sqrt(inner(r,r))", iter, sqrt(inner(r,r))
133 uo = B*qo
134 gamma = sqrt(inner(qo,uo))
135 gammai= 1.0/gamma
137 # Do a swap
138 tmp = po
139 po = p
140 p = tmp
141 p *= gammai
143 tmp = qo
144 qo = q
145 q = tmp
146 q *= gammai
148 tmp = uo
149 uo = u
150 u = tmp
151 u *= gammai
153 alpha = inner(s,q)
154 x += alpha*p
155 s -= alpha*u
157 if iter%rit==0:
158 r = b - A*x
159 rho = inner(r,s)
160 else:
161 rho -= alpha*alpha
162 rho = fabs(rho)
164 t = A*u
165 alpha = inner(t,u)
166 beta = inner(t,uo)
168 po *= -beta
169 po -= alpha*p
170 po += u
172 # update qo
173 #po-=beta*po
174 qo *= -beta
175 qo -= alpha*q
176 qo += t
179 residuals.append(sqrt(rho))
180 #print "sqrt(rho) ", sqrt(rho)
182 # Call user provided callback with solution
183 if callable(callback):
184 callback_converged = callback(k=iter, x=x, r=sqrt(rho)) #r=r)
186 iter += 1
187 #print "r",iter,"=",sqrt(inner(r,r))
188 #print "precondMinRes finished, iter: %d, ||e||=%e" % (iter,sqrt(inner(r,r)))
189 return x,residuals, [], []