text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# Don't forget to Import(platforms.sse)

#NOTE:
# Bluestein breaks for SSE_2x64f -- disabled for small for now

# doParSimdDft seems broken
benchSSE := function()
    if LocalConfig.cpuinfo.SIMD().hasSSE2() then
        return rec(
            2x32f := rec(
                wht := rec(
                    small := _defaultSizes(s->doSimdWht(s, SSE_2x32f, rec(propagateNth := true, useDeref := true, verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), [4]),
                    medium := _defaultSizes(s->doSimdWht(s, SSE_2x32f, rec(propagateNth := true, useDeref := true, oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([2..10], i->2^i))
                ),
                1d := rec(
                    dft_sc := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_2x32f, rec(propagateNth := true, useDeref := true, verify:=true, tsplBluestein:=false, interleavedComplex := true, PRDFT:=true, URDFT:= true, cplxVect := true, stdTTensor := false, globalUnrolling:=10000)), [ 2..32 ]),

                    ),
                    dft_ic := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_2x32f, rec(propagateNth := true, useDeref := true, verify:=true, tsplBluestein:=false, interleavedComplex := true, PRDFT:=true, URDFT:= true, cplxVect := true, stdTTensor := false, globalUnrolling:=10000)), [2..32])
                    )
                )
            ),
            2x64f := rec(
                wht := rec(
                    small := _defaultSizes(s->doSimdWht(s, SSE_2x64f, rec(verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), [4]),
                    medium := _defaultSizes(s->doSimdWht(s, SSE_2x64f, rec(oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([2..10], i->2^i))
                ),
                1d := rec(
                    dft_sc := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_2x64f, rec(verify:=true, tsplBluestein:=false, interleavedComplex := false, PRDFT:=true, URDFT:= true, cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                                [ 2..32 ]),
                        medium := _defaultSizes(s->doSimdDft(s, SSE_2x64f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)),
                                _svctSizes(1024, 16, 2)),
    #                    medium := spiral.libgen.doParSimdDft(SSE_2x64f, 1, _svctSizes(1024, 16, 4), false, true, false, false),
    #                    large := spiral.libgen.doParSimdDft(SSE_2x64f, 1, List([4..20], i->2^i), false, true, false, false)

                    ),
                    dft_ic := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_2x64f, rec(verify:=true, tsplBluestein:=false, interleavedComplex := true, PRDFT:=true, URDFT:= true, cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                                [2..32]),
                        medium := _defaultSizes(s->doSimdDft(s, SSE_2x64f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true)),
                                _svctSizes(1024, 16, 2)),
                        medium_cx := _defaultSizes(s->doSimdDft(s, SSE_2x64f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true, cplxVect := true, realVect := false)),
                                _svctSizes(1024, 16, 2)),
    #                    medium := spiral.libgen.doParSimdDft(SSE_2x64f, 1, _svctSizes(1024, 16, 4), false, true, false, true),
    #                    large := spiral.libgen.doParSimdDft(SSE_2x64f, 1, List([4..20], i->2^i), false, true, false, true)
                    ),
                    trdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, SSE_2x64f, rec( verify:=true,
                                interleavedComplex := true, PRDFT:=true, URDFT:= true, tsplBluestein := false, cplxVect := true,
                                realVect := true, propagateNth := true, useDeref := true,
                                globalUnrolling:=10000)), 4*[1..32]),

                    dht := _defaultSizes(s->doSimdSymDFT(TDHT, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    mdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, SSE_2x64f, rec(verify :=true)), 8*[1..12]),
                    imdct := _defaultSizes(s->doSimdSymDFT(TIMDCT, s, SSE_2x64f, rec(verify :=true)), 8*[1..12])
                ),
                2d := rec(
                    dft_ic := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_2x64f, rec(interleavedComplex := true,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    4*List([1..16], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_2x64f, rec(verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dft_sc := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_2x64f, rec(interleavedComplex := false,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    4*List([1..16], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_2x64f, rec(verify:=true, interleavedComplex := false, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dct2 := _defaultSizes(s->doSimdSymMDDFT(DCT2, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16]),
                    dct3 := _defaultSizes(s->doSimdSymMDDFT(DCT3, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16]),
                    dct4 := _defaultSizes(s->doSimdSymMDDFT(DCT4, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16]),
                    dst2 := _defaultSizes(s->doSimdSymMDDFT(DST2, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16]),
                    dst3 := _defaultSizes(s->doSimdSymMDDFT(DST3, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16]),
                    dst4 := _defaultSizes(s->doSimdSymMDDFT(DST4, s, SSE_2x64f, rec(verify := true)), [4, 8, 12, 16])
                )
            ),
            4x32i := rec(
                wht := rec(
                    small := _defaultSizes(s->doSimdWht(s, SSE_4x32i, rec(verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), List([], i->2^i)),
                    medium := _defaultSizes(s->doSimdWht(s, SSE_4x32i, rec(oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([4..10], i->2^i))
                )
            ),
            4x32f := rec(
                wht := rec(
                    small := _defaultSizes(s->doSimdWht(s, SSE_4x32f, rec(verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), List([1..3], i->2^i)),
                    medium := _defaultSizes(s->doSimdWht(s, SSE_4x32f, rec(oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([4..10], i->2^i))
                ),
                1d := rec(
                    dft_sc := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_4x32f, rec(verify:=true, interleavedComplex := false, stdTTensor := false, globalUnrolling:=10000)),
                                [ 2..64 ]),
                        medium := _defaultSizes(s->doSimdDft(s, SSE_4x32f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)),
                                _svctSizes(1024, 16, 4)),
    #                    medium := spiral.libgen.doParSimdDft(SSE_4x32f, 1, _svctSizes(1024, 16, 4), false, true, false, false),
    #                    large := spiral.libgen.doParSimdDft(SSE_4x32f, 1, List([4..20], i->2^i), false, true, false, false)

                    ),
                    dft_ic := rec(
                        small := _defaultSizes(s->doSimdDft(s, SSE_4x32f, rec(verify:=true, interleavedComplex := true, PRDFT:=true, URDFT:= true, cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                                [ 2..64 ]),
                        medium := _defaultSizes(s->doSimdDft(s, SSE_4x32f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true)),
                                _svctSizes(1024, 16, 4)),
                        medium_cx := _defaultSizes(s->doSimdDft(s, SSE_4x32f, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true,
                                cplxVect := true, realVect := false, PRDFT := false, URDFT := true)),
                                _svctSizes(1024, 16, 4))
    #                    medium := spiral.libgen.doParSimdDft(SSE_4x32f, 1, _svctSizes(1024, 16, 4), false, true, false, true),
    #                    large := spiral.libgen.doParSimdDft(SSE_4x32f, 1, List([4..20], i->2^i), false, true, false, true)
                    ),
                    rdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    trdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, SSE_4x32f, rec( verify:=true,
                                interleavedComplex := true, PRDFT:=true, URDFT:= true, tsplBluestein := false, cplxVect := true,
                                stdTTensor := false, realVect := true, propagateNth := true, useDeref := true,
                                globalUnrolling:=10000)), 8*[1..32]),
                    dht := _defaultSizes(s->doSimdSymDFT(TDHT, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    mdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, SSE_4x32f, rec(verify :=true)), 32*[1..12]),
                    imdct := _defaultSizes(s->doSimdSymDFT(TIMDCT, s, SSE_4x32f, rec(verify :=true)), 32*[1..12])
                ),
                2d := rec(
                    dft_ic := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_4x32f, rec(interleavedComplex := true,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    16*List([1..8], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_4x32f, rec(verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dft_sc := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_4x32f, rec(interleavedComplex := false,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    16*List([1..8], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_4x32f, rec(verify:=true, interleavedComplex := false, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dct2 := _defaultSizes(s->doSimdSymMDDFT(DCT2, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dct3 := _defaultSizes(s->doSimdSymMDDFT(DCT3, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dct4 := _defaultSizes(s->doSimdSymMDDFT(DCT4, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst2 := _defaultSizes(s->doSimdSymMDDFT(DST2, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst3 := _defaultSizes(s->doSimdSymMDDFT(DST3, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst4 := _defaultSizes(s->doSimdSymMDDFT(DST4, s, SSE_4x32f, rec(verify := true)), [4, 8, 12, 16]),
				# DO NOT USE BECAUSE IT GENERATES A NON_COMPATIBLE INIT() FUNCTION FOR USE WITH STANDARD TIMER
                #    conv := _defaultSizes(s->doSimdMdConv(s, SSE_4x32f, rec(svct := true, oddSizes := false,
				#					splitComplexTPrm := true, TRCDiag_VRCLR := true, globalUnrolling := 150, measureFinal := false)),
				#					16*[1..32])
                )
            ),
            8x16i := rec(
                wht := rec(
                    small := _defaultSizes(s->doSimdWht(s, SSE_8x16i, rec(verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), List([1..5], i->2^i)),
                    medium := _defaultSizes(s->doSimdWht(s, SSE_8x16i, rec(oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([6..10], i->2^i))
                ),
                1d := rec(
                    dft_ic := rec(
                        medium := _defaultSizes(s->doSimdDft(s, SSE_8x16i, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true, verify := true)), 64*[1..4]),
                        small := _defaultSizes(s->doSimdDft(s, SSE_8x16i, rec(verify:=true, interleavedComplex := true, stdTTensor := false, globalUnrolling:=10000)), [2..16])
                    ),
                    dft_sc := rec(
                        medium := _defaultSizes(s->doSimdDft(s, SSE_8x16i, rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false, verify := true)), 64*[1..4]),
                        small := _defaultSizes(s->doSimdDft(s, SSE_8x16i, rec(verify:=true, interleavedComplex := false, stdTTensor := false, globalUnrolling:=10000)), [2..16])
                    ),
                    rdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, SSE_8x16i, rec(verify := true)), 128*[1..4]),
                    dht := _defaultSizes(s->doSimdSymDFT(TDHT, s, SSE_8x16i, rec(verify :=true)), 128*[1..4]),
                    dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, SSE_8x16i, rec(verify := true, fracbits := 8)), [128]),
                    dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, SSE_8x16i, rec(verify := true, fracbits := 8)), [128]),
                    dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, SSE_8x16i, rec(verify := true)), 128*[1..4]),
                    dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, SSE_8x16i, rec(verify := true)), 128*[1..4]),
                    dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, SSE_8x16i, rec(verify := true)), 128*[1..4]),
                    dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, SSE_8x16i, rec(verify := true)), 128*[1..4]),
                    mdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, SSE_8x16i, rec(verify :=true)), 128*[1..4]),
                    imdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, SSE_8x16i, rec(verify :=true)), 128*[1..4])
                ),
                2d := rec(
                    dft_ic := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_8x16i, rec(interleavedComplex := true,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    64*List([1..2], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_8x16i, rec(verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dft_sc := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, SSE_8x16i, rec(interleavedComplex := false,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    64*List([1..2], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, SSE_8x16i, rec(verify:=true, interleavedComplex := false, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dct2 := _defaultSizes(s->doSimdSymMDDFT(DCT2, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16]),
                    dct3 := _defaultSizes(s->doSimdSymMDDFT(DCT3, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16]),
                    dct4 := _defaultSizes(s->doSimdSymMDDFT(DCT4, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16]),
                    dst2 := _defaultSizes(s->doSimdSymMDDFT(DST2, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16]),
                    dst3 := _defaultSizes(s->doSimdSymMDDFT(DST3, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16]),
                    dst4 := _defaultSizes(s->doSimdSymMDDFT(DST4, s, SSE_8x16i, rec(verify := true, fracbits := 10)), [8, 16])
                )
            )
        );
    else
        return false;
    fi;
end;


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Class(SymSPL, BaseContainer, SumsBase, rec(
    isBlock:=true,
    transpose := self >> self,
    rng:=self>>self._children[1].rng(),
    dmn:=self>>self._children[1].dmn(),
    vcost := self >> self._children[1].vcost()
));

ext: gi


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#D TPrm.getTags := self >> GetTags(self);
#D TPrm.firstTag := self >> let(tags := GetTags(self), When(tags=[], ANoTag, tags[1]));
#D TPrm.withTags := (self, tags) >> AddTag(self, tags);

# Breakdowns for vectorized TPrm( IJ )
#
NewRulesFor(TPrm, rec(
    # IP(N, any perm of size 1) -> I(N)
    TPrm_IP_Base1 := rec(
        applicable := (self,t) >> t.hasTags() and let(v := t.firstTag().v, p := t.params[1],
            ObjId(p) = IP and Rows(p.params[2]) = 1),

        apply := (self, t, C, Nonterms) >> I(Rows(t))
    ),

    TPrm_IJ_Vec := rec(
        requiredFirstTag := AVecReg,
        applicable := (self,t) >> t.hasTags() and let(v := t.firstTag().v, p := t.params[1],
            ObjId(p) = IP and ObjId(p.params[2]) = J and
            (Rows(p.params[2]) mod v) = 0
        ),

        children := (self, t) >> let(tags := t.getTags(), v := tags[1].v,
            [[ TPrm(J(v)).withTags(tags) ]]),

        apply := (self, t, C, Nonterms) >> let(
            v := t.firstTag().v,
            N := Rows(t),
            m := Rows(t.params[1].params[2]),
            k := N/m,
            j := Ind(N/v),

# NOTE: Pulling in breaks for vector code due to some missing comutative rule to pull it all over PushLR() -> VContainer() for now
            VContainer(IDirSum(j,
                COND(fCompose(diagTensor(II(k/2), II(2,0,1), II(m/v)), fBase(j)),
                    VTensor(Prm(fId(1)), v), C[1])) *
                    VTensor(Prm(condIJ(N/v, m/v)), v), t.firstTag().isa)
        )
    ),
    TPrm_format := rec(
        applicable := (self,t) >> t.hasTags() and t.hasTag(AVecReg),
        apply := (self, t, C, Nonterms) >> FormatPrm(t.params[1])
    )
));

# Breakdowns for vectorized TPrm( J )
#
NewRulesFor(TPrm, rec(
    TPrm_J := rec(
        requiredFirstTag := AVecReg,

        applicable := (self,t) >> t.hasTags() and let(p := t.params[1], v := t.firstTag().v, s := Rows(p),
            ObjId(p) = J and IsInt(s/v) and s<>v),

        children := (self, t) >> let(tags := t.getTags(), v := t.firstTag().v,
            [[ TPrm(J(v)).withTags(tags) ]]),

        apply := (self, t, C, Nonterms) >> let(v:=t.firstTag().v, s:=t.params[1].params[1], s2:=s/v,
                                    VTensor(J(s2), v) * Tensor(I(s2), C[1]))
    ),

    TPrm_Jv := rec(
        requiredFirstTag := AVecReg,

        applicable := (self,t) >> t.hasTags() and let(p := t.params[1], v := t.firstTag().v,
            ObjId(p) = J and Rows(p) = v),

        apply := (self, t, C, Nonterms) >> let(
            v := t.firstTag().v, isa := t.firstTag().isa,
            VPerm(J(v), isa.reverse, v, 0))
    ),
));


text:
Class(Dot_X, BaseOperation, rec(
  __call__ := (self, sr, x_n, s1, s2, range) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    start1 := s1,
    start2 := s2,
    range := range,
    x_n := x_n,
    _children := [range]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.start1, ", ", self.start2, ", ", self.range, ")"
  )
));

Class(Accum_X, BaseOperation, rec(
  __call__ := (self, sr, x_n, var, st, range, func) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    var := var,
    start := st,
    range := range,
    func := self >> self._children[1],
    x_n := x_n,
    _children := [func, range]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.var, ", ", self.start, ", ", self.range, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
	)
));

Class(Scat_E, BaseOperation, rec(
  __call__ := (self, sr, x_n, v1, v2, range, svar, func) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    v1 := v1,
    v2 := v2,
    range := range,
    svar := svar,
    func := self >> self._children[1],
    x_n := x_n,
    _children := [func, range]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.v1, ", ", self.v2, ", ", self.range, ", ", self.svar, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(SMP_Scat_E, BaseOperation, rec(
  __call__ := (self, sr, x_n, v1, v2, range, svar, func, nthreads) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    v1 := v1,
    v2 := v2,
    range := range,
    svar := svar,
    func := self >> self._children[1],
    x_n := x_n,
    nthreads := nthreads,
    _children := [func, range]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.v1, ", ", self.v2, ", ", self.range, ", ", self.svar, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(CheckEdgesUpdate, BaseOperation, rec(
  __call__ := (self, edges, var, func) >> SPL(WithBases(self, rec(
    edges := edges,
    var := var,
    func := self >> self._children[1],
    _children := [func]
  ))),
  dims := self >> [self.edges, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.edges, ", ", self.var, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(PruneEdges, BaseOperation, rec(
  __call__ := (self, idx, range, var, x_n, s, k) >> SPL(WithBases(self, rec(
    idx := idx,
    range := range,
    var := var,
    x_n := x_n,
    s := s,
    k := k,
    _children := [range, k]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.idx, ", ", self.range, ", ", self.var, ", ", self.s, ", ", self.k, ")"
  )
));

Class(Accum, BaseOperation, rec(
  __call__ := (self, sr, var, range_min, range_max, func) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    var := var,
    range_min := range_min,
    range_max := range_max,
    func := self >> self._children[1],
    _children := [func, range_min, range_max]
  ))),
  dims := self >> self._children[1].dims(),
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.var, ", ", self.range_min, ", ", self.range_max, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(SMP_Accum, BaseOperation, rec(
  __call__ := (self, sr, var, range_min, range_max, func, nthreads) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    var := var,
    range_min := range_min,
    range_max := range_max,
    func := self >> self._children[1],
    nthreads := nthreads,
    _children := [func, range_min, range_max]
  ))),
  dims := self >> self._children[1].dims(),
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.var, ", ", self.range_min, ", ", self.range_max, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(Assign_Ptr, BaseOperation, rec(
  __call__ := (self, func) >> SPL(WithBases(self, rec(
    func := self >> self._children[1],
    _children := [func]
  ))),
  dims := self >> self._children[1].dims(),
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(Scale_Pt, BaseOperation, rec(
  __call__ := (self, scalar, func) >> SPL(WithBases(self, rec(
    scalar := scalar,
    func := self >> self._children[1],
    _children := [func]
  ))),
  dims := self >> self._children[1].dims(),
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.scalar, ",\n",
        Blanks(i+is), self.func().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(IterateIndefinite, BaseOperation, rec(
  __call__ := (self, op1, op2) >> SPL(WithBases(self, rec(
    op1 := self >> self._children[1],
    op2 := self >> self._children[2],
    _children := [op1, op2]
  ))),
  dims := self >> [V(1), V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(\n",
        Blanks(i+is), self.op1().print(i+is, is), "\n",
        Blanks(i+is), self.op2().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(KTrussWrap, BaseOperation, rec(
  __call__ := (self, s, body) >> SPL(WithBases(self, rec(
    s := s,
    body := self >> self._children[1],
    _children := [body]
  ))),
  dims := self >> [V(1), V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.s, "\n",
        Blanks(i+is), self.body().print(i+is, is), "\n",
        Blanks(i), ")"
  )
));

Class(loopw, Command, rec(
  __call__ := (self, cond, cmd) >> WithBases(self, rec(
    cond := toExpArg(cond),
    cmd := Checked(IsCommand(cmd), cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.cond, self.cmd],
  rSetChild := rSetChildFields("cond", "cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.cond, ",\n",
        Blanks(i+si), self.cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(loopf, Command, rec(
  __call__ := (self, var, range_min, range_max, cmd) >> WithBases(self, rec(
    var := var,
    range_min := range_min,
    range_max := range_max,
    cmd := Checked(IsCommand(cmd), cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.var, self.range_min, self.range_max, self.cmd],
  rSetChild := rSetChildFields("var", "range_min", "range_max", "cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.var, ", ", self.range_min, ", ", self.range_max, ",\n",
        Blanks(i+si), self.cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(smp_reduct_loopf, Command, rec(
  __call__ := (self, sr, var, range_min, range_max, y, nthreads, cmd) >> WithBases(self, rec(
    sr := sr,
    var := var,
    range_min := range_min,
    range_max := range_max,
    y := y,
    nthreads := nthreads,
    cmd := Checked(IsCommand(cmd), cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.sr, self.var, self.range_min, self.range_max, self.y, self.nthreads, self.cmd],
  rSetChild := rSetChildFields("sr", "var", "range_min", "range_max", "y", "nthreads", "cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.var, ", ", self.range_min, ", ", self.range_max, ", ", self.nthreads, ",\n",
        Blanks(i+si), self.cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(smp_loopf, Command, rec(
  __call__ := (self, var, range_min, range_max, nthreads, cmd) >> WithBases(self, rec(
    var := var,
    range_min := range_min,
    range_max := range_max,
    nthreads := nthreads,
    cmd := Checked(IsCommand(cmd), cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.var, self.range_min, self.range_max, self.nthreads, self.cmd],
  rSetChild := rSetChildFields("var", "range_min", "range_max", "nthreads", "cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.var, ", ", self.range_min, ", ", self.range_max, ", ", self.nthreads, ",\n",
        Blanks(i+si), self.cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(if1, Command, rec(
  __call__ := (self, if_cond, if_cmd) >> WithBases(self, rec(
    if_cond := toExpArg(if_cond),
    if_cmd := Checked(IsCommand(if_cmd), if_cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.if_cond, self.if_cmd],
  rSetChild := rSetChildFields("if_cond", "if_cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.if_cond, ",\n",
        Blanks(i+si), self.if_cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(if2, Command, rec(
  __call__ := (self, if_cond, if_cmd, ei_cond, ei_cmd, else_cmd) >> WithBases(self, rec(
    if_cond := toExpArg(if_cond),
    if_cmd := Checked(IsCommand(if_cmd), if_cmd),
    else_cmd := Checked(IsCommand(else_cmd), else_cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.if_cond, self.if_cmd, self.else_cmd],
  rSetChild := rSetChildFields("if_cond", "if_cmd", "else_cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.if_cond, ",\n",
        Blanks(i+si), self.if_cmd.print(i+si, si), ",\n",
        Blanks(i+si), self.else_cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(if3, Command, rec(
  __call__ := (self, if_cond, if_cmd, ei_cond, ei_cmd, else_cmd) >> WithBases(self, rec(
    if_cond := toExpArg(if_cond),
    if_cmd := Checked(IsCommand(if_cmd), if_cmd),
    ei_cond := toExpArg(ei_cond),
    ei_cmd := Checked(IsCommand(ei_cmd), ei_cmd),
    else_cmd := Checked(IsCommand(else_cmd), else_cmd),
    operations := CmdOps
  )),
  rChildren := self >> [self.if_cond, self.if_cmd, self.ei_cond, self.ei_cmd, self.else_cmd],
  rSetChild := rSetChildFields("if_cond", "if_cmd", "ei_cond", "ei_cmd", "else_cmd"),
  print := (self, i, si) >> Print(self.__name__, "(", self.if_cond, ",\n",
        Blanks(i+si), self.if_cmd.print(i+si, si), ",\n",
        Blanks(i+si), self.ei_cond, ",\n",
        Blanks(i+si), self.ei_cmd.print(i+si, si), ",\n",
        Blanks(i+si), self.else_cmd.print(i+si, si), "\n",
        Blanks(i), ")"
  )
));

Class(break_ic, Command, rec(
  __call__ := (self) >> WithBases(self, rec(
    operations := CmdOps
  )),
  rChildren := self >> [],
  rSetChild := rSetChildFields(),
  print := (self, i, si) >> Print(self.__name__, "()")
));

CUnparser.loopw := (self,o,i,is) >> Print(Blanks(i),
    "while (", self(o.cond,i,is), ") {\n",
    self(o.cmd,i+is,is), Blanks(i), "}\n");

CUnparser.loopf := (self,o,i,is) >> Print(Blanks(i),
    "for (int ", o.var, " = ", self(o.range_min,i,is), "; ", o.var, " < ", self(o.range_max,i,is), "; ", o.var, "++) {\n",
    self(o.cmd,i+is,is), Blanks(i), "}\n");

CUnparser.smp_reduct_loopf := (self,o,i,is) >> Print(Blanks(i),
    "#pragma omp parallel for schedule(dynamic) reduction(", o.sr.OpenMP_sum, ":", o.y, ") num_threads(", o.nthreads, ")\n",
    Blanks(i), "for (int ", o.var, " = ", self(o.range_min,i,is), "; ", o.var, " < ", self(o.range_max,i,is), "; ", o.var, "++) {\n",
    self(o.cmd,i+is,is), Blanks(i), "}\n");

CUnparser.smp_loopf := (self,o,i,is) >> Print(Blanks(i),
    "#pragma omp parallel for schedule(dynamic) num_threads(", o.nthreads, ")\n",
    Blanks(i), "for (int ", o.var, " = ", self(o.range_min,i,is), "; ", o.var, " < ", self(o.range_max,i,is), "; ", o.var, "++) {\n",
    self(o.cmd,i+is,is), Blanks(i), "}\n");

CUnparser.if1 := (self,o,i,is) >> Print(Blanks(i),
    "if (", self(o.if_cond,i,is), ") {\n", self(o.if_cmd,i+is,is), Blanks(i), "}\n");

CUnparser.if2 := (self,o,i,is) >> Print(Blanks(i),
    "if (", self(o.if_cond,i,is), ") {\n", self(o.if_cmd,i+is,is), Blanks(i), "}",
    " else {\n", self(o.else_cmd,i+is,is), Blanks(i), "}\n");

CUnparser.if3 := (self,o,i,is) >> Print(Blanks(i),
    "if (", self(o.if_cond,i,is), ") {\n", self(o.if_cmd,i+is,is), Blanks(i), "}",
    " else if (", self(o.ei_cond,i,is), ") {\n", self(o.ei_cmd,i+is,is), Blanks(i), "}",
    " else {\n", self(o.else_cmd,i+is,is), Blanks(i), "}\n");

CUnparser.break_ic := (self,o,i,is) >> Print(Blanks(i),
    "break;\n");

Class(VMV_X, TaggedNonTerminal, rec(
  __call__ := (self, sr, s1, sm, s2, d1, d2, x_n) >> WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    start1 := s1,
    startm := sm,
    start2 := s2,
    d1 := d1,
    d2 := d2,
    x_n := x_n,
    params := []
  )),
  isReal := self >> true,
  dims := self >> [self.x_n, V(1)]
));

Class(TriangleCount, TaggedNonTerminal, rec(
  abbrevs := [ (n) -> [n] ],
  isReal := self >> true,
  dims := self >> [self.params[1], self.params[1]]
));

Class(KTruss, TaggedNonTerminal, rec(
  abbrevs := [ (n, m, k) -> [n, m, k] ],
  isReal := self >> true,
  dims := self >> [self.params[1], self.params[1]]
));

Class(TriangleCountEdges, TaggedNonTerminal, rec(
  abbrevs := [ (n, m, s) -> [n, m, s] ],
  isReal := self >> true,
  dims := self >> [self.params[1], self.params[1]]
));

Class(KTrussPrune, TaggedNonTerminal, rec(
  abbrevs := [ (n, m, k, s) -> [n, m, k, s] ],
  isReal := self >> true,
  dims := self >> [self.params[1], self.params[1]]
));

NewRulesFor(VMV_X, rec(
  Accum_X_Dot := rec(
    info := "Selective accumulation of dot products <Mi,v2> over existing elements in v1",
    maxSize := false,
    applicable := (self, t) >> true,
    apply := (nt, C, ct) -> let(i := Ind(), sr := nt.sr, x_n := nt.x_n,
              Accum_X(sr, x_n, i, nt.start1, nt.d1, Dot_X(sr, x_n, [i, nt.startm[2]], nt.start2, nt.d2)))
  )
));

NewRulesFor(TriangleCount, rec(
  Accum_VMV_FLAME1 := rec(
    info := "FLAME Approach #1",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    children := (self, nt) >> let(i := Ind(), n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              [[ VMV_X(sr, [i, V(0)], [V(0), V(0)], [i, V(0)], i, i, n) ]]),
    apply := (nt, C, ct) -> let(i := C[1]._children[1].rt.node.d1, n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              Assign_Ptr(Scale_Pt(V(0.5), Accum(sr, i, V(1), n, C[1]))))
  ),
  Accum_VMV_FLAME2 := rec(
    info := "FLAME Approach #2",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    children := (self, nt) >> let(i := Ind(), n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              [[ VMV_X(sr, [i, V(0)], [V(0), i+1], [i, i+1], i, n-i-1, n) ]]),
    apply := (nt, C, ct) -> let(i := C[1]._children[1].rt.node.d1, n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              Assign_Ptr(Accum(sr, i, V(1), n-1, C[1])))
  ),
  Accum_VMV_FLAME4 := rec(
    info := "FLAME Approach #4",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    children := (self, nt) >> let(i := Ind(), n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              [[ VMV_X(sr, [i-1, i], [i, i], [i-1, i], n-i, n-i, n) ]]),
    apply := (nt, C, ct) -> let(i := C[1]._children[1].rt.node.startm[1], n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              Assign_Ptr(Scale_Pt(V(0.5), Accum(sr, i, V(1), n, C[1]))))
  ),
  Accum_VMV_FLAME2_SMP := rec(
    info := "FLAME Approach #2",
    maxSize := false,
    applicable := nt -> nt.hasTags(),
    children := (self, nt) >> let(i := Ind(), n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              [[ VMV_X(sr, [i, V(0)], [V(0), i+1], [i, i+1], i, n-i-1, n) ]]),
    apply := (nt, C, ct) -> let(i := C[1]._children[1].rt.node.d1, n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              Assign_Ptr(SMP_Accum(sr, i, V(1), n-1, C[1], nt.firstTag().params[1])))
  )
));

NewRulesFor(KTruss, rec(
  TC_Support_Prune := rec(
    info := "Triangle Count for Support, then Prune. Iterate",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    children := (self, nt) >> let(n := nt.params[1], m := nt.params[2], k := nt.params[3], s := var.fresh_t("s", TArray(TInt, nt.params[2])),
              [[ TriangleCountEdges(n, m, s), KTrussPrune(n, m, k, s) ]]),
    apply := (nt, C, ct) -> Assign_Ptr(Scale_Pt(V(0.5), IterateIndefinite(C[1], C[2])))
  ),
  TC_Support_Prune_SMP := rec(
    info := "Triangle Count for Support, then Prune. Iterate",
    maxSize := false,
    applicable := nt -> nt.hasTags(),
    children := (self, nt) >> let(n := nt.params[1], m := nt.params[2], k := nt.params[3], s := var.fresh_t("s", TArray(TInt, nt.params[2])),
              [[ TriangleCountEdges(n, m, s).withTags(nt.tags), KTrussPrune(n, m, k, s) ]]),
    apply := (nt, C, ct) -> Assign_Ptr(Scale_Pt(V(0.5), IterateIndefinite(C[1], C[2])))
  )
));

NewRulesFor(TriangleCountEdges, rec(
  Endpoint_Row_Dot := rec(
    info := "Dot products of rows corresponding to endpoints of each edge",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    apply := (nt, C, ct) -> let(v1 := var.fresh_t("v", TInt), v2 := var.fresh_t("v", TInt), n := nt.params[1], m := nt.params[2], s := nt.params[3], sr := TSemiring_Arithmetic(n.t),
              Scat_E(sr, n, v1, v2, m, s, Dot_X(sr, n, [v1, 1], [v2, 1], n)))
  ),
  Endpoint_Row_Dot_SMP := rec(
    info := "Dot products of rows corresponding to endpoints of each edge",
    maxSize := false,
    applicable := nt -> nt.hasTags(),
    apply := (nt, C, ct) -> let(v1 := var.fresh_t("v", TInt), v2 := var.fresh_t("v", TInt), n := nt.params[1], m := nt.params[2], s := nt.params[3], sr := TSemiring_Arithmetic(n.t),
              SMP_Scat_E(sr, n, v1, v2, m, s, Dot_X(sr, n, [v1, 1], [v2, 1], n), nt.firstTag().params[1]))
  )
));

NewRulesFor(KTrussPrune, rec(
  Check_Support_Prune := rec(
    info := "Eliminate edges without enough support and break if nothing changes",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    apply := (nt, C, ct) -> let(i := Ind(), v := var.fresh_t("v", TInt), n := nt.params[1], m := nt.params[2], k := nt.params[3], s := nt.params[4],
              CheckEdgesUpdate(m, v, PruneEdges(i, m, v, n, s, k)))
  )
));

DefaultSumsGen.Accum := (self, o, opts) >> Accum(o.sr, o.var, o.range_min, o.range_max, self(o.func(), opts));

DefaultSumsGen.Accum_X := (self, o, opts) >> Accum_X(o.sr, o.x_n, o.var, o.start, o.range, self(o.func(), opts));

DefaultSumsGen.Scat_E := (self, o, opts) >> Scat_E(o.sr, o.x_n, o.v1, o.v2, o.range, o.svar, self(o.func(), opts));

DefaultSumsGen.SMP_Scat_E := (self, o, opts) >> SMP_Scat_E(o.sr, o.x_n, o.v1, o.v2, o.range, o.svar, self(o.func(), opts), o.nthreads);

DefaultSumsGen.Dot_X := (self, o, opts) >> o;

DefaultSumsGen.CheckEdgesUpdate := (self, o, opts) >> CheckEdgesUpdate(o.edges, o.var, self(o.func(), opts));

DefaultSumsGen.PruneEdges := (self, o, opts) >> o;

DefaultSumsGen.SMP_Accum := (self, o, opts) >> SMP_Accum(o.sr, o.var, o.range_min, o.range_max, self(o.func(), opts), o.nthreads);

DefaultSumsGen.Assign_Ptr := (self, o, opts) >> Assign_Ptr(self(o.func(), opts));

DefaultSumsGen.Scale_Pt := (self, o, opts) >> Scale_Pt(o.scalar, self(o.func(), opts));

DefaultSumsGen.IterateIndefinite := (self, o, opts) >> IterateIndefinite(self(o.op1(), opts), self(o.op2(), opts));

DefaultCodegen.Dot_X := meth(self, o, y, x, opts)
  local j1, j1m, j2, j2m;
  if (IsBound(opts.isCSR) and opts.isCSR) then
    j1 := var.fresh_t("j1", x.t);
    j1m := var.fresh_t("j1m", x.t);
    j2 := var.fresh_t("j2", x.t);
    j2m := var.fresh_t("j2m", x.t);
    return decl([j1, j1m, j2, j2m], chain(
      assign(y, o.sr.zero()),
      assign(j1, add(x, o.x_n, V(1), nth(x, o.start1[1]))),
      assign(j1m, add(x, o.x_n, V(1), nth(x, add(o.start1[1], V(1))))),
      assign(j2, add(x, o.x_n, V(1), nth(x, o.start2[1]))),
      assign(j2m, add(x, o.x_n, V(1), nth(x, add(o.start2[1], V(1))))),
      loopw(logic_and(lt(j1, j1m), lt(deref(j1), o.start1[2])),
        assign(j1, add(j1, V(1)))),
      loopw(logic_and(lt(j2, j2m), lt(deref(j2), o.start2[2])),
        assign(j2, add(j2, V(1)))),
      loopw(logic_and(logic_and(lt(j1, j1m), lt(j2, j2m)), logic_and(lt(deref(j1), add(o.start1[2], o.range)), lt(deref(j2), add(o.start2[2], o.range)))),
        if3(lt(deref(j1), deref(j2)),
          assign(j1, add(j1, V(1))),
        lt(deref(j2), deref(j1)),
          assign(j2, add(j2, V(1))),
        chain(
          assign(y, o.sr.sum(y, V(1))),
          assign(j1, add(j1, V(1))),
          assign(j2, add(j2, V(1))))))));
  fi;
end;

DefaultCodegen.Accum_X := meth(self, o, y, x, opts)
  local j, jm, t;
  t := var.fresh_t("t", y.t);
  if (IsBound(opts.isCSR) and opts.isCSR) then
    j := var.fresh_t("j", x.t);
    jm := var.fresh_t("jm", x.t);
    return decl([j, jm], chain(
      assign(y, o.sr.zero()),
      assign(j, add(x, o.x_n, V(1), nth(x, o.start[1]))),
      assign(jm, add(x, o.x_n, V(1), nth(x, add(o.start[1], V(1))))),
      loopw(logic_and(lt(j, jm), lt(deref(j), o.start[2])),
        assign(j, add(j, V(1)))),
      loopw(logic_and(lt(j, jm), lt(deref(j), add(o.start[2], o.range))), decl([o.var, t], chain(
        assign(o.var, deref(j)),
        self(o.func(), t, x, opts),
        assign(y, o.sr.sum(y, t)),
        assign(j, add(j, V(1))))))));
  fi;
end;

DefaultCodegen.Scat_E := meth(self, o, y, x, opts)
  local t, i;
  t := var.fresh_t("t", y.t);
  i := Ind();
  if (IsBound(opts.isCSR) and opts.isCSR) then
    return decl([o.v1, o.svar], chain(
      assign(o.v1, V(0)),
      loopf(i, V(0), o.range, decl([t, o.v2], chain(
        loopw(eq(nth(x, add(o.v1, V(1))), i),
          assign(o.v1, add(o.v1, V(1)))),
        assign(o.v2, nth(x, add(o.x_n, V(1), i))),
        self(o.func(), t, x, opts),
        assign(nth(o.svar, i), t))))));
  fi;
end;

DefaultCodegen.SMP_Scat_E := meth(self, o, y, x, opts)
  local t, i, j;
  t := var.fresh_t("t", y.t);
  i := Ind();
  j := Ind();
  if (IsBound(opts.isCSR) and opts.isCSR) then
    return decl([o.v1, o.svar],
      smp_loopf(i, V(0), o.x_n, o.nthreads,
        loopf(j, nth(x, i), nth(x, add(i, V(1))), decl([o.v1, o.v2, t], chain(
          assign(o.v1, i),
          assign(o.v2, nth(x, add(o.x_n, V(1), j))),
          self(o.func(), t, x, opts),
          assign(nth(o.svar, j), t))))));
  fi;
end;

DefaultCodegen.CheckEdgesUpdate := meth(self, o, y, x, opts)
  return decl([o.var], chain(
    assign(o.var, V(0)),
    assign(y, o.edges),
    self(o.func(), y, x, opts),
    if1(eq(o.edges, o.var),
      break_ic()),
    assign(o.edges, o.var),
    assign(y, o.var)));
end;

DefaultCodegen.PruneEdges := meth(self, o, y, x, opts)
  local t;
  t := var.fresh_t("t", y.t);
  if (IsBound(opts.isCSR) and opts.isCSR) then
    return decl([t], chain(
      assign(t, V(0)),
      loopf(o.idx, V(0), o.range, chain(
        loopw(eq(nth(x, add(t, V(1))), o.idx), chain(
          assign(t, add(t, V(1))),
          assign(nth(x, t), o.var))),
        if1(geq(nth(o.s, o.idx), sub(o.k, V(2))), chain(
          assign(nth(x, add(o.x_n, V(1), o.var)), nth(x, add(o.x_n, V(1), o.idx))),
          assign(o.var, add(o.var, V(1))))))),
      assign(nth(x, add(t, V(1))), o.var)));
  fi;
end;

DefaultCodegen.Accum := meth(self, o, y, x, opts)
  local t;
  t := var.fresh_t("t", y.t);
  return chain(
    assign(y, o.sr.zero()),
    loopf(o.var, o.range_min, o.range_max, decl([t], chain(
      self(o.func(), t, x, opts),
      assign(y, o.sr.sum(y, t))))));
end;

DefaultCodegen.SMP_Accum := meth(self, o, y, x, opts)
  local t;
  t := var.fresh_t("t", y.t);
  return chain(
    assign(y, o.sr.zero()),
    smp_reduct_loopf(o.sr, o.var, o.range_min, o.range_max, y, o.nthreads, decl([t], chain(
      self(o.func(), t, x, opts),
      assign(y, o.sr.sum(y, t))))));
end;

DefaultCodegen.Scale_Pt := meth(self, o, y, x, opts)
  return chain(
    self(o.func(), y, x, opts),
    assign(y, mul(o.scalar, y)));
end;

DefaultCodegen.Assign_Ptr := meth(self, o, y, x, opts)
  local t;
  t := var.fresh_t("t", y.t.t);
  return decl([t], chain(
    self(o._children[1], t, x, opts),
    assign(deref(y), t)));
end;

DefaultCodegen.IterateIndefinite := meth(self, o, y, x, opts)
  return loopw(V(1), chain(
    self(o.op1(), y, x, opts),
    self(o.op2(), y, x, opts)));
end;

GraphIndicesCS := [
    c -> Compile.pullDataDeclsRefs(c),
    c -> Compile.fastScalarize(c),
    c -> UnrollCode(c),
    c -> FlattenCode(c),
    c -> UntangleChain(c),
    #(c, opts) -> CopyPropagate.initial(c, opts),
    #(c, opts) -> HashConsts(c, opts),
    #c -> MarkDefUse(c),
    #(c, opts) -> BinSplit(c, opts),
    #c -> MarkDefUse(c),
    #CopyPropagate, # does CSE
    c -> Compile.declareVars(c)
];

TCDefaults := CopyFields(SpiralDefaults, rec(
  compileStrategy := GraphIndicesCS,
  X := var("IJ", TPtr(TInt)),
  XType := TPtr(TInt),
  arrayDataModifier := "",
  arrayBufModifier := "",
  Y := var("res", TPtr(TInt)),
  YType := TPtr(TInt),
  isCSR := true,
));


text:
#
# walrus: Computational Methods for Finitely Generated Monoids and Groups
#

## Places
InstallGlobalFunction(NewPlace,
function(loc, c, colour)
    local p;
    p := Objectify(PregroupPlaceType, [loc,c,colour]);
    Add(loc![3], p);
    return p;
end);

InstallMethod(Location
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
             p -> p![1]);

InstallMethod(Relator
             , "for a pregroup place"
             , [IsPregroupPlace],
             p -> Relator(p![1]));

InstallMethod(Letter
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
             p -> p![2]);

InstallMethod(Colour
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
             p -> p![3]);

InstallMethod(NextPlaces
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
function(p)
    if not IsBound(p![6]) then
        p![6] := Places(NextLocation(Location(p)));
    fi;
    return p![6];
end);

AddOrUpdate := function(map, key, val)
    local tmp;

    tmp := map[key];
    if tmp = fail or val < tmp then
        map[key] := val;
    fi;
end;

InstallGlobalFunction(OneStepRedCase,
function(P)
    local Q, Ql
          , b, c, d, x, y
          , Lp
          , v, v1, v2
          , res
          , xi1, xi2, binv, l, onv, pres, vg, xi;

    pres := PregroupPresentationOf(Relator(P));
    res := HashMap(1024);
    vg := VertexGraph(pres);

    c := Letter(P);

    for Q in NextPlaces(P) do
        # same relator, one position up
        Ql := Location(Q);
        b := InLetter(Ql); binv := PregroupInverse(b);
        d := OutLetter(Ql);
        x := Letter(Q);
        v := VertexFor(vg, [b, d, 0]);
        onv := OutNeighboursOfVertex(vg, v);

        for y in IntermultMap(binv) do
            v1 := VertexFor(vg, [y, binv, 1]);
            xi1 := Blob(pres, y, binv, c);
            for v2 in onv do
                #T this check is new, and assuming the paper
                #T is accurate correct
                if DigraphVertexLabel(vg, v2)[1][2] = x then
                    xi2 := Vertex(pres, v1, v, v2);

                    AddOrUpdate(res, [ __ID(Q), 1 ], xi1 + xi2);
                fi;
            od;
        od;
    od;
    # Note this list is not necessarily dense atm.
    return res;
end);

InstallGlobalFunction(OneStepGreenCase,
function(P)
 local L, L2, b, c, loc, pls, is_consoledge, v, v1, v2, xi1, xi2,
       R, R2, P2, P2s, P2T, i, j, next, res, len, n, l, pres, vg,
       P2_loc, P2_inletter, P2_outletter, P2_outletterinv, P2_letter;

    pres := PregroupPresentationOf(Relator(P));
    res := HashMap(1024);
    vg := VertexGraph(pres);

    L := Location(P);
    b := OutLetter(L);
    c := Letter(P);

    # Every place that is instanciated with location that is in an instantiation of
    # a place P'.
    # We're interested in consolidated edges between Relator(P) and Relator(L2)
    # from the locations that P and pls are at.
    #
    # Do we have to check that Relator(P) and Relator(pls) don't entirely cancel?
    for L2 in Locations(pres) do
        R2 := Relator(L2);

        # L2 instantiates place on R2, have to test consolidated
        #    edges between R and R2 starting from L/L2 on R/R2
        #    respectively
        if (InLetter(L2) = PregroupInverse(b))
            and (OutLetter(L2) = c) then

            # We compute all consolidated edge places,
            # we could do this incrementally, but I don't
            # currently see a use in doing so
            P2s := ConsolidatedEdgePlaces(L, L2);

            for P2T in P2s do
                P2 := P2T[1];  # Place reachable on R1 by consolidated edge

                P2_loc := Location(P2);
                P2_inletter := InLetter(P2_loc);
                P2_outletter := OutLetter(P2_loc);
                P2_outletterinv := PregroupInverse(P2_outletter);
                P2_letter := Letter(P2);

                # P2T[2] location on R2 reachable by the edge
                len := P2T[3]; # length of consolidated edge
                v1 := VertexFor(vg, [ InLetter(P2T[2]), OutLetter(P2T[2]), 0 ]);
                v := VertexFor(vg, [ P2_inletter, P2_outletter, 0 ]);
                for v2 in OutNeighboursOfVertex(vg, v) do
                    if DigraphVertexLabel(vg, v2)[1][1] = P2_outletterinv and
                       P2_letter = DigraphVertexLabel(vg, v2)[1][2] then

                        xi1 := Vertex(pres, v1, v, v2);
                        if Colour(P2) = "green" then
                            AddOrUpdate(res, [ __ID(P2), len ], xi1);

                        elif Colour(P2) = "red" then
                            next := OneStepRedCase(P2);
                            # testhack
                            for n in Keys(next) do
                                AddOrUpdate( res, [ n[1], len + 1 ]
                                           , xi1 + next[n] );
                            od;
                        else
                            Error("Invalid colour");
                        fi;
                    else

                    fi;
                od;
            od;
        fi;
    od;
    return res;
end);

InstallMethod(OneStepReachablePlaces
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
function(p)
    if not IsBound(p![7]) then
        if Colour(p) = "red" then
            p![7] := OneStepRedCase(p);
        elif Colour(p) = "green" then
            p![7] := OneStepGreenCase(p);
        else
            Error("Invalid colour for place ", p, "\n");
        fi;
    fi;
    return p![7];
end);

InstallMethod(__ID
             , "for a pregroup place"
             , [IsPregroupPlace],
             p -> p![5]);

InstallMethod(ViewString
             , "for a pregroup place"
             , [IsPregroupPlaceRep],
function(p)
    return STRINGIFY("(", ViewString(p![1]),
                     ",", ViewString(p![2]),
                     ",", ViewString(p![3]),
                     ")");
end);

InstallMethod(\=
             , "for two pregroup places"
             , [IsPregroupPlaceRep, IsPregroupPlaceRep],
function(p1, p2)
    return (p1![1] = p2![1])
           and (p1![2] = p2![2])
           and (p1![3] = p2![3]);
end);

max_stars_repo_name: wilfwilson/walrus

max_issues_repo_name: wilfwilson/walrus

max_forks_repo_name: wilfwilson/walrus


text:
#############################################################################
##
#F  MergeAndGroupRecord(<dst>,<rec1>,<rec2>...)  . . . . . . . merge the fields of records. Same fields get grouped into a set.
##
MergeAndGroupRecords := function ( arg )
    local   res,        # merged record, result
            record,     # one of the arguments
            name,       # name of one component of <record>
            set_list;  # keeps track if a set was created for a field

    res := When(arg=[], rec(), arg[1]);
    set_list := [];
    for record  in  arg{[2..Length(arg)]} do
        for name  in UserRecFields( record )  do
            if IsBound(res.(name))  then
                if name in set_list then
                    res.(name) := UnionSet(res.(name), Set([record.(name)]) );
                else
                    res.(name) := UnionSet(Set([res.(name)]), Set([record.(name)]));
                    Add(set_list, name);
                fi;
            else
                res.(name) := record.(name);
            fi;
        od;
    od;
    return res;
end;

Class(CudaCodegen, DefaultCodegen, rec(

    SIMTISum := meth(self, o, y, x, opts)
        local lcode, simtidx, decls, d, v;

        simtidx := o.simt_dim.simt_idx();
        lcode := simt_loop(simtidx, o.var, o.domain,
            self(o.child(1), y, x, opts));

        # If not a regular loop, add simt index to inner array's parallel context
        if simtidx <> Ignore then
            self.simtidx_rec.(o.var.id) := simtidx;
            decls := Collect(lcode, @(1, decl, d -> ForAny(d.vars, v -> IsArrayT(v.t))));
            for d in decls do
                DoForAll(d.vars, v -> When(IsBound(v.parctx), Add(v.parctx, [o.var, simtidx]), v.setAttrTo("parctx", [ [o.var, simtidx] ])) );
            od;
        fi;

        return When( ObjId(o.simt_dim)=ASIMTKernelFlag, simt_block( lcode ), lcode );
    end,

    SIMTSUM := meth(self, o,y,x,opts)
        local i, ch, numch, code_ch, cch, pb, parb;
        ch    := o.children();
        numch := Length(ch);
        code_ch := chain();

        for i in [1..numch] do
            cch := self(ch[i], y, x, opts);

            if o.simt_dim.need_cond(i) then
                Add(code_ch.cmds, IF(o.simt_dim.simt_cond(i), cch, skip()) );
            else
                Add(code_ch.cmds, cch);
            fi;

        od;

        return code_ch;
    end,

    BB := meth(self,o,y,x,opts)
        local res;

        res := Inherited(o, y, x, opts);

        # inner temps are unrolled - drop from pre/post tables (See Compose below)
        if opts.mempool then
            self.tarray_pre  := DropLast(self.tarray_pre,  1);
            self.tarray_post := DropLast(self.tarray_post, 1);
            Add(self.tarray_pre,  rec());
            Add(self.tarray_post, rec());
        fi;

        return res;

    end,

    Compose := meth(self, o,y,x,opts)
        local ch, numch, vecs, code_ch, get_simt_sync, ta_pre, ta_post, ci, cch, simts, parb;
        ch    := o.children();
        numch := Length(ch);

        # propagate x and y arrays and create temporary arrays
        vecs  := self._composePropagate(ch, y, x, c -> TempArray(y,x,c));

        # order them so that first to be evaluated is first in array.
        vecs := Reversed(vecs);
        ch   := Reversed(ch);

        get_simt_sync := (c) -> let(simtss := List( Collect(c, @.cond(e -> IsRec(e) and IsBound(e.simt_dim))), e -> e.simt_dim.simt_sync() ),
                                   When(simtss<>[], Maximum(simtss), skip() )
                                );

        code_ch := [];

        if opts.mempool then
            ta_pre  := Last(self.tarray_pre);
            ta_post := Last(self.tarray_post);
        fi;

        for ci in [1..numch] do
            if opts.mempool then
                Add(self.tarray_pre,  rec());
                Add(self.tarray_post, rec());
            fi;
            cch := self(ch[ci], vecs[ci+1], vecs[ci], opts);
            # Determine sync primitive at ch's boundary
            simts := get_simt_sync(ch[ci]);
            Append(code_ch, When(ObjId(simts)<>skip, [cch, simts], [cch]));

            if opts.mempool then
                # Add new temps to pre/post tables.
                # Composition creates a pre/post relationship btw temps.
                if Last(self.tarray_pre) = rec() then
                    ta_pre  := MergeAndGroupRecords(ta_pre,  rec( (vecs[ci+1].id) := Set([vecs[ci]]) ) );
                    ta_post := MergeAndGroupRecords(ta_post, rec( (vecs[ci].id)   := Set([vecs[ci+1]]) ) );
                else
                    ta_pre  := MergeAndGroupRecords(ta_pre,  Last(self.tarray_pre) );
                    ta_post := MergeAndGroupRecords(ta_post, Last(self.tarray_post) );
                fi;
                self.tarray_pre  := DropLast(self.tarray_pre,  1);
                self.tarray_post := DropLast(self.tarray_post, 1);

            fi;
        od;

        # Wrap code in variable declaration.
        return decl(
                Difference(Flat(vecs{[2..Length(vecs)-1]}), Flat([x,y])),
                chain( code_ch )
            );
    end,

    wrap_gpu_timing := function(ker_list, opts)
        local start, stop, ms;

        if not opts.gpu_timing then
            return ker_list;
        fi;

        Append(opts.includes, ["<stdio.h>", "<helper_cuda.h>"]);

        start := var.fresh_t("start", T_CudaEvent);
        stop := var.fresh_t("stop", T_CudaEvent);
        ms    := var.fresh_t("ms", T_Real(32));

        return decl([start, stop, ms],
                    chain(
                        cu_event_create(addrof(start)),
                        cu_event_create(addrof(stop)),
                        cu_check_errors(cu_event_record(start, false)),
                        ker_list,
                        cu_check_errors(cu_event_record(stop, false)),
                        cu_check_errors(cu_event_synchronize(stop, false)),
                        assign(ms, V(0)),
                        cu_event_elapsed_time(addrof(ms), start, stop),
                        #cprintf("SPIRAL-based CUDA execution [ms]:\\t%f\\n", ms),
                        cu_event_destroy(start),
                        cu_event_destroy(stop)
                    )
            );
    end,

    wrap_kernel := meth(self, ker_desc, io, params, o_dims, opts)
        local cp_in, cp_out;
        cp_out := cu_memcpy(var.table.("h"::io[1].id), io[1], o_dims[1], cudaMemcpyDeviceToHost);
        cp_in  := [cu_memcpy(io[2], var.table.("h"::io[2].id), o_dims[2], cudaMemcpyHostToDevice)]
                ::List(params, p -> cu_memcpy(p, var.table.("h"::p.id), p.t.size, cudaMemcpyHostToDevice));

        return decl(Set(List(ker_desc.cuda_kers, kc ->kc.dim_grid))::Set(List(ker_desc.cuda_kers, kc ->kc.dim_block)),
                chain(
                    cp_in,
                    self.wrap_gpu_timing(
                        List(ker_desc.cuda_kers, kc -> ApplyFunc(cu_call,[kc.cuda_ker.id, kc.dim_grid, kc.dim_block]::kc.cuda_ker.params) ),
                        opts),
                    cp_out
                ));
    end,

    wrap_kernel_devfunc := meth(self, ker_desc, io, params, o_dims, opts)
        return decl(Set(List(ker_desc.cuda_kers, kc ->kc.dim_grid))::Set(List(ker_desc.cuda_kers, kc ->kc.dim_block)),
                chain(
                        List(ker_desc.cuda_kers, kc -> ApplyFunc(cu_call,[kc.cuda_ker.id, kc.dim_grid, kc.dim_block]::kc.cuda_ker.params) )
                ));
    end,

    get_dim_idx := (self, ker, IdxT) >> let( ll := List(Collect(ker, @(1,simt_loop,l->ObjId(l.simt_idx)=IdxT) ), l->l.simt_idx ),
                                         When(Length(ll)>0,
                                            Checked(ForAll(ll, l -> l.size()=ll[1].size()), ll[1].size()),
                                            V(1)
                                        )),

    _reduce_tarray_graph := function(pre, post, ta_list)
        local ta_id, ta, t, new_pre, new_post;
        new_pre  := Copy(pre);
        new_post := Copy(post);
        # Only consider provided arrays
        for ta_id in UnionSet(UserRecFields(pre), UserRecFields(post)) do
            ta := var.table.(ta_id);
            if not ta in ta_list then
                if IsBound(new_pre.(ta_id)) then
                    for t in new_pre.(ta_id) do
                        new_post.(t.id) := UnionSet( When(IsBound(new_post.(t.id)),  new_post.(t.id),  []),
                                                     When(IsBound(new_post.(ta_id)), new_post.(ta_id), []) );
                        new_post.(t.id) := Difference(new_post.(t.id), [ta]);
                    od;
                fi;
                if IsBound(new_post.(ta_id)) then
                    for t in new_post.(ta_id) do
                        new_pre.(t.id) := UnionSet( When(IsBound(new_pre.(t.id)),  new_pre.(t.id),  []),
                                                    When(IsBound(new_pre.(ta_id)), new_pre.(ta_id), []) );
                        new_pre.(t.id) := Difference(new_pre.(t.id), [ta]);
                    od;
                fi;

                if IsBound(new_pre.(ta_id)) then
                    Unbind(new_pre.(ta_id));
                fi;
                if IsBound(new_post.(ta_id)) then
                    Unbind(new_post.(ta_id));
                fi;
            fi;
        od;

        return [new_pre, new_post];
    end,

    _next_color := function(i, e, exclude)
        local c;
        c := i;
        while c<=e and c in exclude do
            c := c+1;
        od;

        return c;
    end,

    _simple_coloring := meth(self, pre, post, taken)
        local q, x, y, ta, locked, res, id_set, id_cnt, col, parb, tas_on_altb;

        id_set := UnionSet(UserRecFields(pre), UserRecFields(post));
        id_cnt := Length(id_set);

        x := List(Filtered(id_set, v -> not IsBound(pre.(v))),  vid -> var.table.(vid) );
        y := List(Filtered(id_set, v -> not IsBound(post.(v))), vid -> var.table.(vid) );

        res := Copy(taken);
        locked := Set(List(UserRecFields(taken), v -> taken.(v)));
        q := Filtered( Flat(List(x, v -> post.(v.id))), v -> not v in y and not IsBound(taken.(v.id)) );

        while q <> [] do
            [ta, q] := SplitAt(q, 1);
            ta := ta[1];
            col := self._next_color(1, id_cnt, locked
                                                ::List(Filtered(pre.(ta.id), t -> IsBound(res.(t.id)) ), v -> res.(v.id))
                        );
            res.(ta.id) := col;

            #Setup next iteration
            if IsBound(post.(ta.id)) and ForAll(Flat( List(post.(ta.id), p -> pre.(p.id) ) ),
                                            v -> IsBound(res.(v.id))) then
                    # push on queue
                    q := Filtered(post.(ta.id), p -> not p in y)::q;
            fi;
            # else do nothing - still need to complete other par. branches

            if Length(pre.(ta.id)) > 1 then
                #Unlock color of previous divergent point
                locked := DropLast(locked, 1);
            fi;
            if IsBound(post.(ta.id)) and Length(post.(ta.id)) > 1 then
                Add(locked, col);
            fi;

        od;

        return res;
    end,

    _map_col2mem := function(coloring, decl_specs)
        local c, p, ta, cols, samecol_tas, ta2mem, mem_list;

        cols := Set(Flat(List(UserRecFields(coloring), v -> coloring.(v))));
        ta2mem := rec();
        mem_list := [];
        for c in cols do
            samecol_tas := List(Filtered(UserRecFields(coloring), v -> c = coloring.(v)), vid -> var.table.(vid));
            p := var.fresh_t("P", TArray(samecol_tas[1].t.t, Maximum(List(samecol_tas, v -> v.t.size))) );
            p.("decl_specs") := decl_specs;
            for ta in samecol_tas do
                ta2mem.(ta.id) := p;
            od;
            Add(mem_list, p);
        od;

        return [ta2mem, mem_list];
    end,

    compact_arrays := meth(self, ker, ta_list, ker_args, decl_specs, opts)
        local cross_ker_tas, pre, post, ck_pre, ck_post, ck_col, coloring, ta2mem;

        #Currently limited to coloring cross-kernel tas
        cross_ker_tas := Filtered(ta_list, v -> IsBound(v.crossker) and v.crossker);
        ta_list := Difference(ta_list, cross_ker_tas);
        [ck_pre, ck_post] := self._reduce_tarray_graph(self.tarray_pre, self.tarray_post, cross_ker_tas::ker_args);
        ck_col := self._simple_coloring(ck_pre, ck_post, rec());
        #assign mem areas based on coloring
        [ta2mem, cross_ker_tas] := self._map_col2mem(ck_col, decl_specs);
        ker := SubstVars(ker, ta2mem);

        return [ker, cross_ker_tas::ta_list];
    end,

    extract_tarrays := meth(self, ker, cond, opts)
        local decls, d, ta, tarrays;

        decls := Collect(ker, @(1, decl, d -> ForAny(d.vars, v -> IsArrayT(v.t)) ) );

        tarrays := [];
        for d in decls do
            # Collect temp arrays based on cond
            ta := Filtered(d.vars, v -> cond(v) );
            d.vars := Difference(d.vars, ta);

            Append(tarrays, ta);
        od;

        return tarrays;
    end,

    to_grid_arrays := meth(self, ker, ta_list, ker_args, opts)
        local ta, idx_list, blk_ids, th_ids, idx, idx_cnt, idx_rng, derefs, nths, v;

        for ta in ta_list do
            ta.("decl_specs") := ["__device__"];
            idx_list := When(IsBound(ta.parctx), ta.parctx, []);
            #parctx format: [ [loopvar1, simtIdx1], ... ]
            blk_ids := Filtered(idx_list, idx -> IsSimtBlockIdx(idx[2])); Sort(blk_ids, (i1, i2) -> i1[2].name < i2[2].name);
            th_ids  := Filtered(idx_list, idx -> IsSimtThreadIdx(idx[2])); Sort(th_ids, (i1, i2) -> i1[2].name < i2[2].name);
            for idx in th_ids::blk_ids do
                idx_cnt := idx[2].count(); # Num of threads associated with index
                if idx_cnt > 1 then
                    idx_rng := idx[2].get_rng(); # Thread range: [<1st tid>, <last tid>]
                    derefs := Collect(ker, @(1, deref, v -> ta in v.free()));
                    nths   := Collect(ker, @(1,   nth, v -> ta in v.free()));

                    for v in derefs do
                        v.loc := add(v.loc, ta.t.size*(idx[2]-idx_rng[1]) );
                    od;
                    for v in nths do
                        v.idx := add(v.idx, ta.t.size*(idx[2]-idx_rng[1]));
                    od;
                    ta.t.size := ta.t.size*idx_cnt;
                fi;
            od;
        od;

        if opts.mempool then
            [ker, ta_list] := self.compact_arrays(ker, ta_list, ker_args, ["__device__"], opts);
        fi;

        return [ker, ta_list];
    end,

    to_block_shmem_arrays := meth(self, ker, ta_list, ker_args, opts)
        local ta, idx_list, blk_ids, th_ids, idx, idx_cnt, idx_rng, derefs, nths, size, padded, v;

        for ta in ta_list do
            padded := false;
            ta.("decl_specs") := ["__shared__"];
            idx_list := When(IsBound(ta.parctx), ta.parctx, []);
            th_ids  := Filtered(idx_list, idx -> IsSimtThreadIdx(idx[2])); Sort(th_ids, (i1, i2) -> i1[2].name < i2[2].name);
            for idx in th_ids do
                idx_cnt := idx[2].count();
                if idx_cnt > 1 then
                    size := ta.t.size;
                    if not padded and Mod(size, 2) = 0 then
                        size := size+1;
                    fi;
                    padded := true;
                    idx_rng := idx[2].get_rng();
                    derefs := Collect(ker, @(1, deref, v -> ta in v.free()));
                    nths   := Collect(ker, @(1,   nth, v -> ta in v.free()));

                    for v in derefs do
                        v.loc := add(v.loc, size*(idx[2]-idx_rng[1]) );
                    od;
                    for v in nths do
                        v.idx := add(v.idx, size*(idx[2]-idx_rng[1]));
                    od;
                    ta.t.size := size*idx_cnt;
                fi;
            od;
        od;

        ker := decl(ta_list, ker);

        return [ker, ta_list];
    end,

    extract_ker_datas := meth(self, cuda_icode, opts)
        local ker_datas, tcuicode, d;

        ker_datas := List(Collect(cuda_icode, data), d -> d.var);

        tcuicode := SubstTopDown(cuda_icode, data, e -> e.cmd);
        while tcuicode <> cuda_icode do
            cuda_icode := tcuicode;
            tcuicode := SubstTopDown(cuda_icode, data, e -> e.cmd);
        od;

        return [cuda_icode, ker_datas];
    end,

    _get_data_loopindex_deps := function(icode, d, opts)
        local ln, fv, la, li, lv, locli;
        ln := Collect(icode, @(1, nth).cond(v -> d in v.loc.free()));
        fv := Set( Flat(List(ln, v -> v.idx.free())) );
        li := Filtered(fv, v -> IsLoopIndex(v));

        #If code in SSA form using lv set to avoid cyclic deps could be avoided
        fv := Difference(fv, li); lv := Copy(fv);
        while(fv <> []) do
            la := Collect(icode, @(1, assign).cond(a -> a.loc in fv));
            fv := Set( Flat(List(la, a -> a.exp.free())) );
            locli := Filtered(fv, v -> IsLoopIndex(v));
            fv := Difference(fv, locli);
            fv := Difference(fv, lv);
            li := UnionSet(li, locli);
            lv := UnionSet(lv, fv);
        od;

        return li;
    end,

    ker_datas_to_device := meth(self, ker_datas, cuda_icode, opts)
        local d, li;

        for d in Reversed(ker_datas) do
            li := self._get_data_loopindex_deps(cuda_icode, d, opts);
            # If any warp-divergent access then allocate in dev mem otherwise constant.
            if ForAny(li, v -> IsBound(self.simtidx_rec.(v.id))
                                and ObjId(self.simtidx_rec.(v.id)) = simtThreadIdxX) then
                d.("decl_specs") := ["__device__"];
            else
                d.("decl_specs") := ["__constant__"];
            fi;
            cuda_icode := data(d, d.value, cuda_icode);
        od;

        return cuda_icode;
    end,

    make_kernels := meth(self, full_kernel, device_data, full_args, opts)
        local kercx, ker_cnt, ker, ker_args, ker_datas, cuda_ker, cuda_desc, cuda_icode,
                dim_grid, dim_block, tbody, tas_grid, tas_block_shmem, cross_ker_tas, ta;

        cuda_desc := rec(grid_tarrays := [], cuda_kers := [] );

        [full_kernel, ker_datas] := self.extract_ker_datas(full_kernel, opts);

        # Mark cross kernel temporaries
        cross_ker_tas := Filtered(Flat(List(Collect(full_kernel, @@(1, decl, (d, ctx) -> not IsBound(ctx.simt_block) or ctx.simt_block = [] )),
                            d -> d.vars)), v -> IsArrayT(v.t));
        DoForAll(cross_ker_tas, v -> v.setAttrTo("crossker", true));
        #Extract all temp arrays that should map to dev mem
        tas_grid := self.extract_tarrays(full_kernel,
                                            When(opts.use_shmem,
                                                    v -> (IsBound(v.crossker) and v.crossker) or (IsBound(v.t.size) and (v.t.size >= opts.max_shmem_size/8/2)),
                                                    True
                                            ), opts);
        #Allocate in dev mem
        [full_kernel, tas_grid] := self.to_grid_arrays(full_kernel, tas_grid, full_args, opts);
        Append(cuda_desc.grid_tarrays, tas_grid);

        kercx := Collect(full_kernel, simt_block);
        kercx := When(kercx = [], [ full_kernel ], kercx);
        ker_cnt := 0;
        for ker in kercx do

            tbody := When(ObjId(ker) = simt_block, chain(ker.cmds), ker);
            ker_args := Intersection(device_data, tbody.free())::Intersection(full_args, tbody.free());
            #Extract remaining shared mem temp arrays if required
            tas_block_shmem := self.extract_tarrays(tbody, v -> opts.use_shmem, opts);
            #Allocate in shared mem
            [tbody, tas_block_shmem] := self.to_block_shmem_arrays(tbody, tas_block_shmem, ker_args, opts);

            dim_grid := var.fresh_t_obj("g", Dim3(), rec( x := self.get_dim_idx(tbody, simtBlockIdxX),
                                                          y := self.get_dim_idx(tbody, simtBlockIdxY),
                                                          z := self.get_dim_idx(tbody, simtBlockIdxZ)
                                                        ) );
            dim_block := var.fresh_t_obj("b", Dim3(), rec( x := self.get_dim_idx(tbody, simtThreadIdxX),
                                                           y := self.get_dim_idx(tbody, simtThreadIdxY),
                                                           z := self.get_dim_idx(tbody, simtThreadIdxZ)
                                                        ) );

            cuda_ker := specifiers_func(["__global__"], TVoid, "ker_code"::StringInt(ker_cnt), Filtered(ker_args, d -> not IsBound(d.("decl_specs")) or not "__constant__" in d.("decl_specs")), tbody );

            Add(cuda_desc.cuda_kers,  rec(dim_grid := dim_grid, dim_block := dim_block, cuda_ker := cuda_ker));

            ker_cnt := ker_cnt + 1;
        od;

        cuda_icode := chain(List(cuda_desc.cuda_kers, ck -> ck.cuda_ker));
        cuda_icode := self.ker_datas_to_device(ker_datas, cuda_icode, opts);

        return [cuda_desc, cuda_icode];

    end,

    make_init := meth(self, initname, initparams, init_datas, device_datas, cuda_desc, io, params, o_dims, opts)
        local d, dd, t_id, t_var, tarrays, init_chain, init_func;

        # Use the generated init for the host arrays
        for d in init_datas do
            t_var := d.var;
            d.var := t_var.clone();
            t_id := d.var.id;
            d.var.id := "h"::t_var.id;
            d.init := SubstVars(d.init, rec((t_var.id) := d.var));
            Unbind(var.table.(t_id));
            var.table.(d.var.id) := d.var;
        od;

        init_chain := chain(List(device_datas, dd ->    When(IsBound(dd.("decl_specs")) and "__constant__" in dd.("decl_specs"),
                                                            chain(SReduce(dd.init, opts),
                                                                cu_memcpy_to_sym(dd, var.table.("h"::dd.id), dd.t.size)
                                                                ),
                                                            chain(cu_allocate(dd, dd.t.t, dd.t.size),
                                                                SReduce(dd.init, opts),
                                                                cu_memcpy(dd, var.table.("h"::dd.id), dd.t.size, cudaMemcpyHostToDevice)
                                                                )
                                                        )
                                )
                                ::When(not (IsBound(opts.devFunc) and opts.devFunc), List(params, dd -> cu_allocate(dd, dd.t.t, dd.t.size) ), [])
                            );
        if not (IsBound(opts.devFunc) and opts.devFunc) then
            for i in [1..Length(io)] do
                Add(init_chain.cmds, cu_allocate(io[i], io[i].t.t, o_dims[i]));
            od;
        fi;

        init_func := func(TVoid, initname, Set(Collect(init_chain, param)), init_chain);
        return init_func;
    end,

    make_destroy := meth(self, destroyname, init_datas, device_datas, cuda_desc, io, params, opts)
        local destroy_chain, destroy_func;

        destroy_chain := chain( List(Filtered(device_datas, dd -> not IsBound(dd.("decl_specs")) or not "__constant__" in dd.("decl_specs") ),
                                    d -> cu_free(d))
                                ::When(not (IsBound(opts.devFunc) and opts.devFunc), List(params, d -> cu_free(d)), [])
                                ::When(IsBound(opts.devFunc) and opts.devFunc, [], List(io, d -> cu_free(d))));

        destroy_func := func(TVoid, destroyname, Set(Collect(destroy_chain, param)), destroy_chain);
        return destroy_func;
    end,

    assemble_prog := meth(self, datas, initparams, params, io, o, icode, opts)
        local hparams, hio, p, hp, init_datas, device_datas, ker_args,
              sub, cuda_sub, initsub, destroysub,
              cuda_desc, cuda_icode, icode_wrap, initcode, destroycode, prog,
              wrapk_f;

        wrapk_f := When(IsBound(opts.devFunc) and opts.devFunc, self.wrap_kernel_devfunc, self.wrap_kernel);
        hparams := [];
        for p in params do
            p.setAttr("decl_as_ptr");
            hp := p.clone();
            Unbind(var.table.(hp.id));
            hp.id := When(IsBound(opts.devFunc) and opts.devFunc, p.id, "h"::p.id);
            var.table.(hp.id) := hp;
            Add(hparams, hp);
        od;

        hio := [];

        if not (IsBound(opts.devFunc) and opts.devFunc) then
            for p in io do
                hp := p.clone();
                Unbind(var.table.(hp.id));
                hp.id := "h"::p.id;
                var.table.(hp.id) := hp;
                Add(hio, hp);
            od;
        fi;

        # Device data here refers to data initialized in init function and therefor residing in dev mem
        init_datas := Filtered(datas, e -> IsBound(e.var.init));
        device_datas := List(init_datas, d -> d.var);

        sub := Cond(IsBound(opts.subName), opts.subName, "transform");
        cuda_sub := "transform_ker";
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        destroysub := Cond(IsBound(opts.subName), Concat("destroy_", opts.subName), "destroy");
        ker_args := Concatenation(io, params);

        [cuda_desc, cuda_icode] := self.make_kernels(icode, device_datas, ker_args, opts);

        icode_wrap := func(TVoid, sub, Concatenation(When(IsBound(opts.devFunc) and opts.devFunc, io, hio), hparams), chain(wrapk_f(self, cuda_desc, io, params, o.dimensions, opts)));
        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            initcode := self.make_init(initsub, initparams, init_datas, device_datas, cuda_desc, io, params, o.dimensions, opts);
            destroycode := self.make_destroy(destroysub, init_datas, device_datas, cuda_desc, io, params, opts);
            prog := program(
                decl(List(datas, x->x.var)::cuda_desc.grid_tarrays::device_datas::params::When(IsBound(opts.devFunc) and opts.devFunc, [], io),
                    chain(
                        initcode,
                        cuda_icode,
                        icode_wrap,
                        destroycode
                    )));
        else
            initcode := self.make_init(initsub, initparams, init_datas, device_datas, cuda_desc, io, params, o.dimensions, opts);
            destroycode := self.make_destroy(destroysub, init_datas, device_datas, cuda_desc, io, params, opts);
            prog := program(
                decl(cuda_desc.grid_tarrays::params::When(IsBound(opts.devFunc) and opts.devFunc, [], io),
                    chain(
                        initcode,
                        cuda_icode,
                        icode_wrap,
                        destroycode
                    )));
        fi;

        return prog;
    end,

    Formula := meth(self, o, y, x, opts)
        local icode, icode1, datas, prog, initparams, params, io, t, o, field;

        o := SumsUnification(o.child(1), opts);

        [x, y] := self.initXY(x, y, opts);

        io := When(x=y, [x], [y, x]);
        io := io::When(IsBound(opts.Xptr), [opts.Xptr], []);
        io := io::When(IsBound(opts.Yptr), [opts.Yptr], []);

        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        if opts.mempool then
            self.tarray_pre  := [rec()];
            self.tarray_post := [rec()];
        fi;
        self.simtidx_rec := rec();

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);

        if opts.mempool then
            self.tarray_pre := self.tarray_pre[1];
            self.tarray_post := self.tarray_post[1];
            for field in Difference(UserRecFields(self.tarray_pre), List(io, v -> v.id)) do
                self.tarray_pre.(field)  := Set(Flat(self.tarray_pre.(field)));
                self.tarray_post.(field) := Set(Flat(self.tarray_post.(field)));
            od;
            for field in List(io, v -> v.id) do
                if IsBound(self.tarray_pre.(field)) then
                    self.tarray_pre.(field)  := Set(Flat(self.tarray_pre.(field)));
                fi;
                if IsBound(self.tarray_post.(field)) then
                    self.tarray_post.(field) := Set(Flat(self.tarray_post.(field)));
                fi;
            od;

        fi;

        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        if IsBound(opts.postCompileStrategy) then
            icode := opts.postCompileStrategy(icode, opts);
        fi;

        initparams := Copy(params);
        if IsBound(opts.symbol) then
          params := Concatenation(params, opts.symbol);
        fi;

        if IsBound(opts.accStrategy) then
          icode.iy := y;
          icode.iy.n := o.dims()[1];
          icode.ix := x;
          icode.ix.n := o.dims()[2];
          icode.ivars := Concatenation(params, List(datas, x->x.var));
          icode := opts.accStrategy(icode);
        fi;

        prog := self.assemble_prog(datas, initparams, params, io, o, icode, opts);
        prog.dimensions := o.dims();
        return prog;
    end,
    ));

max_stars_repo_name: broderickpt/spiral-package-simt

max_issues_repo_name: broderickpt/spiral-package-simt

max_forks_repo_name: broderickpt/spiral-package-simt


text:
##  Copyright (c) 2018-2021, Carnegie Mellon University
##  See LICENSE for details

# FFTX Baseline options

_noT := r -> CopyFields(r, rec(forTransposition := false));

Declare(ParseOpts);

Class(FFTXDefaultConf, LocalConfig, rec(
    getOpts := (self, t) >> ParseOpts(self, t),
    operations := rec(Print := s -> Print("<FFTX Default Configuration>")),
));

Class(FFTXOpts, SpiralDefaults, FFTXGenMixin, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX options record>")),
    breakdownRules := Copy(SpiralDefaults.breakdownRules),
    codegen := Copy(SpiralDefaults.codegen),
    tagIt := (self, t) >> t.withTags(self.tags)
));

FFTXOpts.codegen.GathPtr := fftx.codegen.MultiPtrCodegenMixin.GathPtr;
FFTXOpts.codegen.ScatPtr := fftx.codegen.MultiPtrCodegenMixin.ScatPtr;
FFTXOpts.codegen.OO := (self, o, y, x, opts) >> skip();
FFTXOpts.sumsgen.IterHStack := MultiPtrSumsgenMixin.IterHStack;

Class(FFTXConvOpts, FFTXOpts, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX Conv options record>")),
    breakdownRules := Copy(SpiralDefaults.breakdownRules),
    codegen := CopyFields(MultiPtrCodegenMixin, SpiralDefaults.codegen)
));

Class(FFTXGlobals, rec(
    registerConf := meth(arg)
        local self, key;
        self := arg[1];
        self.(arg[2].defaultName) := arg[2].defaultOpts;
        key := Filtered(RecFields(arg[2].defaultOpts()), i -> not IsSystemRecField(i));
        self.supportedConfs.(key[1]) := arg[2].confHandler;
    end,
    supportedConfs := rec(
        mdRConv :=
            function(arg)
                local opts;
                opts := Copy(FFTXConvOpts);
                opts.useDeref := false;
                opts.breakdownRules.PRDFT := List([ PRDFT1_Base2, PRDFT1_CT, PRDFT_PD], _noT);
                opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD], _noT);
                opts.breakdownRules.IPRDFT2 := List([ IPRDFT2_Base1, IPRDFT2_Base2, IPRDFT2_CT], _noT);
                opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT, PRDFT3_OddToPRDFT1], _noT);
                opts.breakdownRules.URDFT := List([ URDFT1_Base1, URDFT1_Base2, URDFT1_Base4, URDFT1_CT ], _noT);
                opts.breakdownRules.DFT := [DFT_Base, DFT_PD, DFT_CT];
                opts.breakdownRules.IOPrunedMDRConv := [IOPrunedMDRConv_3D_2trip_zyx_freqdata];
                opts.breakdownRules.MDDFT := [ MDDFT_Base, MDDFT_RowCol ];
                opts.breakdownRules.GT := [ GT_NthLoop ];
                opts.breakdownRules.DFT := [ DFT_Base, DFT_CT, DFT_Rader, CopyFields(DFT_GoodThomas, rec(maxSize := 15)), DFT_PD ];
                opts.topTransforms := [DFT, MDDFT, PrunedDFT, PrunedIDFT, IOPrunedConv, PrunedPRDFT, PrunedIPRDFT, TTensorI, IOPrunedMDRConv,
                    PRDFT, IPRDFT, PRDFT1, IPRDFT1, PRDFT2, IPRDFT2, PRDFT3, IPRDFT3, TTensorI ];
                opts.breakdownRules.MDRConv := [ MDRConv_3D_2trip_zyx_freqdata ]; # [MDRConv_Base]
                return opts;
            end,
    ),
    defaultConf := (arg) >> FFTXDefaultConf,
    mdRConv := (arg) >> rec(mdRConv := true),
    getOpts :=
        meth(arg)
            local opts, rfields;
            if Length(arg) >= 2 then
                rfields := Filtered(RecFields(arg[2]), i -> not IsSystemRecField(i));
                if Length(rfields) >= 1 and IsBound(arg[1].supportedConfs.(rfields[1])) then
                    return ApplyFunc(arg[1].supportedConfs.(rfields[1]), Drop(arg, 1));
                fi;
            fi;
            opts := Copy(FFTXOpts);
            opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
            opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_PD, PRDFT_Rader], _noT);
            opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD, IPRDFT_Rader ], _noT);
            opts.breakdownRules.IPRDFT2 := List([ IPRDFT2_Base1, IPRDFT2_Base2, IPRDFT2_CT], _noT);
            opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);
            opts.breakdownRules.MDDFT := [ MDDFT_Base, MDDFT_RowCol ];
            opts.breakdownRules.DFT := [ DFT_Base, DFT_CT, DFT_Rader, CopyFields(DFT_GoodThomas, rec(maxSize := 15)), DFT_PD ];
            opts.breakdownRules.GT := [ GT_NthLoop ];
            opts.breakdownRules.MDRConv := [ MDRConv_3D_2trip_zyx_freqdata ]; # [MDRConv_Base]
            return opts;
        end
));

_promote1 := t >> let(t1 := RulesFFTXPromoteNT(Copy(t)), RulesFFTXPromoteNT_Cleanup(t1));

# this is a first experimental opts-deriving logic. This needs to be done extensible and properly
ParseOpts := function(conf, t)
    local tt, _tt, _conf, _opts;

    if conf = FFTXDefaultConf then
        # promote with default conf rules
        tt := _promote1(Copy(t));

        if ObjId(tt) = TFCall then
            _tt := tt.params[1];
            # check for convolution
            if (ObjId(_tt) = MDRConv) or ((ObjId(_tt) = TTensorI) and (ObjId(_tt.params[1]) = MDRConv)) then
                _conf := FFTXGlobals.mdRConv();
                _opts := FFTXGlobals.getOpts(_conf);
                return _opts;
            fi;
            # check for Hockney. This is for N=128
            if ObjId(_tt) = IOPrunedMDRConv then
                _conf := FFTXGlobals.defaultHockneyConf(rec(globalUnrolling := 16, prunedBasemaxSize := 7));
                _opts := FFTXGlobals.getOpts(_conf);
                return _opts;
            fi;
        fi;

        # check for WarpX
        _conf := FFTXGlobals.defaultWarpXConf();
        _opts := FFTXGlobals.getOpts(_conf);
        tt := _opts.preProcess(Copy(t));
        if ObjId(tt) = TFCall and ObjId(tt.params[1]) = TCompose then
            _tt := tt.params[1].params[1];
            # detect promoted WarpX
            if IsList(_tt) and Length(_tt) = 3 and List(_tt, ObjId) = [ TNoDiagPullinRight, TRC, TNoDiagPullinLeft ] then
                return _opts;
            fi;
        fi;
        # we are doing nothing special
        return FFTXGlobals.getOpts();
    fi;

    # Here we have to handle GPU configs
    Error("Don't know how to derive opts!\n");
end;

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx