text:
Class(NewScalarProduct, BaseOperation, rec(
  __call__ := (self, rv) >> SPL(WithBases(self, rec(
    rv := rv,
    _children := [rv]
  ))),
  dims := self >> [self.rv.element.domain(), V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.rv, ")")
));

Class(DenseScalarProduct, BaseOperation, rec(
  __call__ := (self, rv) >> SPL(WithBases(self, rec(
    rv := rv,
    _children := [rv]
  ))),
  dims := self >> [self.rv.element.domain(), V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.rv, ")")
));

Class(SparseScalarProduct, BaseOperation, rec(
  __call__ := (self, rv) >> SPL(WithBases(self, rec(
    rv := rv,
    _children := [rv]
  ))),
  dims := self >> [self.rv.element.domain(), V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.rv, ")")
));

Declare(IsSparseT);

Class(sparse_nth3, Loc, rec(
    __call__ := (self, list, idx) >> WithBases(self,
        rec(operations := NthOps,
			list := list,
            idx := toExpArg(idx))).setType().cfold(),

    can_fold := self >> self.idx _is funcExp or (IsValue(self.idx) and
                  (IsValue(self.list) or (IsVar(self.list) and IsBound(self.list.var)) or self.list _is apack)),
    cfold := self >> When(self.can_fold(), self.eval(), self),

    rChildren := self >> [self.list, self.idx],
    rSetChild := rSetChildFields("list", "idx"),

	eval := meth(self)
		local result, iterator;
		result := 0;
		iterator := 1;
		while iterator <= Length(self.list.var) and result = 0 do
			if self.idx = self.list.var[iterator].pair.first then
				#Print("if\n");
				result := self.list.var[iterator].pair.second;
				iterator := iterator + 1;
				break();
				#Print(iterator);
			else
				#Print("else\n");
				result := V(0);
				iterator := iterator + 1;
				#Print(iterator);
			fi;
		od;
		return result;
		end,

	#	result := var.fresh_t("result", TInt);
	#	g := var.fresh_t("g", TInt);
	#	decl[result, g], chain(
	#	loopf(g, V(1), Length(fds.var),
	#	chain(IF(eq(x, fds.var[g].pair.first),
	#	assign(result,fds.var[g].pair.second),
	#	assign(result,V(0))))),
	#	return result);
	#  end;

    computeType := self >> Cond(
	IsArrayT(self.list.range()), self.list.range(),
        ObjId(self.list.t) = TSym, TSym("Containee"), #used with C++ container objects (EnvList)
        self.list.t = TUnknown,  self.list.range(),
	Error("Unknown types of 1st argument <self.list> in ", ObjId(self))
    ),
));

Class(TStruct, CompositeTyp, rec(
	__call__ := arg >> let(
	self := arg[1],
	fields := CopyFields(arg[2]),
	WithBases(self,
	rec(fields := fields, operations := TypOps))),
    print := self >> Print(self.name,"(", self.fields.name, ", ", self.fields.key, ", ", self.fields.val, ")"),
));

Class(TPair, TStruct, rec(
   __call__ := arg >> let(
   self := arg[1],
   pair := CopyFields(arg[2]),
   WithBases(self,
   rec(pair := pair, operations := TypOps))),

   print := self >> Print(self.name,"(",self.pair.first,", ", self.pair.second,")"),
));

Class(KVPair, CompositeTyp, rec(
	__call__ := (self, first ,second) >> WithBases(self,rec(
		first := first,
		second := second,
		operations := TypOps)),
		#t := Checked(IsType(t), t),
		first := self >> self.first,
		second := self >> self.second,
	print := self >> Print(self.name,"(",self.first,", ",self.second,")"),
));

SparseMatSPL := function ( S )
  local M, t, v, i;
  t := 0;
  i := 1;
  if IsSparseT(S.data_type) then
	M := NullMat(S.data_type.dims()[1], S.data_type.dims()[2]);
	while t < S.data_type.size()-1 do
		Print(t);
		if t = S.element[i][1] then
			M[1][t+1] := S.element[i][2];
			t := t + 1;
			i := i + 1;
		else
			M[1][t+1] := 0;
			t := t + 1;
		fi;
  	od;
  fi;
  return M;
  end;

#F sparse_nth(<loc>, <idx>) -- symbolic representation of array access
#F
Class(sparse_nth, Loc, rec(
    __call__ := (self, loc, idx) >> WithBases(self,
        rec(operations := NthOps,
            loc := toExpArg(loc),
            idx := toExpArg(idx))).setType().cfold(),

    can_fold := self >> self.idx _is funcExp or (IsValue(self.idx) and
                  (IsValue(self.loc) or (IsVar(self.loc) and IsBound(self.loc.value)) or self.loc _is apack)),
    cfold := self >> When(self.can_fold(), self.eval(), self),

    rChildren := self >> [self.loc, self.idx],
    rSetChild := rSetChildFields("loc", "idx"),

    ev := self >> let(e := self.eval(),
	Cond(IsBound(e.v), e.v, e)),

    eval := self >> let(loc := self.loc.eval(), idx := self.idx.eval(),
        Cond(IsList(loc),
				 Cond(idx.v >= Length(loc), errExp(self.t), Cond(idx.v = loc[1][x][1], V(0), V(loc))),
			 idx _is funcExp,
                 self.t.value(idx.args[1]),
             not IsValue(idx),
                 self,
             idx.v < 0,
                 errExp(self.t),
             loc _is apack,
                 Cond(idx.v >= Length(loc.args), errExp(self.t), loc.args[idx.v+1]),
             IsValue(loc),
                 Cond(idx.v >= Length(loc.v), errExp(self.t), V(loc.v[idx.v+1])),
             IsVar(loc) and IsBound(loc.value),
                 Cond(idx.v >= Length(loc.value.v), errExp(self.t), V(loc.value.v[idx.v+1])),
             self)),

    computeType := self >> Cond(
	IsPtrT(self.loc.t) or IsArrayT(self.loc.t) or IsListT(self.loc.t), self.loc.t.t,
        ObjId(self.loc.t)=TSym, TSym("Containee"), #used with C++ container objects (EnvList)
        self.loc.t = TUnknown,  self.loc.t,
	Error("Unknown types of 1st argument <self.loc> in ", ObjId(self))
    ),

    isExpComposite := true
));

#IsSparseOfs := o -> IsRec(o) and IsBound(o.isSparseOfs) and o.isSparseOfs;

#F FDataOfs(<datavar>, <len>, <ofs>)
#
Class(FDataSparseOfs, Function, rec(
    #__call__ := (self, datavar, len, ofs) >> WithBases(self, rec(
    #var := datavar,
    #operations := PrintOps,
    #ofs := toExpArg(ofs),
    #len := Checked(IsPosIntSym(len), len)
    #)),
   __call__ := arg >> let(
       self := arg[1],
       object := arg[2],
	   len := arg[3],
       ofs := toExpArg(arg[4]),
	   WithBases(self, rec(var := object, len := len, ofs := ofs, operations := PrintOps))),

   print := self >> Print(self.name, "(", self.var, ", ", self.len, ", ", self.ofs, ")"),
   rChildren := self >> [ self.var, self.len, self.ofs ], #DOMAIN AND RANGE BROKEN
   rSetChild := rSetChildFields("var", "len", "ofs" ),

# <-Daniele's changes
#    rChildren := self >> [ self.var, self.len, self.ofs],
#    rSetChild := rSetChildFields("var", "len", "ofs"),

#    rChildren := self >> [ self.var, self.len, self.ofs, self._domain, self._range],
#    rSetChild := rSetChildFields("var", "len", "ofs", "_domain", "_range"),
#    from_rChildren := (self, rch) >> ObjId(self)(rch[1], rch[2], rch[3]).setDomain(rch[4]).setRange(rch[5]),
# ->

    #domain := self >> self.len,

    #at := (self, n) >> When(IsInt(n) and IsValue(self.ofs) and IsBound(self.var.value),
    #   self.var.value.v[n + self.ofs.v + 1],
    #    nth(self.var, n + self.ofs)),
	#
    #tolist := self >> List([0..EvalScalar(self.len-1)], i -> nth(self.var, self.ofs+i)),
    #lambda := self >> let(x := Ind(self.domain()), Lambda(x, nth(self.var, self.ofs+x))),

	at := (self, n) >> When(IsInt(n) and IsValue(self.ofs) and IsBound(self.var.value),
       self.var.value.v[n + self.ofs.v + 1],
        sparse_nth3(self.var, n + self.ofs)),

    tolist := self >> List([0..EvalScalar(self.len-1)], i -> sparse_nth3(self.var, self.ofs+i)),
    lambda := self >> let(x := Ind(self.domain()), Lambda(x, sparse_nth3(self.var, self.ofs+x))),

	range := self >> When(self._range=false, self.var.t, self._range),
    domain := self >> self.len,
    #range := self >> When(self._range=false, self.var.t.t, self._range),
    inline := true,
    free := self >> self.ofs.free()
));

#F FDataSparse(<datavar>) -- symbolic function i -> datavar[i],
#F
#F domain = datavar.range
#F range = datavar.t
#F
#F
Class(FDataSparse, Function, rec(
   __call__ := arg >> let(
       self := arg[1],
       object := arg[2],
       var := Cond(IsList(arg[3]), arg[3], [arg[3]]),
	   WithBases(self, rec(var := var, object := object, operations := PrintOps))),

   print := self >> Print(self.name, "(", self.object, ", ", self.var, ")"),
   rChildren := self >> [ self.var ],
   rSetChild := rSetChildFields("var"),

   at := (self, n) >> self.lambda().at(n),
   tolist := self >> self.lambda().tolist(),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, sparse_nth3(self, x))),

   domain := self >> self.object.size,
   range := self >> self.object.t,

   inline := true,
   free := self >> Set([]),
   part := (self, len, ofs) >> FDataSparseOfs(self.var, len, ofs),
));

Class(FDataSparseMat, Function, rec(
	__call__ := arg >> let(
       self := arg[1],
       object := arg[2],
       var := Cond(IsList(arg[3]), arg[3], [arg[3]]),
	   WithBases(self, rec(var := var, object := object, operations := PrintOps))),

	print := self >> Print(self.name, "(", self.object, ", ", self.var, ")"),
	rChildren := self >> [self.var],
	rSetChild := rSetChildFields("var"),

	at := (self, n) >> self.lambda().at(n),
	tolist := self >> self.lambda().tolist(),
	#lambda := self >> let(k := Ind(object.expr.loc.loc.dims[1]))
	lambda := self >> let(k := Ind(self.domain()), u := Ind(self.domain()), Lambda(k, Lambda(u, self.object.at(k,u)))),

	domain := self >> Length(self.var),
	range := self >> self.object.t.params[1],
));

Class(FDataSparseMatOfs, Function, rec(
	__call__ := arg >> let(
       self := arg[1],
       object := arg[2],
	   len := arg[3],
       ofs := arg[4],
	   WithBases(self, rec(object := object, len := len, ofs := ofs, operations := PrintOps))),

	print := self >> Print(self.name, "(", self.object, ", ", self.len, ", ", self.ofs, ")"),
   	rChildren := self >> [ self.object, self.len, self.ofs ],
   	rSetChild := rSetChildFields("object", "len", "ofs"),

	at := (self, n) >> When(IsInt(n) and IsValue(self.ofs) and IsBound(self.object.value),
       self.object.value.v[n + self.ofs.v + 1],
        nth(self.object, self.lambda().at(n))),

	#tolist := self >> self.lambda().tolist(),
	tolist := self >> List([0..EvalScalar(self.len-1)], i -> nth(self.object, self.ofs+i)),
	lambda := self >> let(k := Ind(self.domain()), u := Ind(self.domain()), Lambda(k, Lambda(u, self.object.at(k,u)))),

	range := self >> self.object.t,
    domain := self >> self.len,
	inline := true,
    free := self >> self.ofs.free()
));

Declare(TSparse);
Declare(Ttrace);

Class(sparse_nth2, BaseOperation, rec(
	__call__ := (self, sa, index) >> SPL(WithBases(self, rec(
	  sa := sa,
	  index := index,
	  _children := [index]
	  ))),
	  dims := self >> [self.index, V(1)],
	  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
	  doNotMarkBB := true,
));

#Class(GathPtr, Gath, rec(
#    rChildren := self >> [self.ptr, self.func],
#    rSetChild := rSetChildFields("ptr", "func"),
#    new := (self, ptr, func) >> SPL(WithBases(self, rec(
#        ptr := ptr,
#      	func := Checked(IsFunction(func) or IsFuncExp(func), func)))).setDims()
#));

Class(TSparse2, SumsBase, BaseMat, rec(
    _short_print := true,
	rChildren := self >> [self.t, self.size1, self.list1, self.size2, self.list2],
	rSetChild := rSetChildFields("t","size1", "list1", "size2", "list2"),
	dims := self >> [self.size1, self.size2],

	new := (self, t, size1, list1, size2, list2) >> SPL(WithBases(self,
		rec(dimensions := [size1, size2],
		size1 := size1,
		size2 := size2,
		t := t,
		list1 := list1,
		list2 := list2))),

));

Class(struct_nth, nth, rec(
	__call__ := (self, loc, elem, idx) >> WithBases(self, rec(
		loc := loc,
		elem := elem,
		idx := idx,
		t := Cond(IsPtrT(loc) or IsSparseT(loc.t), loc.t.t, loc.t),
		operations := NthOps
	)),
	rChildren := self >> [self.loc, self.elem, self.idx],
	rSetChild := rSetChildFields("loc", "elem", "idx"),

	eval := self >> self,
	#print := (self, i, si) >> Print(self.name, "
	isExpComposite := false
));

spiral.spl.RowVec.traversal := (self, sa, body) >> let(
			itr := var.fresh_t("itr", TInt),
			loopw(neq(struct_nth(sa, "index", itr), V(0)), chain(
				body
			))
	);

Class(RowVec2, RowVec, rec(
	init := self >> let(
		result := var.fresh_t("result", TInt),
		assign(result, V(0))
	),

	traversal := (self, sa, body) >> let(
			itr := var.fresh_t("itr", TInt),
			loopw(neq(struct_nth(sa, "index", itr), NULL(TInt)), chain(
				body
			))
	),
));

Class(DAGNode, TaggedNonTerminal,  rec(
    abbrevs :=  [ (nt, ylist, xlist) -> Checked(IsSPL(nt), [nt, When(IsList(ylist), ylist, [ylist]), When(IsList(xlist), xlist, [xlist])]) ],
    #transpose := self >> ObjId(self)(self.params[1].transpose(), self.params[2], self.params[3]).withTags(self.getTags())
));

