Σ-SPL
Gather and Scatter
Gather and Scatter Matrices
g1 := Gath(fId(2)); # Gath(fId(.)) = I(.)
g2 := Gath(fBase(4, 0)); # Gath(fBase(.,.)) = base vec
g3 := Gath(fTensor(fBase(4, 0), fId(2))); # standard pattern
s1 := Scat(fId(2)); # Scat(fId(.)) = I(.)
s2 := Scat(fBase(4, 2)); # Scat(fBase(.,.)) = base vec
s3 := Scat(fTensor(fId(2), fBase(4, 3))); # standard pattern
Scatter/Kernel/Gather Pattern
A := F(2); j := 0;
# iteration j of Tensor(I(4), F(2))
sag1 := Scat(fTensor(fBase(4, j), fId(2))) * A *
Gath(fTensor(fBase(4, j), fId(2)));
# iteration j of Tensor(F(2), I(4))
sag2 := Scat(fTensor(fId(2), fBase(4, j))) * A *
Gath(fTensor(fId(2), fBase(4, j)));
# iteration j of Tensor(I(4), F(2)) * L(8, 4)
sag3 := Scat(fTensor(fBase(4, j), fId(2))) * A *
Gath(fTensor(fId(2), fBase(4, j)));
pm(sag1); pm(sag2); pm(sag3);
Tensor Product and Gath/Scat/ISum
A := F(2);
j := Ind(4);
# Sigma-SPL for Tensor(I(4), F(2))
sag1 := Scat(fTensor(fBase(j), fId(2))) * A *
Gath(fTensor(fBase(j), fId(2)));
s1 := ISum(j, sag1);
MatSPL(s1) = MatSPL(Tensor(I(4), F(2)));
Other Tensor Patterns
# Sigma-SPL for Tensor(F(2), I(4))
s2 := ISum(j, Scat(fTensor(fId(2), fBase(j))) * A *
Gath(fTensor(fId(2), fBase(j))));
MatSPL(s2) = MatSPL(Tensor(F(2), I(4)));
# Sigma-SPL for Tensor(I(4), F(2)) * L(8, 4)
s3 := ISum(j, Scat(fTensor(fBase(j), fId(2))) * A *
Gath(fTensor(fId(2), fBase(j))));
MatSPL(s3) = MatSPL(Tensor(I(4), F(2)) * L(8, 4));
# Direct sum
ISum(j, Scat(fTensor(fBase(j), fId(2))) * Mat([[j,-j],[j, j]]) *
Gath(fTensor(fBase(j), fId(2))));
Advanced Loops
Tensor Product and Gath/Scat/ISum
A := F(2);
j := Ind(4);
k := Ind(2);
# Sigma-SPL for Tensor(I(4), F(2), I(2))
sag := Scat(fTensor(fBase(j), fId(2), fBase(k))) * A *
Gath(fTensor(fBase(j), fId(2), fBase(k)));
s := ISum(k, ISum(j, sag));
MatSPL(s) = MatSPL(Tensor(I(4), F(2), I(2)));
More Complex Example
i := Ind(8);
j := Ind(4);
k := Ind(2);
A := Mat([[j+1,-2*j],[j*k, j+1]]); B := F(2);
s := ISum(k, ISum(j, Scat(fTensor(fBase(j), fBase(k), fId(2))) * A *
Gath(fTensor(fId(2), fBase(k), fBase(j))))) *
ISum(i, Scat(fTensor(fBase(i), fId(2))) * B
* Gath(fTensor(J(2), fCompose(Z(8,3), fBase(i)))));
MatSPL(s);
Gath/Scat, Diag and Perms
Gather Functions
# use Lambda functions directly in Gath/Scat
i := Ind(8);
f1 := Lambda(i, imod(5*i+7, 32)).setRange(32);
g := Gath(f1);
# Use indirection tables in Gath/Scat. By default not supported
f2 := CopyFields(FData(List([0..7], j->V(Mod(5*j+7, 32)))),
rec(range := self >> 32));
s := Scat(f2);
Diagonals
# the true Twiddle diagonal in DFT(8)
d := Diag(fPrecompute(fCompose(dOmega(8, 1),
diagTensor(dLin(V(4), 1, 0, TInt), dLin(2, 1, 0, TInt)))));
MatSPL(d);
e:= RCDiag(fCompose(FData(List([1..32], u->Value(TReal, u))),
fTensor(fId(4), fBase(i))));
# print out the function values
e.element.tolist();
# access the indirection table
e.element.children()[1].var.value;
Defining Σ-SPL Objects
Gather
# spiral-core\namespaces\spiral\spl\sums.gi
Class(Gath, SumsBase, BaseMat, rec(
rChildren := self >> [self.func],
rSetChild := rSetChildFields("func"),
new := (self, func) >> SPL(WithBases(self, rec(func :=
Checked(IsFunction(func) or IsFuncExp(func), func)))).setDims(),
dims := self >> [self.func.domain(), self.func.range()],
sums := self >> self,
area := self >> Sum(Flat([self.func.domain()])),
isReal := self >> true,
transpose := self >> Scat(self.func),
conjTranspose := self >> self.transpose(),
inverse := self >> self.transpose(),
toAMat := self >> let(
n := EvalScalar(self.func.domain()),
N := EvalScalar(self.func.range()),
func := self.func.lambda(),
AMatMat(List([0..n-1], row -> let(idx:=evalScalar(func.at(row)),
Cond(idx _is funcExp, When(idx.args[1]=0, Replicate(N, 0),
Error("... ")),
BasisVec(N, idx)))))),
));
Σ-SPL Index Mapping Functions
Definition
f := fId(4); # I4->I4: i->i
j := Ind(4); # variable with range
g := fBase(j); # I1->I4: i->j
h := fAdd(4,2,1); # I2->I4: i->i+1
u := L(16, 4); # permutation i -> \ell(16,4)(i)
Operations on Functions
f.at(0); # evaluate function
f.tolist(); # create table for function
f.lambda(); # convert to Lambda function
r := fTensor(f, g); # tensor product of functions
s := fCompose(r, u); # function composition
Function Properties/Fields
f.domain();
f.range();
Σ-SPL Diagonal Functions
Definition
f := fConst(4, 1.1); # I4->R: i->1.1
g := dOmega(8, 2); # f_N,k : N -> C : i -> omega(N, k*i)
h := FList(TReal, [1.1, 1.2, 1.4, 1.4]); # table lookup
u := FData([V(1.1), V(1.2), V(1.3), V(1.4)]); # table lookup w/var
Operations on Functions
g.at(3); # evaluate function
g.at(3).ev(); # simplify the result
f.tolist(); # create table for function
g.lambda(); # convert to Lambda function
r := diagTensor(f, u); # tensor product of functions
r.tolist(); # what does diagTensor do?
s := fCompose(r, L(16,4)); # composition of permutation and
s.tolist(); # tensor product of functions
Function Properties/Fields
f.domain();
f.range();
u.var;
u.var.t;
u.var.value;
Defining Σ-SPL Symbolic Functions
fId and fBase
# spiral-core\namespaces\spiral\spl\perms2.gi
Class(fId, PermClass, rec(
domain := self >> self.params[1],
range := self >> self.params[1],
def := size -> Checked(IsPosInt0Sym(size), rec(size := size)),
lambda := self >> let(i := Ind(self.params[1]), Lambda(i,i)),
transpose := self >> self,
isIdentity := True
));
Class(fBase, FuncClass, rec(
abbrevs := [ var -> Checked(IsVar(var) or
ObjId(var)=ind, [var.range, var]) ],
def := (N, pos) -> rec(),
domain := self >> 1,
range := self >> self.params[1],
print := (self, i, is) >> Print(self.name, "(",
When(ObjId(self.params[2]) in [var, ind],
Print(self.params[2]), PrintCS(self.params)), ")"),
lambda := self >> let(i := Ind(1), Lambda(i, self.params[2]))
));
fCompose
# spiral-core\namespaces\spiral\spl\perms2.gi
Class(fCompose, FuncClassOper, rec(
domain := self >> Last(self._children).domain(),
range := self >> self._children[1].range(),
lambda := self >>
FoldL1(List(self.children(), z->z.lambda()),
_rankedLambdaCompose),
transpose := self >> self.__bases__[1](
List(Reversed(self.children()), c->c.transpose())),
isIdentity := self >> ForAll(self._children, IsIdentity),
));
dOmega and dLin for Twiddle Diagonal
Class(dLin, DiagFunc, rec(
checkParams := (self, params) >> Checked(Length(params)=4,
IsPosInt0Sym(params[1]), IsType(params[4]), params),
lambda := self >> let(i:=Ind(self.params[1]),
a := self.params[2], b := self.params[3],
Lambda(i, a*i+b).setRange(self.params[4])),
range := self >> self.params[4],
domain := self >> self.params[1]
));
Class(dOmega, DiagFunc, rec(
checkParams := (self, params) >> Checked(Length(params)=2,
IsPosIntSym(params[1]), IsIntSym(params[2]), params),
lambda := self >> let(i:=Ind(), Lambda(i,
omega(self.params[1], self.params[2]*i)).setRange(TComplex)),
range := self >> TComplex,
domain := self >> TInt
));