A wise man once began, I have this matrix A
... (known here as "W")
use LinearAlgebra,
LayoutCS,
LinearAlgebra.Sparse;
var nv: int = 8,
D = {1..nv, 1..nv},
SD: sparse subdomain(D) dmapped CS(),
W: [SD] real;
SD += (1,2); W[1,2] = 1.0;
SD += (1,3); W[1,3] = 1.0;
SD += (1,4); W[1,4] = 1.0;
SD += (2,2); W[2,2] = 1.0;
SD += (2,4); W[2,4] = 1.0;
SD += (3,4); W[3,4] = 1.0;
SD += (4,5); W[4,5] = 1.0;
SD += (5,6); W[5,6] = 1.0;
SD += (6,7); W[6,7] = 1.0;
SD += (6,8); W[6,8] = 1.0;
SD += (7,8); W[7,8] = 1.0;
const a: real = 0.5;
I would like to build the polynomial
const P = aW + a^2W^2 + .. + a^kW^k
However, it appears that the function .dot
can't be chained. Is there a clear way to do this would building the intermediate elements?
Here's one way to achieve this, but it could be improved by splitting the polynomial computation up such that each matrix power is only computed once:
use LinearAlgebra,
LayoutCS,
LinearAlgebra.Sparse;
var nv: int = 8,
D = {1..nv, 1..nv},
SD: sparse subdomain(D) dmapped CS(),
W: [SD] real;
SD += (1,2); W[1,2] = 1.0;
SD += (1,3); W[1,3] = 1.0;
SD += (1,4); W[1,4] = 1.0;
SD += (2,2); W[2,2] = 1.0;
SD += (2,4); W[2,4] = 1.0;
SD += (3,4); W[3,4] = 1.0;
SD += (4,5); W[4,5] = 1.0;
SD += (5,6); W[5,6] = 1.0;
SD += (6,7); W[6,7] = 1.0;
SD += (6,8); W[6,8] = 1.0;
SD += (7,8); W[7,8] = 1.0;
const a: real = 0.5;
const polynomial = dot(a, W).plus(dot(a**2, W.dot(W))).plus(dot(a**3, W.dot(W).dot(W)));