Class(DAG, TCompose, rec(
    terminate := self >> Error("Not yet implemented."),

    from_rChildren := (self, rch) >> let(
        len := Length(rch),
        transposed := rch[len-1],
        tags := rch[len],
        t := ApplyFunc(ObjId(self), [rch{[1..len-2]}]),
        tt := When(transposed, t.transpose(), t),
        attrTakeA(tt.withTags(tags), self)
    ),

    rChildren := self >>
        Concatenation(self.params[1], [self.transposed, self.tags]),

    rSetChild := meth(self, n, newChild)
        local l;
        l := Length(self.params[1]);
        if n <= l then
            self.params[1][n] := newChild;
        elif n = l+1 then
            self.transposed := newChild;
        elif n = l+2 then
            self.tags := newChild;
        else Error("<n> must be in [1..", l+2, "]");
        fi;
        # self.canonizeParams(); ??
        self.dimensions := self.dims();
    end

));

INF := arg -> Cond(
   Length(arg)=0, Int(1000000),
   Length(arg)=1, var.fresh("i", TInt, toRange(Int(1000000))),
   Error("Usage: INF() | INF(<range>)")
);

#Class(INF, Value, rec(
#	__call__ := self >> SPL(WithBases(self, rec(
#		t := TInt,
#		v := V(1),
#		operations := ValueOps,
#	))),
#	isSPL := true,
#	dims := self >> self.v,
#	eval := self >> self,
#	print := (self, i, si) >> Print(self.name),
#));

IsMaskT := x ->  IsBound(x.isMaskT) and x.isMaskT;
IsAccumT := x ->  IsBound(x.isAccumT) and x.isAccumT;

Class(Mask, Diag, rec(
	__call__ := (self, func) >> SPL(WithBases(self, rec(
		element := func,
		var := func.var,
	))),
	isMaskT := true,
	print := (self, i, si) >> Print(self.name, "(", self.element, ")"),
	index := (self, i, j, dom) >> nth(self.var, add(mul(dom, i), j)),
	#find_index := (self, i, j, itr) >> decl([mitr], chain(assign(mitr, nth(m, @(2).val.var), loopw(logic_and()))))
));

Class(Accumulate, Diag, rec(
	__call__ := (self, func) >> SPL(WithBases(self, rec(
		element := func,
		var := func.var,
	))),
	isAccumT := true,
	print := (self, i, si) >> Print(self.name, "(", self.element, ")"),
));

Declare(SPLScope);

Class(SPLScope, Buf, rec(

	__call__ := (self, spl, scope) >> SPL(WithBases(self, rec(
		spl := Checked(IsSPL(spl), spl),
		scope := scope,
		_children := [spl],
		dimensions := spl.dimensions,
	))),
	dims := self >> Dimensions(self.spl),
	rChildren := self >> [self.spl, self.scope],
	rSetChild := rSetChildFields("spl", "scope"),
	#append := (self, L) >> Union(self.list, L),

	#freshScope := (self, inVar, outVar) >> SPL(WithBases(self, rec(
	#		input := inVar,
	#		output := outVar,
	#		_children := [],
	#	))),

	freshScope := (self) >> var.fresh_t("t", TPtr(TReal)),
));

Class(SpContainer, SumsBase, BaseContainer, rec(
	new := (self, spl) >> SPL(WithBases(self, rec(
		_children := [spl]))).setDims(),

	dims := self  >> self.child(1).dims()
));

Declare(TSparse_Matrix);

Class(PreDstruct, BaseMat, rec(
	__call__ := (self, count) >> SPL(WithBases(self, rec(
		count := count,
		dimensions := [count*count, count*count],
		_children := [count],
	))),
	isPermutaiton := False,
	isReal := True,
	rChildren := self >> [self.count],
	rSetChild := rSetChildFields("count"),
	doNotMarkBB := true,
	dims := self >> self.dimensions,
	print :=  (self, i, si) >> Print(self.name, "(", self.count,")"),
));
Class(PostDstruct, BaseMat, rec(
	__call__ := (self, count) >> SPL(WithBases(self, rec(
		count := count,
		dimensions := [count*count, count*count],
		_children := [count],
	))),
	isPermutaiton := False,
	isReal := True,
	rChildren := self >> [self.count],
	rSetChild := rSetChildFields("count"),
	doNotMarkBB := true,
	dims := self >> self.dimensions,
	print :=  (self, i, si) >> Print(self.name, "(", self.count,")"),
));

Class(HyperSprase, BaseOperation, rec(
	__call__ := arg >> let(
    self := arg[1],
    vector := arg[2],
	matrix := arg[3],
	WithBases(self, rec(vector := vector, matrix := matrix, _children := [vector, matrix], operations := PrintOps))),

   	isSparseT := true,
   	print := self >> Print(self.name, "(", self.vector, ", ", self.matrix, ")"),
   	dims := self >> self.vector.dims(),

));

Class(MatMul, BaseMat, rec(
	__call__ := (self, dimX, dimY) >> SPL(WithBases(self, rec(
		dimX := dimX,
		dimY := dimY,
		dimensions := [dimX*dimY, dimY*dimX],
	))),
	isPermutaiton := False,
	isReal := True,

	rChildren := self >> [self.dimX, self.dimY],
	rSetChild := rSetChildFields("dimX", "dimY"),

	dims := self >> self.dimensions,
	print :=  (self, i, si) >> Print(self.name, "(", self.dimX, ", ", self.dimY,")"),

));

Class(MakeDiag, Diag, rec(
	dims := self >> [self.element.domain(), self.element.domain()],
));

Class(Reduce, BaseMat, rec(
	__call__ := (self, traversal, size, t) >> SPL(WithBases(self, rec(
		traversal := Checked(IsString(traversal), traversal),
		size := size,
		t := Checked(IsType(t), t),
		dimensions := Cond(traversal = "Col", [1, size], [size, 1]),
	))),

	rChildren := self >> [self.traversal, self.t],
	rSetChild := rSetChildFields("traversal", "t"),

	dims := self >> self.dimensions,

	print :=  (self, i, si) >> Print(self.name, "(", self.traversal, ", ", self.size, ", ", self.t, ")"),

	isPermutaiton := False,
	isReal := True,

));

Class(SparseBlk, SumsBase, Mat, rec(
	new := (self, M) >> SPL(WithBases(self, rec(
			data_type := M.object,
            element := M.var,
            TType   := Cond( # NOTE: add checks to M
                            IsList(M.var),     UnifyTypes(List(Flat(M.var), InferType)),
                             IsValue(M.var),    M.t.t,
                           IsSymbolic(M.var), M.t.t),
                       ))).setDims(),
    area := self >> Length(Filtered(Flat(self.element), k -> k<>0)),
    new  := (self, M) >> SPL(WithBases(self, rec(element := M))).setDims(),
    dims := self >> Dimensions(self.element)
));

#Class(Trianglecount, TaggedNonTerminal, rec(
#	abbrevs := [(n) -> [n]],
#	isReal := self>>true,
#	dims := self >> [self.params[2]
#))

Declare(RulesSPLScope);
Declare(RulesMask);

Class(SpMV, TaggedNonTerminal, rec(
	abbrevs := [(n) -> [n]],
	isReal := self >> true,
	dims := self >> self.params[1],
	print := (self, i , si) >> Print(self.name, "(", self.params[1], ")")
));
Class(TVStack, TaggedNonTerminal, rec(
	abbrevs := [(n) -> [n]],
));
Class(TRowVec, TaggedNonTerminal, rec(
	abbrevs := [(n) -> [n]],
));

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

Class(Trace, TaggedNonTerminal, rec(
	abbrevs := [(n) -> [n]],
	isReal := self >> true,
	dims := self >> self.params[1].dims()[1] * self.params[1].dims()[2],
	print := (self, i , si) >> Print(self.name, "(", self.params[1], ")")
));

Class(MatrixMultiply, TaggedNonTerminal, rec(
	abbrevs := [(n,n, p) -> [n,n, p]],
	isReal := self >> true,
	dims := self >> self.params[1] * self.params[2],
	print := (self, i , si) >> Print(self.name, "(", self.params[1], ", ", self.params[2], Checked(self.params[3] <> "", Print(", ", self.params[3])), ")")
));

NewRulesFor(SpMV, rec(
	column_reduce := rec(
			info := "Column reduction based spmv",
			maxSize := false,
			applicable := (self, nt) >> true,
			apply := (nt, c, cnt) -> let(i := Ind(3),
					ts := TSparse(TArray(TInt, 3), TSemiring_Arithmetic(TInt)),
					fdataofs := FDataOfs(ts, 3, 0),
					md := MakeDiag(fdataofs),
					eye := I(i),
					sid := SUM(eye, md),
					matmul := MatMul(3,3),
					spmv := matmul * sid, spmv)
	)
));

#NewRulesFor(MatrixMultiply, rec(
#	 matmul := rec(
#		 info := "Matrix Multplication",
#		 maxSize := false,
#		 applicable := (self, nt) >> Cond(nt.params[3] = "", true, false),
#		 apply := (nt, c, cnt) -> let(i := Ind(3), j := Ind(3), scat := Scat(fTensor(fBase(i), fBase(j))),
#		 gath := Gath(fStack(fTensor(fBase(i), fId(3)), fTensor(fId(3), fBase(j)))), scp := SPLScope.freshScope(),
#		 kernel := scat * SPLScope(RowVec(FDataOfs(scp,6,V(3))), scp)  * gath, kernel2 := Rewrite(kernel, RulesSPLScope, opts), ISum(i, 3, ISum(j, 3, kernel2)))
#	 ),
#	 maskedmxm := rec(
#		 info := "Masked Matrix Multplication",
#		 maxSize := false,
#		 applicable := (self, nt) >> Cond(IsSPL(nt.params[3]), IsMaskT(nt.params[3]), false),
#		 apply := (nt, c, cnt) -> let(
#								i := Ind(3),
#								j := Ind(3),
#								scat := Scat(fTensor(fBase(i), fBase(j))),
#								gath := Gath(fStack(fTensor(fBase(i), fId(3)), fTensor(fId(3), fBase(j)))),
#								scp := SPLScope.freshScope(),
#								kernel := scat * SPLScope(RowVec(FDataOfs(scp,6,V(3))), scp) * gath,
#								mxm := ISum(i, 3, ISum(j, 3, kernel)),
#								maskedmxm := Mask(FDataOfs(M, 9, 0)) * mxm,
#								maskedmxm2 := Rewrite(maskedmxm, RulesMask, opts),
#								maskedmxm3 := Rewrite(maskedmxm2, RulesSPLScope, opts), maskedmxm3)
#	 )
#));

NewRulesFor(Trace, rec(
	csr_trace := rec(
		info := "Trace of CSR matrix",
		maxSize := false,
		applicable := (self, nt) >> true,#nt.params[1].dims()[1] =1,
		#apply := (nt, c, cnt) -> let(i := Ind(4), ISumAcc(Ind(), 4, ScatAcc(fId(1)) * Gath(fTensor(fBase(i), fBase(i))))),
		apply := (nt, c, cnt) -> let(i := Ind(5), ISumAcc(i, 5, Gath(fTensor(fBase(i), fBase(i))))),
	)
));

NewRulesFor(TSparse_Mat, rec(
	csr := rec(
		info := "Vertical stack of Row Vectors",
		maxSize := false,
		applicable := (self, nt) >> nt.params[1].dims()[1] = 1,
		#apply := (nt, c, cnt) -> VStack(nt.params[1]),
		apply := (nt, c, cnt) -> TSparse_Matrix(nt.params[1]),
	),
	csc := rec(
		info := "Horizontal stack of Column Vectors",
		maxSize := false,
		applicable := (self, nt) >> nt.params[1].dims()[2] = 1,
		apply := (nt, c, cnt) -> HStack(nt.params[1]),
	)
));

Class(DStructKernel, Grp);
Class(DStructKernelDot, Grp);
Class(GrpKernel, Grp);
Class(GrpKernel2, Grp);
Class(MaskKernelAxpy, Grp);
Class(MaskKernelDot, Grp);
Class(AccumKernel, Grp, rec(doNotMarkBB := true));
Class(TriCountKernel, Grp);
Class(PushBFSKernel, Grp);
Class(PullBFSKernel, Grp);
Class(BFSKernel, Grp);

TArrayBase.__call__ := (self, t, size) >>
        WithBases(self, rec(
        t    := Cond(Checked(IsType(t), IsSemiring(t)), t.t, t),
		ring := Cond(IsSemiring(t), t, "Error no ring"),
        size := Checked(IsPosInt0Sym(size), size),
        operations := TypOps));

#TArrayBase.rChildren := self >> Cond(IsBound(self.ring) and self.ring <> "Error no ring", [self.t, self.ring, self.size], [self.t, self.size]);
#TArrayBase.rSetChild := self >> Cond(IsBound(self.ring) and self.ring <> "Error no ring", rSetChildFields("t", "ring", "size"), rSetChildFields("t", "size"));
#TArrayBase.from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), rch);

Class(Ttrace, BaseOperation, rec(
	__call__ := (self, object) >> WithBases(self, rec(
        object := object,
        t := object.t,
        operations := TypOps)),

	print := self >> Print(self.name, "(", self.object, ")"),
	#size := self >> self.t.size,
	dims := self >> [1, 1],
	rChildren := self >> [self.object],
	rSetChild := rSetChildFields("object"),
));

IsSparseT := x -> IsType(x) and IsBound(x.isSparseT) and x.isSparseT;

Class(TSparse_Matrix, TArrayBase, rec(
   #__call__ := arg >> let(
	#    size := arg[2],
   #     element := TArray(arg[1], arg[2]),
   #     t := arg[1].t,
	#	prop := Cond(IsList(arg[3]), arg[3], [arg[3]]),
   #     operations := TypOps,WithBases(self, rec()),

	__call__ := arg >> let(
       self := arg[1],
       size := arg[2].size(),
	   element := TArray(arg[2], size),
       prop := Cond(IsList(arg[3]), arg[3], [arg[3]]),
	   t := element.t,
	   WithBases(self, rec(size := size, element := element, prop := prop, t := t, _children := [arg[2], prop], operations := PrintOps))),

   	isSparseT := true,
   	print := self >> Print(self.name, "(", self.element, ", ", self.prop, ")"),
   	dims := self >> [self.size, self.element.size],

	rChildren := self >> [self.element.t, self.size],
	rSetChild := rSetChildFields("t", "size"),

	range := self >> self.element.t.t,

	get_var := self >> let(var.fresh_t("spr_mat", TPtr(TInt))),
	traverse_outer := (self, i, j, var, body) >>  loopf(j, nth(var,i), nth(var, add(i, V(1))), chain(body)),
	index_row := (self, var, n, i) >> nth(var, add(n, add(V(1), i))),
	index_val := (self, var, n, j) >> nth(var, add(n, add(nth(var, n), add(j, V(1))))),
));

Declare(inref);

