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

1from __future__ import print_function 

2from builtins import range 

3from .block_mat import block_mat 

4from .block_compose import block_mul 

5 

6def block_kronecker(A, B): 

7 """Create the Kronecker (tensor) product of two matrices. The result is 

8 returned as a product of two block matrices, (A x Ib) * (Ia x B), where Ia 

9 and Ib are the appropriate identities (A=A*Ia, B=Ib*B). This will often 

10 limit the number of repeated applications of the inner operators. 

11 

12 Note: If the A operator is a DOLFIN matrix and B is not a block_mat, A will 

13 be expanded into a block_mat (one block per matrix entry). This will not 

14 work in parallel. Apart from that, no blocks are expanded, i.e., only the 

15 block matrices are expanded in the Kronecker product. 

16 

17 To form the Kronecker sum, you can extract (A x Ib) and (Ia x B) like this: 

18 C,D = block_kronecker(A,B); sum=C+D 

19 Similarly, it may be wise to do the inverse separately: 

20 C,D = block_kronecker(A,B); inverse = some_invert(D)*ConjGrad(C) 

21 """ 

22 from .block_util import isscalar 

23 import dolfin 

24 

25 if isinstance(A, dolfin.Matrix) and not isinstance(B, block_mat): 

26 A = block_mat(A.array()) 

27 assert isinstance(A, block_mat) or isinstance(B, block_mat) 

28 

29 if isinstance(B, block_mat): 

30 p,q = B.blocks.shape 

31 C = block_mat.diag(A, n=p) 

32 else: 

33 C = A 

34 

35 if isinstance(A, block_mat): 

36 m,n = A.blocks.shape 

37 if isinstance(B, block_mat): 

38 print ("Note: block_kronecker with two block matrices is probably buggy") 

39 D = block_mat(n,n) 

40 for i in range(n): 

41 for j in range(n): 

42 # A scalar can represent the scaled identity of any 

43 # dimension, so no need to diagonal-expand it here. We 

44 # don't do this check on the outer diagonal expansions, 

45 # because it is clearer to always return two block matrices 

46 # of equal dimension rather than sometimes a scalar. 

47 b = B[i,j] 

48 D[i,j] = b if isscalar(b) else block_mat.diag(b, n=m) 

49 else: 

50 D = block_mat.diag(B, n=m) 

51 else: 

52 D = B 

53 

54 return block_mul(C,D) 

55 

56def block_simplify(expr): 

57 """Return a simplified (if possible) representation of a block matrix or 

58 block composition. The simplification does the following basic steps: 

59 - Convert scaled identity matrices to scalars 

60 - Combine scalar terms in compositions (2*2==4) 

61 - Eliminate additive and multiplicative identities (A+0=A, A*1=A) 

62 """ 

63 if hasattr(expr, 'block_simplify'): 

64 return expr.block_simplify() 

65 else: 

66 return expr 

67 

68 

69def block_collapse(expr): 

70 """Turn a composition /inside out/, i.e., turn a composition of block 

71 matrices into a block matrix of compositions. 

72 """ 

73 if hasattr(expr, 'block_collapse'): 

74 return expr.block_collapse() 

75 else: 

76 return expr