CODEWITHSUNDEEP
sparse-matrixchapel
Brian Dolan
Brian Dolan

Reputation: 3136

How to build a matrix polynomial in Chapel

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?

Upvotes: 1

Views: 29

Answers (1)

ben-albrecht
ben-albrecht

Reputation: 1865

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)));

Upvotes: 1

Related Questions