Class(TSparse, TArrayBase, rec(
	__call__ := (self, t, ring) >>
			WithBases(self, rec(
			t    := Checked(IsType(t), t),
			ring := ring,
			size2 := Checked(IsPosIntSym(t.size), t.size),
			operations := TypOps)),
    isSparseT := true,
    vtype := (self, v) >> TSparse(self.t.vtype(v), self.size2/v),
    toPtrType := self >> TPtr(self.t),
    doHashValues := true,
    #dims := self >> Cond(
    #    ObjId(self.t)=TSparse, [self.size] :: self.t.dims(),
    #    [self.size]),
	size := self >> self.size2,
	dims := self >> [1, self.size2],
	rChildren := self >> [self.t, self.ring],
	rSetChild := rSetChildFields("t", "ring"),
#	print := self >> Print(self.__name__, "(", self.t, ", ", self.size1, ", ", self.list1, ", ", self.size2, ", ", self.list2, ", ", self.sa, ")"),
	range := self >> self.t,
	print := self >> Print(self.__name__, "(", self.t, ", ", self.ring, ")"),

    get_var := self >> let(var.fresh_t("spr_arr", TPtr(TSparse(TArray(TInt,5), TSemiring_Arithmetic(TInt))))),

	length := (self, x) >> struct_nth(inref(x), "length", ""),

	get_elem_index := (self, x, idx) >> struct_nth(inref(x), "index", idx),

	get_elem_value := (self, x, idx) >> struct_nth(inref(x), "value", idx),

	traversal := (self, i, low, high, body) >> loopf(i, low, high, chain(body)),

	#num_nonzeros := (self, x) >> let(
	#	itr := var.fresh_t("itr", TInt),
	#	result := var.fresh_t("result", TInt),
	#	decl([itr, result], chain(
	#	assign(result, V(0)),
	#	loopw(neq(struct_nth(x,"value", itr), NULL(TInt)), chain(
	#	if1(neq(struct_nth(x,"value", itr), V(0)), chain(
	#		assign(result, add(result, V(1)))))))))),
));

Class(fItrStack, FuncClassOper, rec(
	__call__ := meth(arg)
	 local self, children, lkup, res, h;
        self := arg[1];
        children := Flat(Drop(arg, 1));
        if self.skipOneChild and Length(children)=1 then return children[1]; fi;

        h := self.hash;
        if h<>false then
            lkup := h.objLookup(self, children);
            if lkup[1]<>false then return lkup[1]; fi;
        fi;
        res := WithBases(self, rec(operations := RewritableObjectOps, _children := children));
        if h<>false then return h.objAdd(res, lkup[2]);
        else return res;
        fi;
    end,

	domain := self >> self._children[2].domain() * self._children[1].range,
    range := self >> self._children[2].range(),
    subdomainsDivisibleBy := (self, n) >> ForAll(self._children, x -> x.domain() mod n = 0),

    #lambda := meth(self)
    #    #local vals;
    #    #vals := [];
	#	for i in [1..self.domain()] do
	#	Add(vals, self.child(2).lambda().at(add(idiv(sub(i,V(1)), add(self.range(),V(1))), mul(self.range(), imod(sub(i,V(1)), add(self.range(), V(1)))))));
	#	od;
	#	return Lambda(self._children[1], vals);
    #end

	lambda := meth(self)
		local v, j;
		v := Ind(self.range());
		#return Lambda(v, self.child(2).lambda().at(add(idiv(v, idiv(self.range(), self.domain())), mul(self.range(), imod(v, idiv(self.range(), self.domain()))))));
		return let(l := Lambda(v, self.child(2).lambda().at(imod(v, self._children[2].domain()))),
					#c := Collect(l, @(1,var, e -> e.range = self._children[2].domain())),
					#c := Collect(l, @(1, add, e -> Error())),
					s := SubstVars(Copy(l), rec((self.child(1).id) := idiv(v, self.child(2).domain()))),
					s);
	end
));

Class(RulesSPLScope, RuleSet);
Class(RulesINF, RuleSet);
Class(RulesSymbol, RuleSet);
Class(RulesTrace, RuleSet);
Class(RulesMG, RuleSet);
Class(RulesMR, RuleSet);
Class(RulesMR2, RuleSet);
Class(RulesMask, RuleSet);
Class(RulesSparseMaskAxpy, RuleSet);
Class(RulesSparseMaskDot, RuleSet);
Class(RulesIDiag, RuleSet);
Class(RulesMatMul, RuleSet);
Class(RulesSigSPMV1, RuleSet);
Class(RulesSigSPMV2, RuleSet);
Class(RulesSigSPMV3, RuleSet);
Class(RulesSigSPGEMM1, RuleSet);
Class(RulesSigSPGEMM2, RuleSet);
Class(RulesAccumulate, RuleSet);
Class(RulesTriCount, RuleSet);
Class(RulesScatRow, RuleSet);

SPMVStrategy := [RulesIDiag, RulesMatMul, RulesSigSPMV1, RulesSigSPMV2];
SPGEMMStrategy := [RulesSigSPGEMM1];

Class(Tcsr, AtomicTyp);

#														fdata := @(6).val.element.params[3],  spm := fdata.get_var(), inside := @(3).val._children[1], v1 := @(1).val.var,
#														v2 := @(2).val.var, v3 := @(3).val.var, n := @(1).val.domain,
#result := ISum(v1, n, ISum(v2, fdata.index_row(spm, n, v1), ISum(v3, nth(X, add(nth(spm, add(n, add(V(1), v2))), V(1))), inside))),
#														result

RewriteRules(RulesTriCount, rec(
	add_loop_bounds := Rule([@@(1,ISum),[@(2, ISum, e-> IsValue(e.domain) = true or IsInt(e.domain)), @(3,RowVec)]], (e,cx)->let(fdata := Collect(e, FDataSparseOfs)[1], v1 := e.var,v2 := @(2).val.var,
																			TriCountKernel(ISum(v1, cx.opts.symbol[1], ISum(v2, sub(opts.sparse_mat.index_row(X, cx.opts.symbol[1], nth(X,add(v1, V(1)))), opts.sparse_mat.index_row(X, cx.opts.symbol[1], nth(X,v1))), @(3).val))))),
));

RewriteRules(RulesSigSPGEMM1, rec(
	add_icode_snippets := Rule([@@(1,ISum),[@(2,ISum), [@(3, ISum), [@(4, Compose), @(5,Scat), @(6, Diag, e->Length(Collect(e, FDataOfs)) > 0), @(7, Gath)]]]], (e,cx)->let(fdata := @(6).val.element.params[3],  spm := fdata.var.get_var(), Append(cx.opts.symbol, [spm]),
																															v1 := e.var,v2 := @(2).val.var, v3 := @(3).val.var, n := cx.opts.symbol[1],
																															newfdata := FDataSparseOfs(fdata.var, fdata.len,fdata.var.index_val(spm, cx.opts.symbol[1], v2)),
																															#d := Diag(fConst(v2.t, 1, fdata.var.index_val(spm, n, v2))),
																															d := Diag(fConst(v2.t, 1, newfdata)),
																															func := fTensor(fBase(v3), fBase(v2)),
																															sa := ScatAcc(func), g := Gath(func),
																															res := Cond(cx.opts.YType.t = Tcsr, GrpKernel2(ISum(v1, n, ISum(v2, sub(nth(spm, add(v1, V(1))), nth(spm, v1)), ISum(v3, sub(nth(X, add(nth(spm, add(n, add(V(1), v2))), V(1))), nth(X, nth(spm, add(n, add(V(1), v2))))), sa * d * g)))),
																															GrpKernel(ISum(v1, n, ISum(v2, sub(nth(spm, add(v1, V(1))), nth(spm, v1)), ISum(v3, sub(nth(X, add(nth(spm, add(n, add(V(1), v2))), V(1))), nth(X, nth(spm, add(n, add(V(1), v2))))), sa * d * g))))),
																															res)),
));

RewriteRules(RulesSigSPGEMM2, rec(
	add_grpkernel := Rule([@(1,ISum),[@(2,ISum), [@(3, ISum), [@(4, Compose), @(5,ScatAcc), @(6, Diag), @(7, Gath)]]]], e -> ISum(@(1).val.var, @(1).val.domain, GrpKernel(@(2).val))),
));

RewriteRules(RulesIDiag, rec(
	convert_I := ARule(SUM, [@(1,I), @(2, MakeDiag)], e-> [let(j := Ind(@(1).val.obj.size), ivs := IterVStack(j, @(2).val.element.len, Gath(fTensor(fId(@(2).val.element.len), fBase(j)))), SUM(ivs, @(2).val))]),
));

RewriteRules(RulesMatMul, rec(
	consume_matmul := ARule(Compose, [@(1,MatMul), [@(2,SUM), @(3,IterVStack), @(4,MakeDiag)]], e->[let(f := @(4).val.element, i := @(3).val.var, g := @(3).val._children[1],
																					IterVStack(i, f.len, Diag(fConst(i.t, f.len, f)) * g.toloop(1)))]),
));

RewriteRules(RulesSigSPMV1, rec(
	remove_extra_isum := Rule([@(1,Compose), @(2, ISum), @(3,ISum)], e->@(3).val),
));

RewriteRules(RulesSigSPMV2, rec(
		rebuild_with_var := Rule([@@(1,ISum), [@(2,ISum), [@(3,Compose), @(4,Scat), @(5,Diag, e-> Length(Collect(e, FDataOfs)) > 0 and Length(Collect(e, TSparse)) > 0),@(6,Gath)]]], (e,cx)-> let(f := @(5).val.element.params[3], v1 := e.var, v2 := @(2).val.var,
																												spa := f.var.get_var(),  Append(cx.opts.symbol, [spa]),
																												newfdata := FDataSparseOfs(f.var, f.len, f.var.get_elem_value(spa, v1)),
																												d := Diag(fConst(v2.t, 1, newfdata)),
																												func := fTensor(fBase(v2), fBase(v1)),
																												sa := ScatAcc(func), g := Gath(func),
																												result := GrpKernel(ISum(v1, f.var.length(spa), ISum(v2, sub(nth(X, add(f.var.get_elem_index(spa, v1), V(1))), nth(X, f.var.get_elem_index(spa, v1))), sa * d * g))),
																												result)),
		rebuild_dense := Rule([@@(1,ISum), [@(2,ISum), [@(3,Compose), @(4,Scat), @(5,Diag, e-> Length(Collect(e, FDataOfs)) > 0 and Length(Collect(e, TSparse)) = 0),@(6,Gath)]]], (e,cx)-> let(f := @(5).val.element.params[3], v1 := e.var, v2 := @(2).val.var,
																												spa := f.var,  Append(cx.opts.symbol, [spa]),
																												d := Diag(fConst(v2.t, 1, nth(spa, v1))),
																												func := fTensor(fBase(v2), fBase(v1)),
																												sa := ScatAcc(func), g := Gath(func),
																												result := GrpKernel(ISum(v1, f.len, ISum(v2, sub(nth(X, add(v1, V(1))), nth(X, v1)), sa * d * g))),
																												result)),
));

RewriteRules(RulesSigSPMV3, rec(
		add_grpkernel := Rule([@(1,ISum), [@(2,ISum), [@(3,Compose), @(4,ScatAcc), @(5,Diag),@(6,Gath)]]], e-> GrpKernel(@(1).val)),
));

RewriteRules(RulesScatRow, rec(
		remove_scat := Rule([@(1,Compose), @(2,Scat), @(3, RowVec, e-> IsBound(e.element.var.t.ring) and e.element.var.t.ring <> "Error no ring")], e-> @(3).val),
));

#																								func := fTensor(fId(1), fBase(@(1).val.var)),
#																								sa := ScatAcc(func), g := Gath(func),
#comments for rulesmxmtrace
#cis := Rule([ISum, [ISum, @(1)]], e -> [let(rec(result := ISumAcc(e.var, e.domain, @(1).val)), result)]),
#remove_trace := ARule(Compose, [@(1, RowVec), @(2, Gath), @(3,ISum)], e->[@(3).val]),
#change_isum_isumacc := ARule(ISum, [@(1, SPLScope)], e -> [let(updom := e.var, inputs := e._children[1], isa := ISumAcc(updom, updom.range, inputs), isa)]),
#change_move_scat := ARule(ISum, [@(1, SPLScope)], e -> [Error()]),
#remove_trace := ARule(Compose, [@(1, RowVec), @(2, Gath), @(3,ISum)], e->[@(3).val]),
#remove_trace := ARule(Compose, [@(1, Gath), [@(2, ISum), @(3, ISum, e -> let(s := Collect(Copy(e._children[1]._children), Scat), s[1].func.__name__ = "fTensor"
	#																			and Length(s[1].func._children) = 2 and s[1].func._children[1].__name__ = "fBase"
	#																			and s[1].func._children[2].__name__ = "fBase"))]],
	#																			e->[let(v1 := @(2).val.var, v2 := @(3).val.var, f := fBase(@(2).val.var), c := @(3).val._children[1], two := @(3).val._children[1]._children[2] * @(3).val._children[1]._children[3],
	#																			ISum(v1, v1.range, ISum(v2, v2.range, COND(fBase(v1), Scat(fBase(v1)) * two, O(c.dims()[1],1)))))]),
RewriteRules(RulesMG, rec(
	move_gath := ARule(Compose, [@(1, Gath), [@(2, ISum), @(3, ISum, e -> let(s := Collect(Copy(e._children[1]._children), Scat), Length(s) = 1))]],
																							e -> [let(exp := @(1).val * @(3).val._children[1],
																									ISum(@(2).val.var, @(2).val.var.range,
																									ISum(@(3).val.var, @(3).val.var.range,
																								  	exp)))]),

	gath_scat_cond := ARule(Compose, [@(1, Gath), @(2, Scat)], e->[let(v1 := @(2).val.func._children[1].params[2], v2 := @(2).val.func._children[2].params[2],
																		COND(eq(v1, v2),
																			Scat(@(2).val.func._children[1]),
																			O(@(2).val.func._children[1].params[1], 1)))]),
));

RewriteRules(RulesMR, rec(
	collapse_loop_cond := Rule([@(1,ISum), [@(2,ISum), [@(3, Compose), @(4,COND),...]]], e->let(v1 := @(1).val.var, v2 := @(2).val.var, com := @(3).val,
																								com2 := SubstTopDown(Copy(com), @(5, COND), g -> @(5).val._children[1]),
																								com3 :=  SubstVars(Copy(com2), rec((v2.id) := v1)),
																								ISum(v1, v1.range, com3))),
));

RewriteRules(RulesMR2, rec(
	move_rowvec := ARule(Compose, [@(1, RowVec), [@(2, ISum), [@(3, Compose), @(4,Scat), ...]]], e -> [let(v1 := @(2).val.var, s := @(3).val._children[2], f := @(3).val._children[3],
																											rv := RowVec(fCompose(@(1).val.element, @(4).val.func)),
																											ISumAcc(v1, v1.range, rv * s * f))]),
));

