mm410-- source --
-- Matrix Multiplication by PDC, Mou 1991 Algo.4.10
-- Consider A as a stack of two submatrices, B as a side-by-side of two
-- then C = [ [ A0 B0, A0 B1 ], [ A1 B0, A1 B1 ] ] (2x2 blocks)
--
-- Definitions:
--
N = 4
d = outer @ ( d_lr(0), d_lr(1))
c = 'c_lr^2'
reduce = PDC(d_lr,id,id,[+],atomq,id)
mm_b = wrap @ reduce @ [*] @ (first,first) @ (id,transpose)
baseq = AND @ (isrow,iscol)
mm410 = PDC(d,c,id,id,baseq,mm_b)
--
-- Application
--
A = randMat(N,N,fourth_roots(),seed=42)
B = randMat(N,N,fourth_roots(),seed=43)
trace(mm410, A, B)
-- end source --
-- trace: mm410([[ 1, 1, -1, 1.0j] [[ 1, -1, 1.0j, -1.0j]
[ 1.0j, 1.0j, 1, 1] [ -1, 1, -1.0j, -1.0j]
[-1.0j, 1, 1, 1] [ 1, -1.0j, -1, -1.0j]
[ 1.0j, 1.0j, 1, 1.0j]] [-1.0j, 1.0j, 1, 1]]) --
mm410([[ 1, 1, -1, 1.0j] [[ 1, -1, 1.0j, -1.0j]
[ 1.0j, 1.0j, 1, 1] [ -1, 1, -1.0j, -1.0j]
[-1.0j, 1, 1, 1] [ 1, -1.0j, -1, -1.0j]
[ 1.0j, 1.0j, 1, 1.0j]] [-1.0j, 1.0j, 1, 1]])
divide d -> (
( [[ 1, 1, -1, 1.0j] [[ 1, -1]
[1.0j, 1.0j, 1, 1]] [ -1, 1]
[ 1, -1.0j]
[-1.0j, 1.0j]] ),
( [[ 1, 1, -1, 1.0j] [[ 1.0j, -1.0j]
[1.0j, 1.0j, 1, 1]] [-1.0j, -1.0j]
[ -1, -1.0j]
[ 1, 1]] ),
( [[-1.0j, 1, 1, 1] [[ 1, -1]
[ 1.0j, 1.0j, 1, 1.0j]] [ -1, 1]
[ 1, -1.0j]
[-1.0j, 1.0j]] ),
( [[-1.0j, 1, 1, 1] [[ 1.0j, -1.0j]
[ 1.0j, 1.0j, 1, 1.0j]] [-1.0j, -1.0j]
[ -1, -1.0j]
[ 1, 1]] )
)
mm410([[ 1, 1, -1, 1.0j] [[ 1, -1]
[1.0j, 1.0j, 1, 1]] [ -1, 1]
[ 1, -1.0j]
[-1.0j, 1.0j]])
divide d -> (
( [[1, 1, -1, 1.0j]] [[ 1]
[ -1]
[ 1]
[-1.0j]] ),
( [[1, 1, -1, 1.0j]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]] ),
( [[1.0j, 1.0j, 1, 1]] [[ 1]
[ -1]
[ 1]
[-1.0j]] ),
( [[1.0j, 1.0j, 1, 1]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]] )
)
mm410([[1, 1, -1, 1.0j]] [[ 1]
[ -1]
[ 1]
[-1.0j]])
⇣ atom; basef -> [[0.0]]
mm410([[1, 1, -1, 1.0j]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]])
⇣ atom; basef -> [[-1.0+1.0j]]
mm410([[1.0j, 1.0j, 1, 1]] [[ 1]
[ -1]
[ 1]
[-1.0j]])
⇣ atom; basef -> [[1.0-1.0j]]
mm410([[1.0j, 1.0j, 1, 1]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]])
⇣ atom; basef -> [[0.0]]
combine c_lr_sq -> [[ 0.0, -1.0+1.0j]
[1.0-1.0j, 0.0]]
mm410([[ 1, 1, -1, 1.0j] [[ 1.0j, -1.0j]
[1.0j, 1.0j, 1, 1]] [-1.0j, -1.0j]
[ -1, -1.0j]
[ 1, 1]])
divide d -> (
( [[1, 1, -1, 1.0j]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]] ),
( [[1, 1, -1, 1.0j]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]] ),
( [[1.0j, 1.0j, 1, 1]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]] ),
( [[1.0j, 1.0j, 1, 1]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]] )
)
mm410([[1, 1, -1, 1.0j]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]])
⇣ atom; basef -> [[1.0+1.0j]]
mm410([[1, 1, -1, 1.0j]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]])
⇣ atom; basef -> [[0.0]]
mm410([[1.0j, 1.0j, 1, 1]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]])
⇣ atom; basef -> [[0.0]]
mm410([[1.0j, 1.0j, 1, 1]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]])
⇣ atom; basef -> [[3.0-1.0j]]
combine c_lr_sq -> [[1.0+1.0j, 0.0]
[ 0.0, 3.0-1.0j]]
mm410([[-1.0j, 1, 1, 1] [[ 1, -1]
[ 1.0j, 1.0j, 1, 1.0j]] [ -1, 1]
[ 1, -1.0j]
[-1.0j, 1.0j]])
divide d -> (
( [[-1.0j, 1, 1, 1]] [[ 1]
[ -1]
[ 1]
[-1.0j]] ),
( [[-1.0j, 1, 1, 1]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]] ),
( [[1.0j, 1.0j, 1, 1.0j]] [[ 1]
[ -1]
[ 1]
[-1.0j]] ),
( [[1.0j, 1.0j, 1, 1.0j]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]] )
)
mm410([[-1.0j, 1, 1, 1]] [[ 1]
[ -1]
[ 1]
[-1.0j]])
⇣ atom; basef -> [[-2.0j]]
mm410([[-1.0j, 1, 1, 1]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]])
⇣ atom; basef -> [[1.0+1.0j]]
mm410([[1.0j, 1.0j, 1, 1.0j]] [[ 1]
[ -1]
[ 1]
[-1.0j]])
⇣ atom; basef -> [[2.0]]
mm410([[1.0j, 1.0j, 1, 1.0j]] [[ -1]
[ 1]
[-1.0j]
[ 1.0j]])
⇣ atom; basef -> [[-1.0-1.0j]]
combine c_lr_sq -> [[-2.0j, 1.0+1.0j]
[ 2.0, -1.0-1.0j]]
mm410([[-1.0j, 1, 1, 1] [[ 1.0j, -1.0j]
[ 1.0j, 1.0j, 1, 1.0j]] [-1.0j, -1.0j]
[ -1, -1.0j]
[ 1, 1]])
divide d -> (
( [[-1.0j, 1, 1, 1]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]] ),
( [[-1.0j, 1, 1, 1]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]] ),
( [[1.0j, 1.0j, 1, 1.0j]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]] ),
( [[1.0j, 1.0j, 1, 1.0j]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]] )
)
mm410([[-1.0j, 1, 1, 1]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]])
⇣ atom; basef -> [[1.0-1.0j]]
mm410([[-1.0j, 1, 1, 1]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]])
⇣ atom; basef -> [[-2.0j]]
mm410([[1.0j, 1.0j, 1, 1.0j]] [[ 1.0j]
[-1.0j]
[ -1]
[ 1]])
⇣ atom; basef -> [[-1.0+1.0j]]
mm410([[1.0j, 1.0j, 1, 1.0j]] [[-1.0j]
[-1.0j]
[-1.0j]
[ 1]])
⇣ atom; basef -> [[2.0]]
combine c_lr_sq -> [[ 1.0-1.0j, -2.0j]
[-1.0+1.0j, 2.0]]
combine c_lr_sq -> [[ 0.0, -1.0+1.0j, 1.0+1.0j, 0.0]
[1.0-1.0j, 0.0, 0.0, 3.0-1.0j]
[ -2.0j, 1.0+1.0j, 1.0-1.0j, -2.0j]
[ 2.0, -1.0-1.0j, -1.0+1.0j, 2.0]]
-- result: [[ 0.0, -1.0+1.0j, 1.0+1.0j, 0.0]
[1.0-1.0j, 0.0, 0.0, 3.0-1.0j]
[ -2.0j, 1.0+1.0j, 1.0-1.0j, -2.0j]
[ 2.0, -1.0-1.0j, -1.0+1.0j, 2.0]]
--
|