Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/block_mat.py: 60%
130 statements
« prev ^ index » next coverage.py v7.2.1, created at 2023-03-20 13:03 +0000
1from __future__ import division
2from builtins import range
3from .block_base import block_container
4from .block_vec import block_vec
6class block_mat(block_container):
7 """Block of matrices or other operators. Empty blocks doesn't need to be set
8 (they may be None or zero), but each row must have at least one non-empty block.
10 As a special case, the diagonal may contain scalars which act as the
11 (scaled) identity operators."""
13 def __init__(self, m, n=None, blocks=0):
14 if n is None:
15 blocks = m
16 m = len(blocks)
17 n = len(blocks[0]) if m else 0
18 block_container.__init__(self, (m,n), blocks)
20 def matvec(self, x):
21 from dolfin import GenericVector, Matrix, PETScMatrix
22 from .block_util import isscalar
23 import numpy
24 m,n = self.blocks.shape
25 y = block_vec(m)
27 for i in range(m):
28 for j in range(n):
29 if isinstance(self[i,j], (numpy.ndarray, numpy.matrix)): 29 ↛ 30line 29 didn't jump to line 30, because the condition on line 29 was never true
30 z = numpy.matrix(self[i,j]) * numpy.matrix(x[j].get_local()).transpose()
31 z = numpy.array(z).flatten()
32 if y[i] is None:
33 y[i] = x[j].copy()
34 y[i].resize(len(z))
35 y[i][:] += z[:]
36 continue
37 if self[i,j] is None or self[i,j]==0 \ 37 ↛ 40line 37 didn't jump to line 40, because the condition on line 37 was never true
38 or x[j] is None or (isscalar(x[j]) and x[j]==0):
39 # Skip multiply if zero
40 continue
41 if self[i,j] == 1: 41 ↛ 43line 41 didn't jump to line 43, because the condition on line 41 was never true
42 # Skip multiply if identity
43 z = x[j]
44 else:
45 # Do the block multiply
46 if isinstance(self[i,j], PETScMatrix): 46 ↛ 47line 46 didn't jump to line 47, because the condition on line 46 was never true
47 z = self[i,j].create_vec(dim=0)
48 self[i,j].mult(x[j], z)
49 else:
50 z = self[i,j] * x[j]
51 if isinstance(z, type(NotImplemented)): return NotImplemented 51 ↛ exitline 51 didn't return from function 'matvec', because the return on line 51 wasn't executed
52 if not isinstance(z, (GenericVector, block_vec)): 52 ↛ 61line 52 didn't jump to line 61, because the condition on line 52 was never true
53 # Typically, this happens when for example a
54 # block_vec([0,0]) is used without calling allocate() or
55 # setting BCs. The result is a Matrix*scalar=Matrix. One
56 # way to fix this issue is to convert all scalars to a
57 # proxy class in block_vec.__init__, and let this proxy
58 # class have a __rmul__ that converts to vector on demand.
59 # (must also stop conversion anything->blockcompose for
60 # this case)
61 raise RuntimeError(
62 'unexpected result in matvec, %s\n-- possibly because a block_vec contains scalars ' \
63 'instead of vectors, use create_vec() or allocate()' % type(z))
64 if y[i] is None:
65 y[i] = z
66 else:
67 if len(y[i]) != len(z): 67 ↛ 68line 67 didn't jump to line 68, because the condition on line 67 was never true
68 raise RuntimeError(
69 'incompatible dimensions in block (%d,%d) -- %d, was %d'%(i,j,len(z),len(y[i])))
70 y[i] += z
71 return y
73 def transpmult(self, x):
74 import numpy
75 from dolfin import GenericVector
77 m,n = self.blocks.shape
78 y = block_vec(self.blocks.shape[0])
80 for i in range(n):
81 for j in range(m):
82 if self[j,i] is None or self[j,i]==0:
83 # Skip multiply if zero
84 continue
85 if self[i,j] == 1: 85 ↛ 87line 85 didn't jump to line 87, because the condition on line 85 was never true
86 # Skip multiply if identity
87 z = x[j]
88 elif numpy.isscalar(self[j,i]):
89 # mult==transpmult
90 z = self[j,i]*x[j]
91 else:
92 # Do the block multiply
93 z = self[j,i].transpmult(x[j])
94 if not isinstance(z, (GenericVector, block_vec)): 94 ↛ 96line 94 didn't jump to line 96, because the condition on line 94 was never true
95 # see comment in matvec
96 raise RuntimeError(
97 'unexpected result in matvec, %s\n-- possibly because RHS contains scalars ' \
98 'instead of vectors, use create_vec() or allocate()' % type(z))
99 if y[i] is None:
100 y[i] = z
101 else:
102 if len(y[i]) != len(z): 102 ↛ 103line 102 didn't jump to line 103, because the condition on line 102 was never true
103 raise RuntimeError(
104 'incompatible dimensions in block (%d,%d) -- %d, was %d'%(i,j,len(z),len(y[i])))
105 y[i] += z
106 return y
108 def copy(self):
109 from . import block_util
110 m,n = self.blocks.shape
111 y = block_mat(m,n)
112 for i in range(m):
113 for j in range(n):
114 y[i,j] = block_util.copy(self[i,j])
115 return y
117 def create_vec(self, dim=1):
118 """Create a vector of the correct size and layout for b (if dim==0) or
119 x (if dim==1)."""
120 xx = block_vec(self.blocks.shape[dim])
121 xx.allocate(self, dim)
122 return xx
124 def scheme(self, name, inverse=None, **kwargs):
125 """Return a block scheme (block relaxation method). The input block_mat
126 is normally expected to be defined with inverses (or preconditioners)
127 on the diagonal, and the normal blocks off-diagonal. For example, given
128 a coefficient matrix
130 AA = block_mat([[A, B],
131 [C, D]]),
133 the Gauss-Seidel block preconditioner may be defined as for example
135 AApre = block_mat([[ML(A), B],
136 [C, ML(D)]]).scheme('gs')
138 However, for the common case of using the same preconditioner for each
139 block, the block preconditioner may be specified in a simpler way as
141 AApre = AA.scheme('gs', inverse=ML)
143 where "inverse" defines a class or method that may be applied to each
144 diagonal block to form its (approximate) inverse.
146 The returned block operator may be a block_mat itself (this is the case
147 for Jacobi), or a procedural operator (this is the case for the
148 Gauss-Seidel variants).
150 The scheme may be one of the following (aliases in brackets):
152 jacobi [jac]
153 gauss-seidel [gs]
154 symmetric gauss-seidel [sgs]
155 truncated gauss-seidel [tgs]
156 sor
157 ssor
158 tsor
160 and they may take keyword arguments. For the Gauss-Seidel-like methods,
161 'reverse' and 'w' (weight) are supported, as well as 'truncated' and
162 'symmetric'.
163 """
164 from .block_scheme import blockscheme
165 if inverse is not None: 165 ↛ 166line 165 didn't jump to line 166, because the condition on line 165 was never true
166 m,n = self.blocks.shape
167 mat = block_mat(m,n)
168 mat[:] = self[:]
169 for i in range(m):
170 mat[i,i] = inverse(mat[i,i])
171 else:
172 mat = self
173 return blockscheme(mat, name, **kwargs)
175 @staticmethod
176 def diag(A, n=0):
177 """Create a diagonal block matrix, where the entries on the diagonal
178 are either the entries of the block_vec or list A (if n==0), or n
179 copies of A (if n>0). For the case of extracting the diagonal of an
180 existing block matrix, use D=A.scheme('jacobi') instead.
182 If n>0, A is replicated n times and may be either an operator or a
183 block_vec/list. If it is a list, the entries of A define a banded
184 diagonal. In this case, len(A) must be odd, with the diagonal in the
185 middle.
186 """
187 if n==0:
188 n = len(A)
189 mat = block_mat(n,n)
190 for i in range(n):
191 mat[i,i] = A[i]
192 else:
193 if isinstance(A, (list, tuple, block_vec)):
194 if len(A)%2 != 1:
195 raise ValueError('The number of entries in the banded diagonal must be odd')
196 else:
197 A = [A]
198 mat = block_mat(n,n)
199 m = len(A)
200 for i in range(n):
201 for k in range(m):
202 j = i-m//2+k
203 if not 0 <= j < n:
204 continue
205 mat[i,j] = A[k]
206 return mat
208 def block_simplify(self):
209 """Try to convert identities to scalars, recursively. A fuller
210 explanation is found in block_transform.block_simplify.
211 """
212 from .block_util import isscalar
213 from .block_transform import block_simplify
214 m,n = self.blocks.shape
215 res = block_mat(m,n)
216 # Recursive call
217 for i in range(m):
218 for j in range(n):
219 res[i,j] = block_simplify(self[i,j])
220 # Check if result after recursive conversion is the (scaled) identity
221 v0 = res.blocks[0,0]
222 if m != n: 222 ↛ 223line 222 didn't jump to line 223, because the condition on line 222 was never true
223 return res
224 for i in range(m): 224 ↛ 235line 224 didn't jump to line 235, because the loop on line 224 didn't complete
225 for j in range(n): 225 ↛ 224line 225 didn't jump to line 224, because the loop on line 225 didn't complete
226 block = res.blocks[i,j]
227 if not isscalar(block): 227 ↛ 229line 227 didn't jump to line 229, because the condition on line 227 was never false
228 return res
229 if i==j:
230 if block!=v0:
231 return res
232 else:
233 if block!=0:
234 return res
235 return v0