#e->[let(v1 := @(2).val.var, v2 := @(3).val.var,
#inside := @(3).val._children[1], dom := @(2).val.domain,
#con := COND(eq(@(1).val.index(v1,v2,dom), 1)))]),

RewriteRules(RulesAccumulate, rec(
	add_accumkernel := ARule(Compose, [@(1, Accumulate), @(2, [GrpKernel, ISum])], e->[AccumKernel(@(2).val)]),
));

RewriteRules(RulesSparseMaskAxpy, rec(
	move_mask := ARule(Compose, [@(1,Mask, e->IsSparseT(e.element.var.t)), [@(2, GrpKernel), [@(3,ISum), [@(4,ISum), @(5,ISum)]]]], (e,cx)-> [let(exp := @(1).val * @(5).val._children[1],
																																		ISum(@(3).val.var, @(3).val.domain,
																																		ISum(@(4).val.var, @(4).val.domain,
																																		exp)))]),

	move_mask2 := ARule(Compose, [@(1,Mask, e->IsSparseT(e.element.var.t)), [@(2, GrpKernel2), [@(3,ISum), [@(4,ISum), @(5,ISum)]]]], (e,cx)-> [let(exp := @(1).val * @(5).val._children[1],
																																		GrpKernel2(ISum(@(3).val.var, @(3).val.domain,
																																		ISum(@(4).val.var, @(4).val.domain,
																																		exp))))]),
	#self.symbol := Concat(self.symbol, [self.sparse_mat.get_var()]), cx.opts.symbol := Concat(cx.opts.symbol,[spm2]),
	add_mask_kernel := ARule(Compose, [@(1,Mask, e->IsSparseT(e.element.var.t)), @(2, ScatAcc), @(3,Diag), @(4, Gath)], e->[let(mask2 := Mask(FDataOfs(@(1).val.element.var.get_var(), @(1).val.element.len, @(1).val.element.ofs)), MaskKernelAxpy(mask2 * @(2).val * @(3).val * @(4).val))]),
));

RewriteRules(RulesSparseMaskDot, rec(
	convert_mask := ARule(Compose, [@(1,Mask, e->IsSparseT(e.element.var.t)), [@(2,ISum), @(3,ISum)]], e->[let(mask2 := Mask(FDataOfs(@(1).val.element.var.get_var(), @(1).val.element.len, add(@(2).val.var, @(3).val.var))),
																											exp := mask2 * @(3).val._children[1],
																											ISum(@(2).val.var, @(2).val.domain,
																											ISum(@(3).val.var, @(3).val.domain,
																											MaskKernelDot(exp))))]),
));

RewriteRules(RulesMask, rec(
	move_mask := ARule(Compose, [@(1,Mask), [@(2, ISum), @(3, ISum)]], e -> [let(exp := @(1).val * @(3).val._children[1],
																			ISum(@(2).val.var, @(2).val.var.range,
																			ISum(@(3).val.var, @(3).val.var.range,
																			exp)))]),
	replace_mask := Rule([@(1,ISum), [@(2,ISum), [@(3, Compose), @(4,Mask),...]]], e -> let(v1 := @(1).val.var, v2 := @(2).val.var, dom := @(2).val.domain,
																	inside := @(3).val, inside2 := SubstTopDown(Copy(inside), @(5, Mask), g -> I(@(5).val.element.len)),
																	con := COND(eq(@(4).val.index(v1,v2,dom), 1), inside2, Scat(fTensor(fBase(v1), fBase(v2))) * O(1,1)), ISum(v1, v1.range, ISum(v2, v2.range, con)))),
));

#comments for rulestrace
#change_fbase := ARule(fCompose, [@(1,itrfStack, e -> let(one := e._children[2]._children[1].domain(), two := e._children[2]._children[2].domain(), one = two)), @(2, fBase)], e -> [let(c := Collect(e, var), v2 := c[1], SubstTopDown(e, @(1,var, e-> var.range <> v2.range), e->v2))]),#[fCompose(@(1).val, @(1).val._children[2]._children[1])]),
RewriteRules(RulesTrace, rec(
	change_fbase := ARule(Gath, [@(1, fCompose, e -> let(length := Length(e._children), length = 2))],
									 e -> [let(c := Collect(e, var),
												v2 := c[Length(c)],
												result := SubstTopDown(Copy(e), @(1, var, e-> var.range <> v2), e -> v2),
												Error())]),
));

#comments for RulesScope
#newspl2 := SubstTopDown(Copy(newspl), @(1, var, p -> p = s.scope), e -> X),
#rowvec_gath := ARule(Compose, [@(1, RowVec), @(2, Gath)],e -> [let(i := Ind(@(1).val.element._children[1].len), ISum(i, e._children[2].element._children[1].len, diagMul(e._children[2].element,e._children[3].func)))]),
#[let(scat := e._children[1], s := e._children[2], g := e._children[3].func, newscope := SPLScope(RowVec(diagMul(fCompose(s.spl.element._children[1], g._children[2]))), s.scope), newscope * Gath(g._children[1]))]),
RewriteRules(RulesSPLScope, rec(
	move_inside_gath := ARule(Compose, [@(1, SPLScope, e->
		let(check1 := Collect(e.spl, FDataOfs), check2 := Collect(e.spl, @@(1, var, (s, cx) -> s = e.scope and cx.FDataOfs[1].var <> s)), Length(check1) > 0 and Length(check2) = 0)),
		[@(2,Gath, e-> Length(e.func._children) = 2 and e.func._children[1].domain() = e.func._children[2].domain()), @(3, fStack)]],
			e -> [let(scat := e._children[1],
					s := e._children[2],
					g := e._children[3].func,
					fdataofs := FDataOfs(s.spl.element.var, s.spl.element.len/2, V(0)),
					newspl := RowVec(fCompose(fdataofs, g._children[2])),
					newscope := SPLScope(newspl * Gath(g._children[1]), s.scope),
					newscope)]),

	move_inside_scat := ARule(Compose, [@(1, Scat), @(2, SPLScope)], e-> [SPLScope(@(1).val * @(2).val.spl, @(2).val.scope)]),
));

RewriteRules(RulesStrengthReduce, rec(
	add := ARule(add, [@1, INF()], e -> INF),
	fTen := ARule(fTensor, [@(1, Ind(INF()))], e -> V(1)),
	fBas := ARule(fBase, [@(1, Ind(INF()))], e -> V(1)),
	fI := ARule(fId, [@(1, Ind(INF()))], e -> V(1)),
	fSta := ARule(fStack, [@(1, Ind(INF()))], e -> V(1)),
	#Ind := ARule(Ind, [@(1, INF)], e -> V(1)),
	#fdo := ARule(FDataofs, [@(1, INF)], e -> V(1))
));

#RewriteRules(RulesSymbol, rec(
#	addsymbol := Cond(IsBound(opts.symbol), SubstTopDown(Copy(cs), [@(1,nth), ..., @(2,var, e-> e=X), ...],
#    e -> nth(opts.symbol[1], e.idx)), cs);
#));

Class(T_Sparse, T_Type, rec(
       t := TInt,
));

Class(NULL, Command, rec(
__call__ := (self, t) >> WithBases(self, rec(
	t    := Checked(IsType(t), t),
    operations := CmdOps
  )),
 	rChildren := self >> [],
	print := (self, i, si) >> Print(self.__name__)
));

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

Class(inref, nth, rec(
	__call__ := (self, loc) >> Inherited(loc, TInt.value(0)),
    rChildren := self >> [self.loc],
    rSetChild := rSetChildFields("loc"),
));

Class(if4, Command, rec(
  __call__ := (self, if_cond, if_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(add_to_sparse, Command, rec(
	__call__ := (self, yy, val, middle, idx, row) >> WithBases(self, rec(
		yy := yy,
		val := val,
		middle := middle,
		idx := idx,
		row := row,
		operations := CmdOps
	)),
	rChildren := self >> [self.yy, self.val, self.middle, self.idx, self.row],
  	rSetChild := rSetChildFields("yy", "val", "middle", "idx", "row"),
	print := (self, i, si) >> Print(self.__name__, "(", self.yy, ", ", self.val, ", ", self.middle, ", ", self.idx, ", ", self.row, ")")
));

#For the data structure
Class(Tdstruct, AtomicTyp);
Class(Ttree, AtomicTyp);
Class(next_ptr, Command, rec(
	__call__ := (self, var, value) >> WithBases(self, rec(
		var := var,
		value := value,
		operations := CmdOps
	)),
	rChildren := self >> [self.var, self.value],
  	rSetChild := rSetChildFields("var", "value"),
	print := (self, i, si) >> Print(self.__name__, "(", self.var, ", ", self.value, ")")
));
Class(allocate_ds, Command, rec(
	__call__ := (self, count) >> WithBases(self, rec(
		count := count,
		operations := CmdOps
	)),
	rChildren := self >> [self.count],
  	rSetChild := rSetChildFields("count"),
	print := (self, i, si) >> Print(self.__name__, "(", self.count,")")
));
Class(insert_ds, Command, rec(
	__call__ := (self, dstruct, loc, val, root) >> WithBases(self, rec(
		dstruct := dstruct,
		loc := Cond(IsValue(loc), loc, V(loc)),
		val := Cond(IsValue(val), val, V(val)),
		root := root,
		operations := CmdOps
	)),
	rChildren := self >> [self.dstruct, self.loc, self.val, self.root],
  	rSetChild := rSetChildFields("dstruct", "loc", "val", "root"),
	print := (self, i, si) >> Print(self.__name__, "(", self.dstruct, ", ",
        self.loc, ", ", self.val, ", ", self.root, Blanks(i), ")")
));
Class(ds_to_csr, Command, rec(
	__call__ := (self, dstruct) >> WithBases(self, rec(
		dstruct := dstruct,
		operations := CmdOps
	)),
	rChildren := self >> [self.dstruct],
  	rSetChild := rSetChildFields("dstruct"),
	print := (self, i, si) >> Print(self.__name__, "(", self.dstruct, Blanks(i), ")")
));
Class(dstructToTree, Command, rec(
	__call__ := (self, dstruct, count) >> WithBases(self, rec(
		dstruct := dstruct,
		count := count,
		operations := CmdOps
	)),
	rChildren := self >> [self.dstruct, self.count],
  	rSetChild := rSetChildFields("dstruct","count"),
	print := (self, i, si) >> Print(self.__name__, "(", self.count, ", ",
        self.dstruct,")")
));
Class(destory_ds, Command, rec(
	__call__ := (self, dstruct, root) >> WithBases(self, rec(
		dstruct := dstruct,
		root := root,
		operations := CmdOps
	)),
	rChildren := self >> [self.dstruct, self.root],
  	rSetChild := rSetChildFields("dstruct", "root"),
	print := (self, i, si) >> Print(self.__name__, "(", self.dstruct, ", ", self.root, Blanks(i), ")")
));

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

#only supports types with .ctype TInt TReal
CUnparser.malloc := (self, o, i, is) >> Print("(", o.t.ctype, "*)malloc(", o.size, "*", self(sizeof(o.t), i, is), ")");

#Unparser for data structure
CUnparser.allocate_ds := (self, o, i, is) >> Print("allocate_ds(", self(o.count, i, is), ")");
CUnparser.insert_ds := (self, o, i, is) >> Print(Blanks(i), "insert_ds(", o.dstruct, ", ", self(o.loc, i, is), ", ", self(o.val, i, is), ", ", o.root, ")", ";\n");
CUnparser.ds_to_csr := (self, o, i, is) >> Print("ds_to_csr(", o.dstruct, ")");
CUnparser.destory_ds := (self, o, i, is) >> Print(Blanks(i), "destory_ds(", o.dstruct, ", ", o.root, ")");
CUnparser.dstructToTree := (self, o, i, is) >> Print("dstructToTree(",self(o.dstruct, i, is), ", ", self(o.count, i, is), ")");
CUnparser.Tdstruct := (self,t, vars,i,is) >> Print("dstruct ", self.infix(vars, ", ",i+is));
CUnparser.Ttree := (self,t, vars,i,is) >> Print("tree ", self.infix(vars, ", ",i+is));
CUnparser.Tcsr := (self, t, vars, i, is) >> Print("struct csr",self.infix(vars, ", ", i + is));

CUnparser.inref := (self, o, i, is) >> Print("(*", self(o.loc, i, is), ")");

CUnparser.TSparse := (self,t, vars,i,is) >> Print(Blanks(i),
	"struct sparse_arr", self.infix(vars, ", ", i + is));

#CUnparser.sparse_nth3 := (self, o, i, is) >> o.list.get_elem_value(o.list.get_var(), o.idx);

CUnparser.next_ptr := (self, o, i, is) >> Print(o.var, "->", o.value);

CUnparser.add_to_sparse := (self, o, i, is) >> Print(Blanks(i), "add_to_sparse(", self(o.yy, i, is), ", ", o.val, ", ", o.middle, ", ", o.idx, ", ", o.row, ");\n");

CUnparser.struct_nth := (self,o,i,is) >> Cond(o.idx <> "", Print(
	self(o.loc,i,is), ".", o.elem, "[", o.idx, "]"),
	Print(self(o.loc,i,is), ".", o.elem));

CUnparser.if4 := (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.NULL := (self, o, i, is) >> Print(Blanks(i), "NULL");

CUnparser.loopn := (self, o, i, is) >> Cond(IsBound(o.range.args) and Length(o.range.args) > 1 and Length(Collect(o.range, sub)) > 0, let(v := o.var.id, lo := o.range.args[2], hi := o.range.args[1],
        Print(Blanks(i), "for(",self.declare(o.var.t, o.var, i, is), " = ", self(lo,i,is), "; ", v, " < ", self(hi,i,is), "; ", v, "++) {\n",
          self(o.cmd,i+is,is),
          Blanks(i), "}\n")), let(v := o.var.id, lo := 0, hi := o.range,
        Print(Blanks(i), "for(int ", v, " = ", lo, "; ", v, " < ", self(hi,i,is), "; ", v, "++) {\n",
          self(o.cmd,i+is,is),
          Blanks(i), "}\n")));

#Class(sparse_nth, Command, rec(
#  __call__ := (self, loc, idx) >> WithBases(self, rec(
#	loc := loc,
#    idx := idx,
#	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), ")"
#  )
#));
#
#CUnparser.sparse_nth := (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");

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

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

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

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

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

DefaultSumsGen.sparse_nth3 := (self, o, opts) >> o.element.get_elem_value(o.element.get_var(), o.idx);

#DefaultSumsGen.RowVec := (self, o, opts) >> let(n := Ind(o.element._children[1].len), Error(), ISumAcc(n, o.element._children[1].len, ScatAcc(o.element)));
#let(c := Collect(self, @(1,var, e -> IsSparseT(e))), Length(c) > 0)
DefaultSumsGen.RowVec := (self, o, opts) >> Cond(opts.name = "SparseOpts" or Length(Collect(o, TSparse)) > 0 or IsBound(o.element.var) and IsBound(o.element.var.t.t) and IsSparseT(o.element.var.t.t), o, let(i := Ind(o.element.domain()), ISumAcc(i, o.element.domain(), Scat(fId(1)) * Blk1(o.element.at(i))) * Gath(fBase(i))));

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

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

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

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

DefaultSumsGen.Reduce := (self, o ,opts) >> RowVec(fConst(o.t, o.size, 1));

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

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

Class(DataStructureCodegenMixin, DefaultCodegen, rec(
	Scat := meth(self, o, y, x, opts)
		local index_i, index_j, index_k, loc;
		index_i := o.func._children[1].params[1];
		index_j := o.func._children[2].params[2];
		index_k := o.func._children[1].params[2];
		#first := add(index_j, nth(opts.symbol[4], index_i));
		#second := add(index_k, nth(X, opts.sparse_mat.index_row(opts.symbol[4], opts.symbol[1], first)));
		loc := opts.sparse_mat.index_row(X, opts.symbol[1], index_k);
		return insert_ds(y, loc, nth(x,0), opts.symbol[3]);
	end
));

#DefaultSumsGen.IterVStack := (self, o, opts) >> Cond(let(c := Collect(o, TSparse), Length(c) > 0), o, let(
#	bkcols := Cols(o.child(1)),
#	bkrows := Rows(o.child(1)),
#	nblocks := o.domain,
#	cols := Cols(o), rows := Rows(o),
#	ISum(o.var, o.domain,
#	    Scat(fTensor(fBase(nblocks, o.var), fId(bkrows))) *
#	    self(o.child(1), opts) *
#	    Gath(fId(bkcols)no))
#    ));

DefaultSumsGen.SPLScope := (self, o, opts) >> SPLScope(self(o.child(1), opts), o.scope);

DefaultCodegen.SPLScope := (self, o, y, x, opts) >>
	decl([o.scope], chain(assign(o.scope, x), self(o.child(1), y, x, opts)));

	#decl([o.scope], chain(assign(o.scope, x), self(o.child(1), o.scope, x, opts))); #return decl([o.scope], self(o.child(1)._children[1], o.scope.output, x, opts))); #(self(o.child(1), y, x, opts));
	#return decl([o.scope], chain(self(o.child(1)._children[2], o.scope, x, opts), self(o.child(1)._children[1], y, o.scope, opts))); #(self(o.child(1), y, x, opts));

#Overridden Codegen from SpiralDefaults
DefaultCodegen.ISumAcc := (self, o, y, x, opts) >> let(ii := Ind(), chain(loop(ii, Rows(o), assign(nth(y, ii), V(0))), loopf(o.var, V(0), o.domain, self._acc(self(o.child(1), y, x, opts), y))));

DefaultCodegen.Diag := (self, o, y, x, opts) >> Cond(Length(Collect(o, TSparse)) > 0,
   #let(elem := fConst(o.element.params[3].ofs.loc.t, 1, o.element.params[3].var.index_val(o.element.params[3].ofs.loc, opts.symbol[1], o.element.params[3].ofs.idx)),
   let(elem := fConst(o.element.params[1], o.element.params[2], o.element.params[3].ofs),
   ring := Collect(o.element, TSparse)[1].ring, i := Ind(),
   elt := elem.lambda(),
   loop(i, elt.domain(), assign(nth(y, i), ring.product(elt.at(i),nth(x, i))))),
   Cond(Length(Collect(o, TSparse)) = 0 and IsBound(o.element.params[3].loc.t.ring),let(
   ring := o.element.params[3].loc.t.ring,
   elt := o.element.lambda(),
   loop(i, elt.domain(), assign(nth(y, i), ring.product(elt.at(i),nth(x, i))))),
   let(i := Ind(),
   elt := o.element.lambda(),
   loop(i, elt.domain(), assign(nth(y, i), elt.at(i) * nth(x, i))))));

DefaultCodegen.IterVStack := meth(self, o, y, x, opts)
	local v, i, j, itr;
	v := o._children[1].element.var.get_var();
	itr := var.fresh_t("itr", TInt);
	i := Ind();
	j := Ind();
	return decl([v,i,j,itr], chain(assign(itr, V(0)), o._children[1].element.var.traversal(i, V(0), opts.symbol[1],
				if1(eq(o._children[1].element.var.get_elem_index(v, itr), i), chain(
					opts.sparse_mat.traverse_outer(i,j,x,assign(nth(y, nth(x, add(opts.symbol[1], add(V(1), j)))), add(nth(y, nth(x, add(opts.symbol[1], add(V(1), j)))),
					mul(opts.sparse_mat.index_val(x,opts.symbol[1], j), o._children[1].element.var.get_elem_value(v, itr))))),
					assign(itr, add(itr, V(1))))))));
end;

DefaultCodegen.RowVec := (self, o, y, x, opts) >> Cond(let(c := Collect(o, TSparse), Length(c) > 0) and IsBound(opts.vec) and opts.vec = true and IsBound(opts.mat) and opts.mat = false,
	let(i := Ind(), j := Ind(),v := o.element.var.get_var(), Append(opts.symbol, [v]), t := var.fresh_t("t", TInt),
	decl([i, j, t], chain(assign(t, V(0)), assign(i, V(0)), assign(j, V(0)), loopw(logic_and(lt(i, o.element.var.length(v)), lt(j, o.element.var.length(x))),
	if3(lt(o.element.var.get_elem_index(v, i), o.element.var.get_elem_index(x, j)),
			assign(i, add(i, V(1))),
			lt(o.element.var.get_elem_index(x, j), o.element.var.get_elem_index(v, i)),
			assign(j, add(j, V(1))),
			chain(assign(t, mul(o.element.var.get_elem_value(v, i), o.element.var.get_elem_value(x, j))),
			assign(i, add(i, V(1))), assign(j, add(j, V(1)))))), assign(nth(y, 0), t)))),
let(c := Collect(o, TSparse), Length(c) > 0) and IsBound(opts.vec) and opts.vec = false and IsBound(opts.mat) and opts.mat = false,
	let(i := Ind(), v := o.element.var.get_var(), t := var.fresh_t("t", TInt),
		decl([v, i, t], chain(assign(t, V(0)), assign(i, V(0)), loopw(lt(i, o.element.var.length(v)),
		chain(assign(t, mul(o.element.var.get_elem_value(v, i), nth(x,o.element.var.get_elem_index(v, i)))),
		assign(i, add(i, V(1))))), assign(nth(y, 0), t)))),
let(c := Collect(o, TSparse), Length(c) = 0) and IsBound(opts.vec) and opts.vec = true and IsBound(opts.mat) and opts.mat = false,
	let(i := Ind(), Cond((o.element.var in opts.symbol) = false, Append(opts.symbol, [o.element.var]), skip()), ring := Cond(IsBound(o.element.var.t.ring) and o.element.var.t.ring <> "Error no ring", o.element.var.t.ring, TSemiring_Arithmetic(o.element.var.t)),
		decl([], chain(loopn(i, sub(nth(x, add(o.element.ofs,V(1))), nth(x,o.element.ofs)), chain(assign(deref(y),
		ring.sum(deref(y), ring.product(opts.sparse_mat.index_val(x, opts.symbol[1], i),
		nth(o.element.var,opts.sparse_mat.index_row(x,opts.symbol[1],i)))))))))),
let(c := Collect(o, add), Length(c) > 0) or let(c2 := Collect(o, TSparse), Length(c2) > 0)  and IsBound(opts.mat) and opts.mat = true and IsBound(opts.vec) and opts.vec = false,
	let(itr1 := var.fresh_t("itr", TPtr(TInt)), itr2 := var.fresh_t("itr", TPtr(TInt)), col1 := var.fresh_t("xcol", TInt),
		col2 := var.fresh_t("spmcol", TInt), m2 := o.element.var.get_var(), When(y.id <> "tempalg1", Append(opts.symbol, [m2])), i := o.element.ofs.args[1],
		j := o.element.ofs.args[2], decl([itr1, itr2], chain(assign(itr1, add(x,i)), assign(itr2, add(m2, j)),
		loopw(logic_and(neq(itr1, add(x, add(i,V(1)))), neq(itr2,add(m2, add(j,V(1))))), decl([col1, col2], chain(assign(col1, opts.sparse_mat.index_row(x,opts.symbol[1], deref(itr1))),
		assign(col2,opts.sparse_mat.index_row(m2,opts.symbol[1], deref(itr2))), if3(gt(col1, col2), assign(itr2, add(itr2, V(1))), lt(col1, col2), assign(itr1, add(itr1, V(1))),
		chain(assign(nth(y, add(mul(i, opts.symbol[1]), j)), o.element.var.t.ring.sum(nth(y, add(mul(i, opts.symbol[1]), j)), o.element.var.t.ring.product(opts.sparse_mat.index_val(x, opts.symbol[1], deref(itr1)), opts.sparse_mat.index_val(m2, opts.symbol[1], deref(itr2))))),
		assign(itr1, add(itr1, V(1))), assign(itr2, add(itr2, V(1))))))))))),
let(i := Ind(),
   	func := o.element.lambda(),
   	t := TempVar(x.t.t),
  	chain(assign(t, 0), loop(i, func.domain(), assign(t, add(t, mul(func.at(i), nth(x, i))))), assign(nth(y, 0), t)))
);

#DefaultCodegen.RowVec2 := (self, o, y, x, opts) >>

#	v := var.fresh_t("v", TArray(TInt, o.element.domain())),
#	t := TempVar(x.t.t),
#	decl([v,t], chain(assign(t, 0),
#	loopw(lt(t, o.element.domain()),
#	chain(if1(logic_and(neq(nth(v, i), V(0)), neq(nth(x,i), V(0))), assign(t, add(t, mul(nth(v,i), nth(x,i))))) , assign(nth(y,0), t)))))),
#	let(i := Ind(),
#   func := o.element.lambda(),
#   	t := TempVar(x.t.t),
#  	chain(assign(t, 0), loop(i, func.domain(), assign(t, add(t, mul(func.at(i), nth(x, i))))), assign(nth(y, 0), t))));

DefaultCodegen.TSparse_Matrix := meth(self, o, y, x, opts)
    local i, j, k, n, b;
    i := Ind();
    j := Ind();
    k := Ind();
    n := var.fresh_t("n", TInt);
    b := var.fresh_t("b", TPtr(TInt));
    return decl([], chain(
            loopf(i, V(0), n, decl([], chain(
                    loopf(j, deref(add(x,i)), deref(add(x, add(i, V(1)))), decl([], chain(
                            loopf(k, nth(b, deref(add(x, add(n, add(V(1),j))))), nth(b, add(deref(add(x, add(n, add(V(1),j))),V(1)))), decl([], chain(
                                assign(y, V(0))
                                #assign(nth(y, add(nth(b,add(n,add(V(1),k))), mul(i, n))), V(0)),
                                #assign(nth(y, add(nth(b,add(n,add(V(1),k))), mul(i, n))), add(nth(y, add(nth(b,add(n,add(V(1),k))), mul(i, n))), mul(deref(add(x, add(n, add(nth(x,n), add(V(1), j))))), deref(add(b, add(n, add(nth(b,n), add(V(1), k))))))))
                            )))
                    )))
            )))
    ));
    end;

DefaultCodegen.Ttrace := meth(self, o, y, x, opts)
	local i, j, n;
	i := Ind();
	j := Ind();
	n := var.fresh_t("n", TInt);
	return decl([], chain(
		loopf(i, V(0), sub(n,V(1)), decl([], chain(
			loopf(j, deref(add(x, i)), deref(add(x, add(i, V(1)))), decl([], chain(
				if1(eq(i, deref(add(x,add(n,j)))), assign(y, add(deref(y), deref(add(x, add(n, add(nth(x,n), add(V(1), j))))))))
			)))
		)))
	));
	end;

DefaultCodegen.VStack := meth(self, o, y, x, opts)
	local iA, jA, v, i, j, n;
	iA := var.fresh_t("iA", TPtr(TInt));
	jA := var.fresh_t("jA", TPtr(TInt));
	#val := var.fresh_t("val", TPtr(TInt));
	v := var.fresh_t("v", TPtr(TInt));
	i := Ind();
	j := Ind();
	n := var.fresh_t("n", TInt);
	return decl([iA, jA], chain(
		assign(iA, x),
		assign(n, V(10)),
		assign(jA, add(x, add(n, V(1)))),
		#assign(val, add(IJ), add(n, add(1, nth(iA, n))))
		loopf(i, V(0), n, decl([v], chain(
			loopf(j, deref(add(iA, i)), deref(add(iA, add(i,V(1)))), chain(
				#assign(deref(add(y,i)), V(0))
				assign(deref(add(y,i)), add(deref(add(y,i)), mul(deref(add(jA, add(tcast(TInt, deref(add(iA,n))),j))), deref(add(v, deref(add(jA, j)))))))
				#assign(deref(add(y,i)), deref(add(y,i)), mul(deref(add(jA, add(nth(iA,n), j))), deref(add(v, deref(add(jA, j)))))))
			)))
		))
	));
	end;

#DefaultCodegen.Gath := meth(self, o, y, x, opts)
#	local i, index_j;
#	i := Ind();
#	index_j := o.func._children[1].params[2];
#	return loop(i, o.func.domain(), assign(nth(y, i), nth(x, add(opts.symbol[1], add(nth(x,opts.symbol[1]), add(index_j, V(1)))))));
#end;
#DefaultCodegen.Scat := meth(self, o, y, x, opts)
#	local i, index_j;
#	i := Ind();
#	index_j := o.func._children[1].params[2];
#	return loop(i, o.func.domain(), assign(nth(y,nth(X, add(opts.symbol[1], add(V(1), index_j)))), nth(x, i)));
#end;

DefaultCodegen.PreDstruct := (self, o, y, x, opts) >>
	Cond(IsBound(opts.init), opts.init(opts.symbol[2], opts.symbol[3], o.count), Error("Wrong opts\n"));

DefaultCodegen.PostDstruct := (self, o, y, x, opts) >>
	Cond(IsBound(opts.to_csr), opts.to_csr(opts.symbol[2]), Error("Wrong opts\n"));

DefaultCodegen.DStructKernelAxpy := (self, o, y, x, opts) >> Cond(IsBound(opts.init) and IsBound(opts.to_csr), let(curr := var.fresh_t("curr", TPtr(Tdstruct)),
													c1 := decl([curr], chain(opts.init(opts.symbol[2], opts.symbol[3], opts.symbol[1]), assign(curr, opts.symbol[2]))),
													c3 := opts.to_csr(y, opts.symbol[2]), c2 := self(o.child(1), curr, x, opts), newcmd := chain(c2.cmd, assign(curr, next_ptr(curr, "next"))),
													newc := loopn(c2.var, c2.range, newcmd), chain(c1,newc, c3)),  Error("Wrong opts\n"));

DefaultCodegen.DStructKernelDot := (self, o, y, x, opts) >> Cond(IsBound(opts.init) and IsBound(opts.to_csr), let(temp := var.fresh_t("temp", opts.YType), curr := var.fresh_t("curr", TPtr(Tdstruct)),
													c1 := decl([curr], chain(opts.init(opts.symbol[2], opts.symbol[3], opts.symbol[1]), assign(curr, opts.symbol[2]))),
													c3 := opts.to_csr(y, opts.symbol[2]), c2 := self(o.child(1), temp, x, opts),
													beg := Collect(c2, @(1,nth, e->e.loc = temp))[1], val := Collect(c2, @(1,mul, e->IsNth(e.args[1]) and e.args[1].loc = x))[1],
													result := SubstTopDown(c2, @(1,assign, e->e.loc = beg), e-> insert_ds(curr, beg.idx, val, opts.symbol[3])),
													newcmd := chain(result.cmd, assign(curr, next_ptr(curr, "next"))),
													newc := loopn(c2.var, c2.range, newcmd), chain(c1,newc, c3)),  Error("Wrong opts\n"));

DefaultCodegen.GrpKernel := (self, o, y, x, opts) >> Cond(opts.vec = true,
									Cond(Length(Collect(o, TSparse)) > 0,
										let(temp1 := var.fresh_t("temp", opts.YType), temp2 := var.fresh_t("temp", opts.XType), ring := Collect(o, TSparse)[1].ring,
									c := self(o.child(1), temp1, temp2, opts), new_val := o.child(1)._children[1].var,
									result := SubstTopDown(c, @(1,nth, e->e.loc = temp2), e-> nth(x, add(opts.symbol[1], add(nth(x,opts.symbol[1]), add(new_val, V(1)))))),
									result2_5 := SubstTopDown(result, @(1, assign_acc, e->IsNth(e.loc) and e.loc.loc = temp1), e->assign(e.loc, ring.sum(e.loc, e.exp))),
									result2 := SubstTopDown(result2_5, @(1,nth, e->e.loc = temp1), e-> nth(y,nth(x, add(opts.symbol[1], add(V(1), new_val))))), result2),
										let(temp1 := var.fresh_t("temp", opts.YType), temp2 := var.fresh_t("temp", opts.XType), ring := Collect(o, fConst)[1].params[3].loc.t.ring,
									c := self(o.child(1), temp1, temp2, opts), new_val := o.child(1)._children[1].var,
									result := SubstTopDown(c, @(1,nth, e->e.loc = temp2), e-> nth(x, add(opts.symbol[1], add(nth(x,opts.symbol[1]), add(new_val, V(1)))))),
									result2_5 := SubstTopDown(result, @(1, assign_acc, e->IsNth(e.loc) and e.loc.loc = temp1), e->assign(e.loc, ring.sum(e.loc, e.exp))),
									result2 := SubstTopDown(result2_5, @(1,nth, e->e.loc = temp1), e-> nth(y,nth(x, add(opts.symbol[1], add(V(1), new_val))))), result2)),
									opts.mat = true,
									let(temp1 := var.fresh_t("temp", opts.YType), temp2 := var.fresh_t("temp", opts.XType), ring := Collect(o, TSparse)[1].ring,
									c := Cond(IsBound(opts.ds) and opts.ds = true, self(o.child(1), y, temp2, opts), self(o.child(1), temp1, temp2, opts)), sa := Collect(o, ScatAcc),
									result := SubstTopDown(c, @(1,nth, e->e.loc = temp2), e-> opts.sparse_mat.index_val(x, opts.symbol[1], sa[1].func._children[1].params[2])),
									result2_5 := SubstTopDown(result, @(1, assign_acc, e->IsNth(e.loc) and e.loc.loc = temp1), e->assign(e.loc, ring.sum(e.loc, e.exp))),
									result2 := SubstTopDown(result2_5, @(1,nth, e->e.loc = temp1), e-> nth(y, add(mul(opts.sparse_mat.index_row(x, opts.symbol[1], sa[1].func._children[1].params[2]), opts.symbol[1]), o.child(1).var))), result2),
									self(o.child(1), y, x, opts));

DefaultCodegen.GrpKernel2 := (self, o, y, x, opts) >> let(ring := Collect(o, TSparse)[1].ring, temp1 := var.fresh_t("temp", opts.XType), temp2 := var.fresh_t("temp", opts.XType), c := self(o.child(1), temp1, temp2, opts),
							idx := var.fresh_t("idx", TInt), mid := var.fresh_t("mid", opts.XType), lv := Ind(),
							result := SubstTopDown(c, @(1,nth, e->e.loc = temp2 and Length(Collect(e, opts.symbol[1])) = 0 and Length(Collect(e, @(1,var, g-> g <> temp2))) > 1), e-> nth(x, add(opts.symbol[1], add(nth(x,opts.symbol[1]), add(o.child(1)._children[1]._children[1].var, V(1)))))),
							resultt := SubstTopDown(result, @(1, var, e->e = temp2), e->x),
							varr := Cond(Length(Collect(resultt, loopn)) = 2, Collect(resultt, @(1, assign, e-> IsNth(e.exp) and Length(Collect(e, add)) = 2 and e.exp.loc = x))[1].loc, o.child(1)._children[1]._children[1].var),
							result2 := Cond(Length(Collect(resultt, loopn)) = 2, SubstTopDown(resultt, @(1,nth, e->e.loc = temp1), e-> nth(temp1,varr)), SubstTopDown(resultt, @(1,nth, e->e.loc = temp1), e-> nth(temp1,opts.sparse_mat.index_row(x, opts.symbol[1], varr)))),
							result2_5 := SubstTopDown(result2, @(1, assign_acc, e-> IsNth(e.loc) and e.loc.loc = temp1), e->assign(e.loc, ring.sum(e.loc, e.exp))),
							c2 := Cond(Length(Collect(resultt, loopn)) = 2, chain(assign(nth(mid, varr), V(1))),  chain(assign(nth(mid, opts.sparse_mat.index_row(x, opts.symbol[1], varr)), V(1)))),
							result3 := SubstTopDown(result2_5, @(1, chain, e-> Length(Collect(e, assign)) >= 2 and Length(Collect(e, assign)) < 4), e-> chain(e.cmds, c2.cmds)),
							decl([mid], chain(assign(mid, malloc(mid.t.t, opts.symbol[1])), loopf(lv, V(0), opts.symbol[1], assign(nth(mid, lv), V(0))), loopf(c.var, V(0), c.range, decl([temp1], chain(assign(temp1, malloc(temp1.t.t, opts.symbol[1])), result3.cmd, chain(add_to_sparse(addrof(y), mid, temp1, opts.symbol[1], c.var), loopf(lv, V(0), opts.symbol[1], assign(nth(mid, lv), V(0))))))))));

DefaultCodegen.MaskKernelAxpy := (self, o, y, x, opts) >> let(itr1 := var.fresh_t("itr", TInt), itr2 := var.fresh_t("itr", TInt), col1 := var.fresh_t("xcol", TInt),
		col2 := var.fresh_t("spmcol", TInt), m2 := o.child(1)._children[1].element.var, Append(opts.symbol, [m2]), i := o.child(1)._children[2].func._children[1].params[1], x2 := o.child(1)._children[3].element.params[3].ofs.loc,
		j := o.child(1)._children[3].element.params[3].ofs.idx, decl([itr1, itr2], chain(assign(itr1, nth(m2,i)), assign(itr2, nth(x, nth(x2, add(V(1), add(j, opts.symbol[1]))))), loopw(logic_and(neq(itr1, nth(m2, add(i,V(1)))), neq(itr2,nth(x, nth(x2, add(add(V(1), add(j, opts.symbol[1])),V(1)))))), decl([col1, col2], chain(assign(col1, opts.sparse_mat.index_row(m2,opts.symbol[1], itr1)),
		assign(col2,opts.sparse_mat.index_row(x,opts.symbol[1], itr2)), if3(gt(col1, col2), assign(itr2, add(itr2, V(1))), lt(col1, col2), assign(itr1, add(itr1, V(1))),
		chain(assign(nth(y, add(mul(opts.sparse_mat.index_row(x, opts.symbol[1], itr2), opts.symbol[1]),i)), add(nth(y, add(mul(opts.sparse_mat.index_row(x, opts.symbol[1], itr2), opts.symbol[1]), i)), mul(opts.sparse_mat.index_val(x, opts.symbol[1], itr2), opts.sparse_mat.index_val(x2, opts.symbol[1], j)))),
		assign(itr1, add(itr1, V(1))), assign(itr2, add(itr2, V(1)))))))))));

DefaultCodegen.MaskKernelDot := (self, o, y, x, opts) >> let(c := self(o.child(1)._children[2], y, x, opts), mitr := var.fresh_t("mitr", TInt), m2 := o.child(1)._children[1].element.var, Append(opts.symbol, [m2]), i := o.child(1)._children[1].element.ofs.args[1],
			j := o.child(1)._children[1].element.ofs.args[2], decl([mitr], chain(assign(mitr, nth(m2,i)), loopw(logic_and(lt(mitr, nth(m2, add(i, V(1)))), lt(opts.sparse_mat.index_row(m2, opts.symbol[1], mitr), j)),
			chain(assign(mitr, add(mitr, V(1))))), if1(logic_and(eq(opts.sparse_mat.index_row(m2, opts.symbol[1], mitr), j), eq(opts.sparse_mat.index_val(m2, opts.symbol[1], mitr), V(1))), c))));

DefaultCodegen.AccumKernel := (self, o, y, x, opts) >> let(ring := Collect(o, TSparse)[1].ring, i := Ind(), j := Ind(), accum := var.fresh_t("accum", TPtr(TInt)), create := decl([accum], chain(assign(accum, malloc(accum.t.t, opts.symbol[1])),
	loopf(i, V(0), opts.symbol[1], loopf(j, V(0), opts.symbol[1], assign(nth(accum, add(mul(i,opts.symbol[1]), j)), ring.zero()))))), c:= self(o.child(1), accum, x, opts),
	copy := decl([], chain(loopf(i, V(0), opts.symbol[1], loopf(j, V(0), opts.symbol[1], assign_acc(nth(y, add(mul(i*opts.symbol[1]), j)), nth(accum, add(mul(i*opts.symbol[1]), j))))))), res := chain(create, c, copy), res);

DefaultCodegen.TriCountKernel := (self, o, y, x, opts) >> let(tempalg := var.fresh_t("tempalg", opts.YType), c:= self(o.child(1), tempalg, x, opts), v2 := o.child(1)._children[1].var, ring := Collect(o, TSparse)[1].ring, itr1 := Collect(c, @(1, var, e->e.id = "itr1"))[1],itr2 := Collect(c, @(1, var, e->e.id = "itr2"))[1],
	xcol1 := Collect(c, @(1, var, e->e.id = "xcol1"))[1], spmcol1 := Collect(c, @(1, var, e->e.id = "spmcol1"))[1], res := SubstTopDown(c, @(1, assign, e-> e.loc = xcol1 and IsNth(e.exp)), e->skip()),
	res2 := SubstTopDown(res, @(1, assign, e-> e.loc = spmcol1 and IsNth(e.exp)), e->skip()), result := SubstVars(res2, rec((xcol1.id) := var.table.itr1)), result2 := SubstVars(result, rec((spmcol1.id) := var.table.itr2)),
	s := var.fresh_t("start", opts.XType), ee := var.fresh_t("end", opts.XType), se := chain(assign(s, add(x, add(opts.symbol[1], add(V(1), add(nth(x, o.child(1).var)))))), assign(ee, add(x, add(opts.symbol[1], add(V(1), add(nth(x, add(o.child(1).var, V(1))))))))),
	result2_5 := SubstTopDown(result2, @(1, loopn, e->e.range=opts.symbol[1]), e-> loopn(e.var, e.range, decl([s,ee],chain(se, e.cmd)))), itr := var.fresh_t("i", TPtr(TInt)),
	result2_6 := SubstTopDown(result2_5, @(1, loopn, e->e.range<>opts.symbol[1]), e-> loopn(itr, sub(ee, s), e.cmd)), result2_7 := SubstBottomUp(result2_6,@(1,var, e->e.id=v2.id), e->itr),
	result3 := SubstTopDown(result2_7, @(1, assign, e-> e.loc = itr2 and Length(Collect(e.exp,itr2)) = 0), e->assign(itr2, add(X, add(1, add(opts.symbol[1], nth(x,deref(itr))))))),
	result3_5 := SubstTopDown(result3, @(1, assign, e-> e.loc = itr1 and Length(Collect(e.exp,itr1)) = 0), e->assign(itr1, s)),
	en2 := add(x, add(V(1), add(opts.symbol[1], nth(x,add(deref(itr), V(1)))))),
	result4_1:= SubstTopDown(result3_5, @(1, gt), e-> gt(deref(e.args[1]), deref(e.args[2]))), result4_2:= SubstTopDown(result4_1, @(1, lt), e-> lt(deref(e.args[1]), deref(e.args[2]))),
	result4 := SubstTopDown(result4_2, @(1, loopw, e->Length(Collect(e.cond, neq)) > 0), e-> loopw(logic_and(lt(itr1, ee), lt(itr2, en2)), e.cmd.cmd)),
	result5 := SubstTopDown(result4, @(1, assign, e-> IsNth(e.loc) and e.loc.loc = tempalg), e-> assign_acc(deref(y), V(1))), result5);

	#let(c:= self(o.child(1), y, x, opts), ring := Collect(o, TSparse)[1].ring, itr1 := Collect(c, @(1, var, e->e.id = "itr1"))[1],itr2 := Collect(c, @(1, var, e->e.id = "itr2"))[1],
	#out1 := var.fresh_t("out", itr1.t), out2 := var.fresh_t("out", itr1.t), res := SubstTopDown(c, @(1, assign, e-> e.loc = itr1 and IsNth(e.exp)), e->assign(out1, e.exp)), v2 := o.child(1)._children[1].var,
	#res2 := SubstTopDown(res, @(1, assign, e-> e.loc = itr2 and IsNth(e.exp)), e->assign(out2, e.exp)),result := SubstVars(res2, rec((itr1.id) := var.table.xcol1)), result2 := SubstVars(result, rec((itr2.id) := var.table.spmcol1)),
	#result3 := SubstTopDown(result2, @(1, assign, e->IsBound(e.loc.id) and e.loc.id = "xcol1" and IsNth(e.exp)), e-> assign(e.loc, v2)), result3); #if3chain := chain(Collect(result2, if3)[1]),result3 := SubstTopDown(result2, @(1, chain, e-> Length(Collect(e, if3)) = 1 and Length(Collect(e, loopw)) = 0), e-> if3chain), result3);,

DefaultCodegen.PushBFSKernel := (self, o, y, x, opts) >> let(front := Filtered(Collect(o, var), x-> x = opts.symbol[2])[1], c:= self(o.child(1), front, x, opts), expr := Collect(c, @(1,nth, e-> e.loc = x and Length(Collect(e.idx, add)) = 2))[1],
			if_cond := if1(logic_and(eq(nth(front, expr), V(1)), eq(nth(y, expr), add(opts.symbol[1], V(1)))), assign(nth(y, expr), c.var)),
			result := SubstTopDown(c, @(1, chain, e-> Length(Collect(e, assign)) = 3), e-> chain(e.cmds, if_cond)), result);
			#result1 := SubstTopDown(c, @(1, assign, e-> IsNth(e.loc) and e.loc.loc = y and Length(Collect(e.exp, logic_or)) > 0), e-> SubstVars(e, rec(y := front))),

DefaultCodegen.PullBFSKernel := (self, o, y, x, opts) >> let(c:= self(o.child(1), y, x, opts), front := Collect(o, FDataOfs)[1].var,
			result1 := SubstTopDown(c, @(1, loopn, e-> Length(Collect(e.range, sub)) = 0), e->loopn(e.var, e.range, if1(eq(nth(y, c.var), add(opts.symbol[1], V(1))), e.cmd))),
			result2 := SubstTopDown(result1, @(1, loopn, e1-> Length(Collect(e1.range, sub)) > 0), e-> loopn(e.var, e.range, if1(neq(nth(y, opts.sparse_mat.index_row(x, opts.symbol[1], e.var)), add(opts.symbol[1], V(1))), chain(assign(nth(front,opts.sparse_mat.index_row(x, opts.symbol[1], e.var)),V(1)), assign(nth(y, c.var), opts.sparse_mat.index_row(x, opts.symbol[1], e.var)), break())))),
			result2);

DefaultCodegen.BFSKernel := (self, o, y, x, opts) >> let(c:= self(o.child(1), y, x, opts), source := var.fresh_t("source", TInt), Append(opts.symbol, [source]), front := Filtered(Collect(o, var), x-> x = opts.symbol[2])[1], i := Ind(), sum := var.fresh_t("sum", TInt), rsum := var.fresh_t("rsum", TInt),
			decorator := decl([sum, rsum], chain(assign(nth(y,source), source), assign(nth(front,source), V(1)), assign(sum, V(1)), assign(rsum, V(0)), loopw(neq(sum, rsum), chain(assign(sum, rsum), assign(rsum, V(0)),c, loopf(i, V(0), opts.symbol[1], chain(assign_acc(rsum, nth(front, i)))))))), decorator);
#DefaultCodegen.AccumKernel := (self, o, y, x, opts) >> let(ring := Collect(o, TSparse)[1].ring, i := Ind(), j := Ind(), accum := var.fresh_t("accum", TPtr(TInt)), create := decl([accum], chain(assign(accum, malloc(accum.t.t, opts.symbol[1])),
#	loopf(i, V(0), opts.symbol[1], loopf(j, V(0), opts.symbol[1], assign(nth(accum, add(mul(i,opts.symbol[1]), j)), ring.zero()))), self(o.child(1), accum, x, opts), decl([], chain(loopf(i, V(0), opts.symbol[1], loopf(j, V(0), opts.symbol[1],
#	assign_acc(nth(y, add(mul(i,opts.symbol[1]), j)), nth(accum, add(mul(i,opts.symbol[1]), j))))))))), create);

#DefaultCodegen.Gath := meth(self, o, y, x, opts)
#	local i, index_i, index_j;
#	i := Ind();
#	index_j := o.func._children[2].params[2];
#	return decl([], chain(
#		loopf(i, nth(x, index_j), nth(x, add(index_j, V(1))), chain(
#			assign(nth(y, sub(i, nth(x, index_j))), nth(x, add(opts.symbol[1], add(nth(x,opts.symbol[1]), add(i, V(1))))))
#		))
#	));
#end;
#
#DefaultCodegen.Scat := meth(self, o, y, x, opts)
#	local i, index_i, index_j;
#	i := Ind();
#	index_j := o.func._children[2].params[2];
#	return decl([], chain(
#		loopf(i, nth(x, index_j), nth(x, add(index_j, V(1))), chain(
#			assign(nth(y,nth(x, add(opts.symbol[1], add(V(1), i)))), nth(x, sub(i, nth(x, index_j))))
#		))
#	));
#end;
#
#DefaultCodegen.RowVec2 := meth(self, o, y, x, opts)
#	local body, result, sa;
#	sa := var.fresh_t("sa", TArray(TInt, o.element.object.size));
#	result := var.fresh_t("result", TInt);
#	body := assign(y , itr);
#	return decl([result, sa], chain(
#		assign(result, V(0)),
#		o.traversal(sa, body)
#	));
#	end;
#
DefaultCodegen.NewScalarProduct := meth(self, o, y, x, opts)
	local itr, num_nz;
	itr := o.rv.element.object.get_iterator();
	num_nz := o.rv.element.object.num_nonzeros("sa2");
	#return decl([itr], chain(
	return decl([], chain(
		assign(y, V(0)),
		assign(itr, V(0)),
		loopw(lt(itr, o.rv.element.object.size), chain(
		if1(logic_and(neq(struct_nth(x, "value", itr), V(0)),neq(struct_nth("sa2", "value", itr), V(0))), chain(
			assign(y, add(deref(y), mul(struct_nth(x, "value", itr), struct_nth("sa2", "value", itr)))),
			assign(itr, add(itr, V(1)))))))));
	end;

#if1(neq(struct_nth(x, "value", itr), V(0)), chain(
#assign(y, add(deref(y), mul(

DefaultCodegen.DenseScalarProduct := meth(self, o, y, x, opts)
	local arr, itr, size;
	arr := var.fresh_t("arr", TArray(o.rv.element.object.t.t, o.rv.element.object.size));
	itr := var.fresh_t("itr", TInt);
	return data(arr, V(o.rv.element.tolist()), decl([], chain(
		assign(y, V(0)),
		assign(itr, V(0)),
		loopw(lt(itr, o.rv.element.object.size),
		chain(
			if1(neq(nth(arr, itr), V(0)), chain(
			assign(y, add(deref(y), mul(nth(arr, itr), nth(x, itr)))),
			assign(itr, add(itr, V(1))))))))));
	end;

DefaultCodegen.SparseScalarProduct := meth(self, o, y, x, opts)
	local arr, itr, size;
	arr := var.fresh_t("arr", TArray(o.rv.element.object.t.t, o.rv.element.object.size));
	itr := var.fresh_t("itr", TInt);
	return data(arr, V(o.rv.element.tolist()), decl([], chain(
		assign(y, V(0)),
		assign(itr, V(0)),
		loopw(lt(itr, o.rv.element.object.size),
		chain(
			if1(logic_and(neq(nth(arr, itr), V(0)),neq(nth(x, itr), V(0))), chain(
			assign(y, add(deref(y), mul(nth(arr, itr), nth(x, itr)))),
			assign(itr, add(itr, V(1))))))))));
	end;

DefaultCodegen.TSparse := meth(self, o, y, x, opts)
	local row, values;
	row := var.fresh_t("row", TArray(o.t, o.size1));
	values := var.fresh_t("values", TArray(o.t, o.size2));
	return data(row, o.list1, skip());
	end;

DefaultCodegen.sparse_nth2 := meth(self, o, y, x, opts)
	local row, index;
	row := var.fresh_t("row", TArray(o.sa.t, o.sa.size));
	index := var.fresh_t("index", TInt);
	return decl([row, index], chain(
	if4(eq(deref(row + index), V(0)), assign(y, V(0)), assign(y, deref(row+index)))));
	end;

DefaultCodegen.Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, destroysub, io, t, initcode, initparams;

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

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

        #o :=  Process_fPrecompute(o, opts);

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

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);
        #[icode,t] := UTimedAction(ESReduce(icode, opts)); #PrintLine("ESReduce ", t);
        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        #PrintLine("---compile--");
        #DoForAll([1..Length(Compile.times)], i -> PrintLine(i, " ", Compile.times[i], " ",
        #        let(f:=opts.compileStrategy[i], When(IsFunc(f) or IsMeth(f), "---", f))));

        # icode := PowerOpt(icode);
        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        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;

        io := When(x=y, [x], [y, x]);
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        destroysub := Cond(IsBound(opts.subName), Concat("destroy_", opts.subName), "destroy");
        icode := func(TVoid, sub, Concatenation(io, params), Compile(icode, opts));

        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
	    initcode := chain(List(Filtered(datas, e -> IsBound(e.var.init)), x -> SReduce(x.var.init, opts)));
            prog := program(
                decl(List(datas, x->x.var),
                    chain(
                        func(TVoid, initsub, initparams :: Set(Collect(initcode, param)), initcode),
                        icode,
                        func(TVoid, destroysub, [], skip())
                    )));
        else
            prog := program( func(TVoid, initsub, params, chain()), icode);
        fi;
        prog.dimensions := o.dims();
        return prog;
    end;

#CUnparser.TSparse := (self, o, i, is) >> Print(Blanks(i),
#	o.t, "row[", o.size1, "];\n", o.t, "values[", o.size2, "];\n");

SparseDefaults := CopyFields(SpiralDefaults, rec(
  compileStrategy := GraphIndicesCS,
  X := var("X", TPtr(TInt)),
  XType := TPtr(TInt),
  arrayDataModifier := "",
  arrayBufModifier := "",
  Y := var("Y", TPtr(TInt)),
  matrix := TSparse_Matrix(TSparse(TArray(TInt, 5), TSemiring_Arithmetic(TInt)), []),
  YType := TPtr(TInt),
#  isCSR := true,
#  includes := ["<sparse.h>"],
#  #symbol := ["struct_array sa2"]
));

Class(SparseOpts, SpiralDefaults, rec(
    tags := [],
    operations := rec(Print := s -> Print("<Sparse options record>")),
    tagIt := (self, t) >> t.withTags(self.tags),
	name := "SparseOpts",
	arrayDataModifier := "",
  	arrayBufModifier := "",
	vec := false,
	mat := false,
	sparse_mat := TSparse_Matrix(TSparse(TArray(TInt, 5), TSemiring_Arithmetic(TInt)), []),
	sparse_vec := TSparse(TArray(TInt, 5), TSemiring_Arithmetic(TInt)),
	#createRuleTree := (self, t) >> let(t2 := ApplyStrategy(t, self.rewrite, UntilDone, self), RandomRuleTree(t2, self)),
    search := (self, t) >> RandomRuleTree(t, self),
    sumsRuleTree := (self, rt) >> let(spl := SPLRuleTree(rt), spl2 := ApplyStrategy(spl, self.rewrite, UntilDone, self), rt2 := RandomRuleTree(spl2, self),
									 srt := ApplyStrategy(SumsRuleTree(rt2, self), self.rewrite, UntilDone, self), srt),
    codeSums := meth(self, ss)
        local c, X, Y, opts;
        opts := self;
        X := var("X", self.XType);
        Y := var("Y", self.YType);
		#if self.vec = true then
		#	self.symbol := Concat(self.symbol, [self.sparse_vec.get_var()]);
		#fi;
		#if self.mat = true then
		#	self.symbol := Concat(self.symbol, [self.sparse_mat.get_var()]);
		#fi;
		#if IsBound(self.algo) and self.algo = true then
		#	self.symbol := Sublist(self.symbol, [1..Length(self.symbol)-1]);
		#fi;
        c := opts.codegen(Formula(ss), Y, X, opts);
        c.cmds := Cond(IsBound(opts.postProcessCode), opts.postProcessCode(c.cmds), c.cmds);
        c.ruletree := Cond(IsBound(ss.ruletree), ss.ruletree, rec());
		if self.vec = true then
			c := SubstTopDown(Copy(c), @(1, var, e-> IsBound(e.t.t) and e.t.t = self.symbol[Length(self.symbol)].t.t), e->self.symbol[Length(self.symbol)]);
		fi;
		#if self.mat = true then
		#	c := SubstTopDown(Copy(c), @(1, var, e->IsBound(e.t.t) and e.t.t = self.sparse_mat.element.t), e->self.symbol[Length(self.symbol)]);
		#fi;
        if IsBound(self.params) then
            c := SubstBottomUp(c, @(1, func, e->e.id = "transform"),
                e ->func(@(1).val.ret, @(1).val.id, @(1).val.params::self.params, @(1).val.cmd));
        fi;
        #if IsBound(opts.useBinSplit) and opts.useBinSplit then c:= BinSplit(c); fi;
        return c;
    end,
    #prettyPrint := (self, c) >> PrintCode(c.ruletree.node.params[2].fname, c, self),
	prettyPrint := (self, c) >> PrintCode(self.name, c, self),
    genSparse := (self, t) >> self.codeSums(self.sumsRuleTree(self.search(t))),
    globalUnrolling := 0
));

BaseIndicesCS := [
    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
];
#
#SparseIndicesCS := Concatenation(BaseIndicesCS, [
#    c -> MarkDefUse(c),
#    (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))),
#    c -> MarkDefUse(c),
#    (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))),
#    (c, opts) -> Cond(opts.finalBinSplit, BinSplit(c, opts), c),
#    (c, opts) -> Cond(IsBound(opts.scheduler), opts.scheduler(c, opts), c),
#    c -> FixValueTypes(c),
#    c -> Compile.declareVars(c)
#]);

Declare(SparseDefaultConf);
Declare(SparseDataStructConf);
Declare(SparseStructOutConf);

Class(SparseDefaultConf, LocalConfig, rec(
    __call__ :=  arg -> CopyFields(SparseDefaultConf,
        rec(
        )),
    getOpts := meth(self, t)
        local opts, n;
        opts := Copy(SparseOpts);
		opts.compileStrategy := [c -> Compile.pullDataDeclsRefs(c), c -> Compile.fastScalarize(c), c -> FlattenCode(c),
    	c -> UntangleChain(c), c -> MarkDefUse(c), c -> FixValueTypes(c), c -> Compile.declareVars(c)];
		opts.rewrite := [];
        opts.XType := TPtr(TInt);
        opts.YType := TPtr(TInt);
        opts.useDeref := false;
		n := var.fresh_t("n", TInt);
		opts.symbol := [n];
		opts.globalUnrolling := 0;
        return opts;
    end,
    operations := rec(Print := s -> Print("<Sparse Default Configuration>")),
));

Class(SparseStructOutConf, LocalConfig, rec(
    __call__ :=  arg -> CopyFields(SparseDefaultConf,
        rec(
        )),
    getOpts := meth(self, t)
        local opts, n;
        opts := Copy(SparseOpts);
		opts.compileStrategy := [c -> Compile.pullDataDeclsRefs(c), c -> Compile.fastScalarize(c), c -> FlattenCode(c),
    	c -> UntangleChain(c), c -> MarkDefUse(c), c -> FixValueTypes(c), c -> Compile.declareVars(c)];
		opts.rewrite := [];
        opts.XType := TPtr(TInt);
        opts.YType := TPtr(Tcsr);
        opts.useDeref := false;
		n := var.fresh_t("n", TInt);
		opts.symbol := [n];
		opts.globalUnrolling := 0;
		opts.postProcessCode := (self, c) >> let(t1 := var.fresh_t("iA", TPtr(TInt)), t2 := var.fresh_t("jA", TPtr(TInt)),
			t3 := var.fresh_t("val", TPtr(TInt)), n := var.fresh_t("n", TInt), ts := T_Struct("csr", [t1,t2,t3,n]), Append(c, [define([ts])]), Reversed(c));
        return opts;
    end,
    operations := rec(Print := s -> Print("<Sparse Struct Out Configuration>")),
));

Class(SparseDataStructConf, SparseDefaultConf, rec(
	__call__ :=  arg -> CopyFields(SparseDefaultConf,
        rec(
        )),
		getOpts := meth(self, t)
        local opts, n, ds, root;
        opts := Copy(SparseOpts);
        opts.XType := TPtr(TInt);
        opts.YType := TPtr(TPtr(TInt));
        opts.useDeref := false;
		opts.compileStrategy := [c -> Compile.pullDataDeclsRefs(c), c -> Compile.fastScalarize(c), c -> FlattenCode(c),
    	c -> UntangleChain(c), c -> MarkDefUse(c), c -> FixValueTypes(c), c -> Compile.declareVars(c)]; #No copy prop support for linked lists
		opts.rewrite := [];
		opts.ds := true;
		opts.init := (self, ds, tree, value) >> let(itr := var.fresh_t("itr", TPtr(Tdstruct)),
		 decl([itr], chain(assign(ds, allocate_ds(value)), assign(itr, ds), assign(tree, dstructToTree(addrof(itr), value)))));
		opts.to_csr := (self, out, ds) >> chain(assign(deref(out), ds_to_csr(ds)));
		opts.destory_ds := (self, ds, root) >> chain(destory_ds(ds, root));
		n := var.fresh_t("n", TInt);
		ds := var.fresh_t("d", TPtr(Tdstruct));
		root := var.fresh_t("root", TPtr(Ttree));
		opts.symbol := [n, ds, root];
        opts.includes := ["\"data_structure.h\""];
		opts.globalUnrolling := 0;
		opts.codegen := CopyFields(SpiralDefaults.codegen, DataStructureCodegenMixin);
        return opts;
    end,
	operations := rec(Print := s -> Print("<Sparse DataStruct Configuration>")),
));


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

#F NthRootCyc(<cyc>, <nroot>) - <nroot>-th root of <cyc>
#F   to get a root we take all powers of cyc of order n,
#F   and make them powers of cyc of order nroot*n (effectively dividing these
#F   powers by nroot). This is done by appending n*(nroot-1) zeros to list of
#F   coeffs and using CoeefsCyc.
NthRootCyc := (cyc, nroot) -> let(n := OrderCyc(cyc),
    coeffs := CoeffsCyc(cyc,n), # -1 must be E(2), not -E(1)
    sum := Sum(coeffs),
    Cond( # GAP normalization of cyc's (to [-j,j], and not [1,-1]) give headache here
	sum =  1, CycList(Concatenation(coeffs, Replicate(n*(nroot-1), 0))),
	sum = -1, CycList(Concatenation(Replicate(n/2, 0), -coeffs{[1..n/2]}, Replicate(n*(nroot-1), 0))),
	Error("E(N) or E(-N) where N is positive integer is expected")));

#F SkewDFT(n, alpha, k) - Fourier transform for algebra C[x]/x^n - omega_1^alpha
#F   k - rotation
#F   matrix [ w_n ^ (r+alpha)c ]_{r,c}
#F
Class(SkewDFT, TaggedNonTerminal, rec(
    abbrevs := [ (n)       -> Checked(IsPosIntSym(n), [n, 1, 1]),
	         (n,alpha) -> Checked(IsPosIntSym(n), IsRatSym(alpha), [n, alpha, 1]),
	         (n,alpha,k) -> Checked(IsPosIntSym(n), IsIntSym(k), AnySyms(n,k) or Gcd(_unwrap(n),_unwrap(k))=1, IsRatSym(alpha), [n, alpha, k]) ],

    dims := self >> [ self.params[1], self.params[1] ],
    terminate := self >> let(
	N := self.params[1], rot := self.params[3], a := self.params[2], j := Ind(N-1),
        mat := DFT(N).terminate() *
               Diag(diagDirsum(fConst(TReal, 1, 1.0), fPrecompute(
                      Lambda(j, cospi(fdiv(2*a*rot*(j+1), N)) + E(4)*sinpi(fdiv(2*a*rot*(j+1), N)))))),
        Cond(self.transposed, mat.transpose(), mat)),

    conjTranspose := self >> ObjId(self)(self.params[1], self.params[2], -self.params[3]).transpose(),

    isReal := False,

    hashAs := self >> let(t:=ObjId(self)(self.params[1], 1/16, 1).withTags(self.getTags()),
        When(self.transposed, t.transpose(), t)),

    normalizedArithCost := self >> let(n := self.params[1],
        floor(5.0 * n * log(n) / log(2.0))),
));

Class(OnlineDiag, Diag, rec(
    isReal := self >> let(t := self.element.range(),
	Cond(IsVecT(t), t.t <> TComplex and ObjId(t.t) <> T_Complex,
	     t <> TComplex and ObjId(t)<>T_Complex)),
));

Global.compiler.DefaultCodegen.OnlineDiag := (self, o, y, x, opts) >> let(
    cc := o.element.compute(),
    decl(cc.values,
        chain(
            cc.code,
            List([1..Length(cc.values)], i -> assign(nth(y, i-1), cc.values[i] * nth(x, i-1)))
        )
    ));

Class(RCOnlineDiag, Diag, rec(
    dims := self >> Replicate(2, 2*self.element.domain())
));

Global.compiler.DefaultCodegen.RCOnlineDiag := (self, o, y, x, opts) >> let(
    cxmul      := (y, x, c) -> assign(y, c*x),
    cxmul_conj := (y, x, c) -> assign(y, conj(c)*x),
    cc         := o.element.computeWithOps(cxmul, cxmul_conj),

    decl(cc.values, chain(
	cc.code,
	List([0..Length(cc.values)-1], i -> chain(
		assign(nth(y, 2*i),   nth(x, 2*i) * re(cc.values[i+1]) - nth(x, 2*i+1) * im(cc.values[i+1])),
		assign(nth(y, 2*i+1), nth(x, 2*i) * im(cc.values[i+1]) + nth(x, 2*i+1) * re(cc.values[i+1]))))
    ))
);

RewriteRules(RulesDiag, rec(
   RC_OnlineDiag := Rule([RC, [OnlineDiag, @(1)]], e-> RCOnlineDiag(@(1).val))
));

#F SkewTwid(<n>, <a>)  -- twiddle factor function II(n-1) -> TComplex
#F   j -> omegapi(2*a*(j+1) / n)
#F
Class(SkewTwid, RewritableObject, Function, rec(
    updateParams := self >> Checked(IsPosIntSym(self.params[1]), IsRatSym(self.params[2]), 0),

    lambda := self >> let(
        n := self.params[1], a := self.params[2], j := Ind(n-1),
        Lambda(j, omegapi(fdiv(2*a*(j+1), n)))),

    range := self >> TComplex,

    domain := self >> self.params[1]-1,
));

#F OnlineSkewTwid(<n>, <a>, <seed_exps>, <seed_twiddle_func>)  -- same as SkewTwid but twiddles are computed
#F   online from seed values given by <seed_twiddle_func>.tolist()
#F
#F Example:
#F   OnlineSkewTwid(13,1/7,[1], fCompose(SkewTwid(13,1/7), FList(TInt, [0]))).compute();
Class(OnlineSkewTwid, SkewTwid, rec(
    # params[2] is only used by super-class's .lambda(), but not here, not in compute()
    updateParams := self >> Checked(IsPosIntSym(self.params[1]), IsRatSym(self.params[2]),
                                    IsList(self.params[3]), ForAll(self.params[3], IsIntSym),
                                    IsFunction(self.params[4])),

    range := self >> self.params[4].range(),

    # compute() returns [ code := ..., values := ... ],
    compute := self >> self.computeWithOps((res,a,b) -> assign(res, a*b),
	                                   (res,a,b) -> assign(res, a*conj(b))),

    computeWithOps := meth(self, cxmul, cxmul_conj)
        local exp, W, n, c, new_exp, new_W, i, a, b, v, val, num, perm;
        n := EvalScalar(self.params[1]);

        exp := List(self.params[3], EvalScalar);
        W := List([1..Length(exp)], i -> var.fresh_t(Concat("W",StringInt(exp[i]), "_"), self.range()));
        c := List([1..Length(exp)], i -> assign(W[i], self.params[4].at(i-1)));

        # NOTE: check if stuck
        while Length(exp) <> (n-1) do
            new_exp := Cartesian(exp, exp);
            new_W := Cartesian(W, W);
            num := Length(exp);
            for i in [1..Length(new_exp)] do
                [a, b] := new_exp[i];
                if (a+b) > 0 and (a+b) < n and not ((a+b) in exp) then
                    v := var.fresh_t(Concat("W",StringInt(a+b), "_"), self.range());
                    Add(exp, a+b);
                    Add(W, v);
                    Add(c, cxmul(v, new_W[i][1], new_W[i][2]));
                fi;
                if (a-b) > 0 and (a-b) < n and not ((a-b) in exp) then
                    v := var.fresh_t(Concat("W",StringInt(a-b), "_"), self.range());
                    Add(exp, a-b);
                    Add(W, v);
                    Add(c, cxmul_conj(v, new_W[i][1], new_W[i][2]));
                fi;
            od;
            if Length(exp) = num then Error("OnlineSkewTwid.compute is stuck, not enough twiddle seed values given"); fi;
        od;
        perm := SortingPerm(exp);
        return rec(values := Permuted(W,perm), exp := Permuted(exp,perm), code := chain(c));
    end
));

Class(ATwidOnline, AGenericTag, rec(isPushTag := true));
Class(ATwidSplit, AGenericTag, rec(isPushTag := true));

NewRulesFor(SkewDFT, rec(
    SkewDFT_Base2 := rec(
	applicable := t -> t.params[1] = 2,
	apply := (t, C, Nonterms) -> let(
            den:=Denominator(t.params[2]), num:=Numerator(t.params[2]), k:=t.params[3],
	    F(2) * Diag(1,E(2*den)^(k*num)))),

    SkewDFT_toDFT := rec(
	applicable := (self, t) >> logic_or(eq(self.a.maxSize, -1), leq(t.params[1], self.a.maxSize)),
        children := t -> [[ DFT(_unwrap(t.params[1]), _unwrap(t.params[3])).withTags(
		    Filtered(t.getTags(), t -> not ObjId(t) in [ATwidOnline, ATwidSplit])) ]],

        a := rec(
            diagMode := "normal", # "normal" | "split" | "online"
            maxSize := -1
        ),

        forTransposition := true,

        onlineSeeds := n -> Cond(n<=4, [1],
                                 n<=8, [1, 3],
                                       [1, 3, n-1]),

	apply := (self, t, C, Nonterms) >> let(
            N := t.params[1], rot := t.params[3], a := t.params[2], j := Ind(N-1),
            tw := SkewTwid(N, a*rot),
            mode := Cond(t.firstTag()=ATwidOnline(),
		         # if other tags are present, do not use "online" mode
		         # this is necessary because we can't vectorize a standalone SkewDFT in this mode
		              Cond(Length(t.getTags())>1, self.a.diagMode, "online"),
                         t.firstTag()=ATwidSplit(), "split",
                         self.a.diagMode),

	    C[1] *
            Cond(mode = "split",
                    DiagCpxSplit(diagDirsum(FList(TReal, [1.0, 0.0]), fPrecompute(RCData(tw)))),
                 mode = "normal" or N=2,
                    Diag(        diagDirsum(fConst(TReal, 1, 1.0),    fPrecompute(tw))), # XXX

                 mode = "online", let(seeds := self.onlineSeeds(N),
                     DirectSum(I(1), OnlineDiag(OnlineSkewTwid(N, a*rot, seeds, fPrecompute(fCompose(tw, FList(TInt, seeds-1))))))),
                 Error("<self>.diagMode must be \"split\" | \"normal\" | \"online\""))
        )
    ),

    #F SkewDFT_Fact : PDFT_2n_a -> L (PDFT_n_r0 dirsum PDFT_n_r1) L
    #F                               (I2 tensor [[1,r0],[1,r1]]) L
    #F Derived using polynomial factorization:
    #F    (x^2n - a) == (x^n - r0) * (x^n - r1)
    #F where r0 and r1 are two different quadratic roots of a
    #F
    SkewDFT_Fact := rec(
	applicable := t -> t.params[1] > 2 and t.params[1] mod 2 = 0,

	children := t -> let(
	    n:=t.params[1]/2, r0:=t.params[2]/2, r1:=r0+1/2, k:=t.params[3],
	    [[ SkewDFT(n, r0, k), SkewDFT(n, r1, k) ]]),

	apply := (t, C, Nonterms) -> let(
	    n:=t.params[1]/2, den:=Denominator(t.params[2]), num:=Numerator(t.params[2]), k:=t.params[3],
	    DirectSum(C[1], C[2]) ^ L(2*n, 2) *
	    Tensor(I(n), F(2)*Diag(1,E(2*den)^(k*num))) *
	    L(2*n, n))
    ),

    #F SkewDFT_CT : SkewDFT_mn(r) -> (SkewDFT_m(.) tensor I_n) (I_m tensor SkewDFT_n(r)) L
    #F
    #F Derived using polynomial decomposition
    #F    x^mn - a == (x^m)^n - a
    #F
    SkewDFT_CT := rec(
        applicable := (self, t) >> logic_and(gt(t.params[1], 2), logic_neg(t.hasTags()), logic_neg(isPrime(t.params[1]))),
        freedoms := t -> [ divisorsIntNonTriv(t.params[1]) ],
        child := (t, fr) -> let(
            N := t.params[1],  m := fr[1],          n := div(N, m),
            a := t.params[2],  rot := t.params[3],  j := Ind(n),

            [ IDirSum(j, SkewDFT(m, fdiv(a+j, n), rot)),
              SkewDFT(n, a, rot) ]),

        apply := (nt, C, Nonterms) -> let(m := Rows(Nonterms[1].child(1)), n := Rows(Nonterms[2]),
            Tr(n, m) * C[1] * Tr(m, n) * Tensor(I(m), C[2]) * Tr(n, m)
        )
    )
));

ext: gi