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