← Back to branch summary

~cbc.block/cbc.block/dev

« back to all changes in this revision

Viewing changes to block/block_transform.py

  • Committer: Joachim Berdal Haga
  • Date: 2013-01-01 20:07:01 UTC
  • Revision ID: jobh@simula.no-20130101200701-eq8ngmpjz07smu4u
Clean up Kronecker

Show diffs side-by-side

added added

removed removed

Lines of Context:
block/block_transform.py
3
3
 
4
4
def block_kronecker(A, B):
5
5
    """Create the Kronecker (tensor) product of two matrices. The result is
6
 
    returned as a product of two block matrices, (A x Ib) * (Ia x B), because
7
 
    this will often limit the number of repeated applications of the inner
8
 
    operators. However, it also means that A and B must be square since
9
 
    otherwise the identities Ia and Ib are not defined.
 
6
    returned as a product of two block matrices, (A x Ib) * (Ia x B), where Ia
 
7
    and Ib are the appropriate identities (A=A*Ia, B=Ib*B). This will often
 
8
    limit the number of repeated applications of the inner operators.
10
9
 
11
 
    Note: If none of the parameters are of type block_mat, the composition A*B
12
 
    is returned instead of two block matrices.
 
10
    Note: If the A operator is a DOLFIN matrix and B is not a block_mat, A will
 
11
    be expanded into a block_mat (one block per matrix entry). This will not
 
12
    work in parallel. Apart from that, no blocks are expanded, i.e., only the
 
13
    block matrices are expanded in the Kronecker product.
13
14
 
14
15
    To form the Kronecker sum, you can extract (A x Ib) and (Ia x B) like this:
15
16
      C,D = block_kronecker(A,B); sum=C+D
17
18
      C,D = block_kronecker(A,B); inverse = some_invert(D)*ConjGrad(C)
18
19
    """
19
20
    from block_util import isscalar
 
21
    import dolfin
 
22
 
 
23
    if isinstance(A, dolfin.GenericMatrix) and not isinstance(B, block_mat):
 
24
        A = block_mat(A.array())
 
25
    assert isinstance(A, block_mat) or isinstance(B, block_mat)
20
26
 
21
27
    if isinstance(B, block_mat):
22
 
        n = len(B.blocks)
23
 
        C = block_mat.diag(A, n=n)
 
28
        p,q = B.blocks.shape
 
29
        C = block_mat.diag(A, n=p)
24
30
    else:
25
31
        C = A
26
32
 
27
33
    if isinstance(A, block_mat):
28
 
        m = len(A.blocks)
 
34
        m,n = A.blocks.shape
29
35
        if isinstance(B, block_mat):
30
 
            D = block_mat(m,m)
31
 
            for i in range(m):
32
 
                for j in range(m):
 
36
            print "Note: block_kronecker with two block matrices is probably buggy"
 
37
            D = block_mat(n,n)
 
38
            for i in range(n):
 
39
                for j in range(n):
33
40
                    # A scalar can represent the scaled identity of any
34
41
                    # dimension, so no need to diagonal-expand it here. We
35
42
                    # don't do this check on the outer diagonal expansions,
44
51
 
45
52
    return block_mul(C,D)
46
53
 
47
 
 
48
54
def block_simplify(expr):
49
55
    """Return a simplified (if possible) representation of a block matrix or
50
56
    block composition. The simplification does the following basic steps: