Coverage for gpaw/tddft/solvers/bicgstab.py: 91%
82 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1# Written by Lauri Lehtovaara, 2008
2"""This module defines BiCGStab-class, which implements biconjugate
3gradient stabilized method. Requires Numpy and GPAW's own BLAS."""
5import numpy as np
7from gpaw.utilities.blas import axpy
8from gpaw.mpi import rank
10from .base import BaseSolver
13class BiCGStab(BaseSolver):
14 """Biconjugate gradient stabilized method
16 This class solves a set of linear equations A.x = b using biconjugate
17 gradient stabilized method (BiCGStab). The matrix A is a general,
18 non-singular matrix, e.g., it can be nonsymmetric, complex, and
19 indefinite. The method requires only access to matrix-vector product
20 A.x = b, which is called A.dot(x). Thus A must provide the member
21 function dot(self,x,b), where x and b are complex arrays
22 (numpy.array([], complex), and x is the known vector, and
23 b is the result.
25 Now x and b are multivectors, i.e., list of vectors.
26 """
28 def solve(self, A, x, b):
29 if self.timer is not None:
30 self.timer.start('BiCGStab')
32 # number of vectors
33 nvec = len(x)
35 # r_0 = b - A x_0
36 r = self.gd.zeros(nvec, dtype=complex)
37 A.dot(-x, r)
38 r += b
40 q = self.gd.zeros(nvec, dtype=complex)
41 q[:] = r
43 p = self.gd.zeros(nvec, dtype=complex)
44 v = self.gd.zeros(nvec, dtype=complex)
45 t = self.gd.zeros(nvec, dtype=complex)
46 m = self.gd.zeros(nvec, dtype=complex)
47 alpha = np.zeros((nvec, ), dtype=complex)
48 rho = np.zeros((nvec, ), dtype=complex)
49 rhop = np.zeros((nvec, ), dtype=complex)
50 omega = np.zeros((nvec, ), dtype=complex)
51 scale = np.zeros((nvec, ), dtype=complex)
52 tmp = np.zeros((nvec, ), dtype=complex)
54 rhop[:] = 1.
55 omega[:] = 1.
57 # Multivector dot product, a^H b, where ^H is conjugate transpose
58 def multi_zdotc(s, x, y, nvec):
59 for i in range(nvec):
60 s[i] = np.vdot(x[i], y[i])
61 self.gd.comm.sum(s)
62 return s
64 # Multivector ZAXPY: a x + y => y
65 def multi_zaxpy(a, x, y, nvec):
66 for i in range(nvec):
67 axpy(a[i] * (1 + 0J), x[i], y[i])
69 # Multiscale: a x => x
70 def multi_scale(a, x, nvec):
71 for i in range(nvec):
72 x[i] *= a[i]
74 # scale = square of the norm of b
75 multi_zdotc(scale, b, b, nvec)
76 scale = np.abs(scale)
78 # if scale < eps, then convergence check breaks down
79 if (scale < self.eps).any():
80 raise RuntimeError(
81 "BigCGStab method detected underflow for squared norm of"
82 " right-hand side (scale = %le < eps = %le)."
83 % (scale, self.eps))
85 # print 'Scale = ', scale
87 slow_convergence_iters = 50
89 for i in range(self.max_iter):
90 # rho_i-1 = q^H r_i-1
91 multi_zdotc(rho, q, r, nvec)
93 # print 'Rho = ', rho
95 # if i=1, p_i = r_i-1
96 # else beta = (rho_i-1 / rho_i-2) (alpha_i-1 / omega_i-1)
97 # p_i = r_i-1 + b_i-1 (p_i-1 - omega_i-1 v_i-1)
98 beta = (rho / rhop) * (alpha / omega)
100 # print 'Beta = ', beta
102 # if abs(beta) / scale < eps, then BiCGStab breaks down
103 if ((i > 0) and ((np.abs(beta) / scale) < self.eps).any()):
104 raise RuntimeError(
105 "Biconjugate gradient stabilized method failed"
106 " (abs(beta)=%le < eps = %le)."
107 % (np.min(np.abs(beta)), self.eps))
109 # p = r + beta * (p - omega * v)
110 # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
111 multi_zaxpy(-omega, v, p, nvec)
112 multi_scale(beta, p, nvec)
113 p += r
114 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
116 # v_i = A.(M^-1.p)
117 A.apply_preconditioner(p, m)
118 A.dot(m, v)
119 # alpha_i = rho_i-1 / (q^H v_i)
120 multi_zdotc(alpha, q, v, nvec)
121 alpha = rho / alpha
122 # s = r_i-1 - alpha_i v_i
123 multi_zaxpy(-alpha, v, r, nvec)
124 # s is denoted by r
126 # print 'Alpha = ', alpha
128 # x_i = x_i-1 + alpha_i (M^-1.p_i) + omega_i (M^-1.s)
129 # next line is x_i = x_i-1 + alpha (M^-1.p_i)
130 multi_zaxpy(alpha, m, x, nvec)
132 # if ( |s|^2 < tol^2 ) done
133 multi_zdotc(tmp, r, r, nvec)
134 if ((np.abs(tmp) / scale) < self.tol * self.tol).all():
135 # print 'R2 of proc #', rank, ' = ' , tmp, \
136 # ' after ', i+1, ' iterations'
137 break
139 # print if slow convergence
140 if ((i + 1) % slow_convergence_iters) == 0:
141 print('Log10 S2 of proc #', rank, ' = ',
142 np.round(np.log10(np.abs(tmp)), 1),
143 ' after ', i + 1, ' iterations')
145 # t = A.(M^-1.s), M = 1
146 A.apply_preconditioner(r, m)
147 A.dot(m, t)
148 # omega_i = t^H s / (t^H t)
149 multi_zdotc(omega, t, r, nvec)
150 multi_zdotc(tmp, t, t, nvec)
151 omega = omega / tmp
153 # print 'Omega = ', omega
155 # x_i = x_i-1 + alpha_i (M^-1.p_i) + omega_i (M^-1.s)
156 # next line is x_i = ... + omega_i (M^-1.s)
157 multi_zaxpy(omega, m, x, nvec)
158 # r_i = s - omega_i * t
159 multi_zaxpy(-omega, t, r, nvec)
160 # s is no longer denoted by r
162 # if ( |r|^2 < tol^2 ) done
163 multi_zdotc(tmp, r, r, nvec)
164 if ((np.abs(tmp) / scale) < self.tol * self.tol).all():
165 # print 'R2 of proc #', rank, ' = ' , tmp, \
166 # ' after ', i+1, ' iterations'
167 break
169 # print if slow convergence
170 if ((i + 1) % slow_convergence_iters) == 0:
171 print('Log10 R2 of proc #', rank, ' = ',
172 np.round(np.log10(np.abs(tmp)), 1),
173 ' after ', i + 1, ' iterations')
175 # if abs(omega) < eps, then BiCGStab breaks down
176 if ((np.abs(omega) / scale) < self.eps).any():
177 raise RuntimeError(
178 "Biconjugate gradient stabilized method failed"
179 " (abs(omega)/scale=%le < eps = %le)."
180 % (np.min(np.abs(omega)) / scale, self.eps))
181 # finally update rho
182 rhop[:] = rho
184 # if max iters reached, raise error
185 if (i >= self.max_iter - 1):
186 raise RuntimeError(
187 "Biconjugate gradient stabilized method failed to converge"
188 " within given number of iterations (= %d)."
189 % self.max_iter)
191 # done
192 self.iterations = i + 1
193 # print 'BiCGStab iterations = ', self.iterations
195 if self.timer is not None:
196 self.timer.stop('BiCGStab')
198 return self.iterations
199 # print self.iterations