Sweeping algorithm

M.P. Mello
e-mail: margarid .at. ime .dot. unicamp .dot. br

Reference: 
Conley's spectral sequence via the sweeping algorithm
M.P. Mello, K. A. de Rezende, M. R. da Silveira

Needs["Splines`"]

Auxilliary routines

Quality control subroutine

(* global variables: ul, lr, colhead (headers of delta matrix), indMorse (Morse index of column), colMorse[[j]] (first column with j-th Morse index)
Optional parameter: z2 = False --- matrix in Z, z2 = True --- matrix in Z2 *)
QualityControl=.;
QualityControl[matrix_List,columnpartition_List,OptionsPattern[z2->False]]:=
Module[{d=Dimensions[matrix][[1]],dd=Dimensions[matrix][[2]],nula,Parar,matrizbin,n,nb,lrtent,i,k,f,errormsg,maxcolumns,maxrows,mincolumns,minrows},
Clear[ul,lr,colhead,indMorse,colMorse];
If[!OptionValue[z2],f[x_]:=x;errormsg="Matrix is not nilpotent.",f[x_]:=Mod[x,2];errormsg="Matrix is not nilpotent mod 2."];
nula=Table[0,{d},{d}];
Parar=False;
(* Calculating the number and the position of the blocks *)
nb=0;k=1;
Do[If[columnpartition[[i]]*columnpartition[[i+1]]≠0,nb+=1;
ul[k]={Sum[columnpartition[[dd]],{dd,1,i-1}]+1,Sum[columnpartition[[dd]],{dd,1,i}]+1};lr[k]=ul[k]+{columnpartition[[i]]-1,columnpartition[[i+1]]-1};k+=1],{i,1,Dimensions[columnpartition][[1]]-1}];
(* Verifying whether the matrix is square, upper triangular, nilpotent, column partition is consistent, that is, the sum of h's in each group is equal to the number of columns of the matrix, first and last entries are positive. *)
If[d≠dd,Parar=True;Print["Error. Stop, matrix is not square."],If[Sum[Total[Abs[matrix[[i,Range[1,i]]]]],{i,1,d}]≠0,Parar=True;Print["Error. Stop, matrix is not upper triangular."],If[!ArrayQ[matrix,_,IntegerQ],Parar=True;Print["Error. Stop, matrix is not integral."],
If[f[matrix.matrix]≠nula,Parar=True;Print[errormsg],If[matrix==nula,Parar=True;Print["Matrix is null. Any decomposition will do.",
If[(Norm[columnpartition,1]≠d)||(columnpartition[[1]]≤0)||(columnpartition[[-1]]≤0),Parar=True;Print["Error. Column partition is inconsistent."]]]]]]]];
(* Creating headers for each column of the initial matrix and building the vector "indMorse", whose i-th entry is the group to which the i-th column belongs. *)
colhead={};indMorse={};k=0;colMorse={1};
Do[If[columnpartition[[i]]≠0,colMorse=Append[colMorse,colMorse[[-1]]+columnpartition[[i]]];indMorse=Flatten[Append[indMorse,Table[i-1,{columnpartition[[i]]}]]];colhead=Flatten[Append[colhead,Array[DisplayForm[SubsuperscriptBox["h",ToString[i-1],ToString[#]]]&,columnpartition[[i]],k+1]]]];k+=columnpartition[[i]],{i,1,Dimensions[columnpartition][[1]]}];
colMorse=Delete[colMorse,-1];
If[nb>0&&Norm[Flatten[matrix],1]==Sum[Norm[Flatten[matrix[[Range[ul[i][[1]],lr[i][[1]]],Range[ul[i][[2]],lr[i][[2]]]]]],1],{i,1,nb}],Print["Column partition is correct."];
,If[nb>0,Print["Column partition is incorrect."]]];
nb]

Change-of-basis matrix update using Mathematica's built-in minimization routine

(* global variables: ul, lr *)
PivotMudancaBase=.;
PivotMudancaBase[matrix_List,changeofbasismatrix_List,primarypivots_List,blocknumber_Integer,{k_,j_}]:=
Module[{d=Dimensions[matrix][[1]],val,sol,firstcolumn,lastcolumn=j,firstrow=k,lastrow,change=changeofbasismatrix},
firstcolumn=ul[blocknumber][[2]];
lastrow=Max[primarypivots[[ul[blocknumber][[2]];;lr[blocknumber][[2]]]]];
{val,sol}=Minimize[{Subscript[x,#]&[lastcolumn],matrix[[Range[firstrow,lastrow],Range[firstcolumn,lastcolumn]]].Array[Subscript[x,#]&,lastcolumn-firstcolumn+1,firstcolumn]==Array[0&,lastrow-firstrow+1],Subscript[x,#]&[lastcolumn]≥1},Array[Subscript[x,#]&,lastcolumn-firstcolumn+1,firstcolumn],Integers];
change[[Range[firstcolumn,lastcolumn],lastcolumn]] = Array[Subscript[x,#]&,lastcolumn-firstcolumn+1,firstcolumn]/.sol;
change
]

Change-of-basis matrix update using reduction to normal form

NormalForm=.;
NormalForm[matrix_List]:=
Module[{numlin=Dimensions[matrix][[1]],numcol=Dimensions[matrix][[2]],escalonada,i,j,ordemcol,indices,guia,complement,colunas,ultimacol,indiceaserretirado,ordemescalonada,alfa,matrizu,vetordeuns},
escalonada=matrix;
matrizu=IdentityMatrix[numcol];
indices=Table[i,{i,1,numcol}];
complement=ordemescalonada={};
vetordeuns=Table[1,{numcol}];
For[i=1;
colunas=Dimensions[indices][[1]];(* colunas sendo escalonadas *),
i≤numlin && colunas≥1,i++,
While[Count[escalonada[[i,indices]],0]<colunas - 1,
guia=escalonada[[i]];
guia[[complement]]=0;
ordemcol=Ordering[-Abs[guia]];
If[Count[guia,0]==0,ultimacol=Dimensions[guia][[1]],
ultimacol=Flatten[Position[guia[[ordemcol]],0,Infinity,1]][[1]]-1];
Do[alfa=escalonada[[i,ordemcol[[j]]]]/escalonada[[i,ordemcol[[j+1]]]];
escalonada[[Range[1,numlin],ordemcol[[j]]]]-=Sign[alfa]Floor[Abs[alfa]] escalonada[[Range[1,numlin],ordemcol[[j+1]]]];
matrizu[[Range[1,numcol],ordemcol[[j]]]]-=Sign[alfa]Floor[Abs[alfa]] matrizu[[Range[1,numcol],ordemcol[[j+1]]]];,{j,1,ultimacol-1}];
];  (* FIM do WHILE *)
guia=escalonada[[i]];
guia[[complement]]=0;
If[Abs[guia].vetordeuns==0,,
indiceaserretirado=Flatten[Position[Abs[guia],n_/;n≥1]][[1]];
ordemescalonada=Append[ordemescalonada,{i,indiceaserretirado}];
complement=Append[complement,indiceaserretirado];
indices=DeleteCases[indices,indiceaserretirado];
];
 (* FIM do IF *)
colunas=Dimensions[indices][[1]];
 (* FIM do FOR *)
];
{escalonada,ordemescalonada,matrizu}
]

colunamudabase[matrix_List,rhs_List]:=
Module[{numlin=Dimensions[matrix][[1]],numcol=Dimensions[matrix][[2]],ladodir=rhs,escalonada,ordemescalonada,matrizu,sol,alfa,suporte,j},
{escalonada,ordemescalonada,matrizu}=NormalForm[matrix];
sol=Table[0,{numcol}];
suporte=Dimensions[ordemescalonada][[1]];
Do[sol[[ordemescalonada[[j,2]]]]=ladodir[[ordemescalonada[[j,1]]]]/escalonada[[ordemescalonada[[j,1]],ordemescalonada[[j,2]]]];
ladodir-=sol[[ordemescalonada[[j,2]]]]*escalonada[[Range[1,numlin],ordemescalonada[[j,2]]]],{j,1,suporte}];
sol=matrizu.sol;
alfa=Apply[LCM,Denominator[sol]];
sol *=alfa;
sol=Append[-sol,alfa];
sol
]

(* global variables: ul, lr *)
PivotMudancaBaseFormaNormal=.;
PivotMudancaBaseFormaNormal[matrix_List,changeofbasismatrix_List,primarypivots_List,blocknumber_Integer,{k_,j_}]:=
Module[{d=Dimensions[matrix][[1]],val,sol,firstcolumn,lastcolumn=j,firstrow=k,lastrow,change=changeofbasismatrix},
firstcolumn=ul[blocknumber][[2]];
lastrow=Max[primarypivots[[ul[blocknumber][[2]];;lr[blocknumber][[2]]]]];
change[[Range[firstcolumn,lastcolumn],lastcolumn]] = colunamudabase[matrix[[Range[firstrow,lastrow],Range[firstcolumn,lastcolumn-1]]],matrix[[Range[firstrow,lastrow],lastcolumn]]];
change
]

Quotient space computation

(* alfa = 0 means system is infeasible
alfa ≥ 1 means system is feasible and alfa is the minimum value *)
CalculoEspacoQuociente=.;
CalculoEspacoQuociente[matrix_List]:=
Module[{numlin=Dimensions[matrix][[1]],numcol=Dimensions[matrix][[2]],matcoef,escalonada,ordemescalonada,matrizu,ladodir,sol,suporte,alfa},
ladodir=matrix[[All,numcol]];
matcoef=matrix[[All,Range[1,numcol-1]]];
{escalonada,ordemescalonada,matrizu}=NormalForm[matcoef];
sol=Table[0,{numcol-1}];
suporte=Dimensions[ordemescalonada][[1]];
Do[sol[[ordemescalonada[[j,2]]]]=ladodir[[ordemescalonada[[j,1]]]]/escalonada[[ordemescalonada[[j,1]],ordemescalonada[[j,2]]]];
ladodir-=sol[[ordemescalonada[[j,2]]]]*escalonada[[Range[1,numlin],ordemescalonada[[j,2]]]],{j,1,suporte}];
If[ladodir.Table[1,{numlin}]==0,alfa=Apply[LCM,Denominator[sol]],alfa=0];
alfa]

Sweeping algorithm over ℤ with quotient space computation

(* global variables: ul, lr, pp, pmd, changeofbasismatrix, colhead, indMorse, colMorse *)
(* Optional parameter Fn: If true, use alternate routine (that employs the integer normal form) to update the change-of-basis matrix
If false, use built-in function Minimize *)
SpectralSequenceConstruction=.;
SpectralSequenceConstruction[matrix_List,columnpartition_List,OptionsPattern[{Nf->False,Erp->False}]]:=
Module[{d=Dimensions[matrix][[1]],Parar,firstcolumn,firstrow,lastrow,lastcolumn,m,n,i,j,k,currentmatrix,iter,marcasprimarias,marcasmudancabase,marcas,pivotsprimarios,nb,col,row,mu,musigma,erp,drp,bordo,matnumerador,matdenominador1,matdenominador2,musigmahead,aux,nMorse,poscolumnpart},
(* Initialization *)
Clear[ul,lr,pp,pmd,colhead,changeofbasismatrix];
Parar=False;
currentmatrix=matrix;
bordo[0]=matrix;
(* Verifying input *)
(* nb = number of blocks *)
nb=QualityControl[currentmatrix,columnpartition];
poscolumnpart=Delete[columnpartition,Position[columnpartition,0]];
nMorse=Dimensions[colMorse][[1]];
(* iter = number of iterations = number of diagonals to be swept *)
If[nb==0,Parar=True,iter=Max[Array[lr[#][[2]]-ul[#][[1]]&,nb]]];
If[Parar,Return["error"]];
marcas={};
Print["Total number of diagonals to be swept = ",iter];
pmd=Table[{},{iter}];
pivotsprimarios=Table[{},{iter}];
(* pp[[j]] = 0, se coluna j ainda nao contem pivot primário
	= i, se col. j contem pivot primario na linha i *)
marcasprimarias={LightYellow,EdgeForm[Thick]};
pp=Table[0,{d}];
Print["i = 0"];
Print[Text[Style[Superscript[Δ,#]&[0]]]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[{Text[colhead[[d+1-col]],{-1,col-.2},{1,1}],Text[colhead[[col]],{col-.5,d+1},{0,0}]},{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->{600,500}]];
Print["---------",0,"---------",0,"---------",0,"---------"];
changeofbasismatrix[0]=IdentityMatrix[d];
(* escalonada[i] contém o vetor (μ^(.,i+1)) *)
mu[0]=Table[1,{d}];
musigma[i_]=changeofbasismatrix[i].DiagonalMatrix[mu[i]];
erp[0]=Table["",{d}];
Do[(* i-th diagonal is swept *)
changeofbasismatrix[i]=changeofbasismatrix[i-1];
marcasmudancabase={LightBlue,EdgeForm[Directive[Thick,Dotted]]};
mu[i]=mu[i-1];
Do[(* entries in block n are examined *)
{firstrow,firstcolumn}=ul[n];
{lastrow,lastcolumn}=lr[n];
Do[(*  entry in column j is examined *)
(* Varrendo elemento (k,j) na diagonal auxiliar i, portanto k = j-i *)
k=j-i;
(* Mark-up step *)
If[currentmatrix[[k,j]]≠0,
If[pp[[j]]>k,(* Do not mark, primary pivot below *),If[MemberQ[pp[[1;;j]],k],(* Primary pivot on the left. Mark entry (k,j) as change-of-basis pivot *)
pmd[[i]]=Append[pmd[[i]],{k,j}];marcasmudancabase=Append[marcasmudancabase,Rectangle[{j-1,d-k},{j,d-k+1}]];
(* Change-of-basis matrix update *)
If[OptionValue[Nf],changeofbasismatrix[i]=PivotMudancaBaseFormaNormal[matrix,changeofbasismatrix[i],pp,n,{k,j}];,changeofbasismatrix[i]=PivotMudancaBase[matrix,changeofbasismatrix[i],pp,n,{k,j}];]
,
aux=mu[i];
aux[[j]]=0;
mu[i]=aux;
pp[[j]]=k;pivotsprimarios[[i]]=Append[pivotsprimarios[[i]],{k,j}];marcas=Append[marcas,{k,j}];
marcasprimarias=Append[marcasprimarias,Rectangle[{j-1,d-k},{j,d-k+1}]];]]],
{j,Max[{firstcolumn,i+firstrow}],Min[{lastcolumn,i+lastrow}]}];,(* END of DO in j (column index) *)
{n,1,nb}];
(* END of DO in n (block index) *)
Print["i = ",i];
Print[Superscript[Δ,#]&[i]," = ",Show[Graphics[{{GrayLevel[.7],Thickness[0.01],Line[{{i,d},{d,i}}]},{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,indMorse[[d+1-row]],ToString[d+1-row]<>","<>ToString[i]],"=",(colhead.changeofbasismatrix[i-1])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,indMorse[[col]],ToString[col]<>","<>ToString[i]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],marcasprimarias,marcasmudancabase,Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->{600,500}]];
bordo[i]=matrix.changeofbasismatrix[i];
(* Computation of (Subscript[E, j]^i) and (Subscript[d, j]^i) --- BEGIN *)
If[OptionValue[Erp],
Clear[aux];aux=Table["{0}",{d}];drp=Table[0,{d}];
Do[(* Do in n from 1 to nMorse, the number of distinct Morse indices *)
Do[(* Do in j, column with n-th Morse index *)
Print[SubsuperscriptBox["E",j-1,i]//DisplayForm];
If[erp[i-1][[j]]==="{0}",,
Clear[matnumerador,matdenominador1,matdenominador2,musigmahead];
matnumerador=Transpose[Table[musigma[Max[i-1-(j-k),0]][[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1],k]],{k,colMorse[[n]],j}]];
musigmahead=Table[DisplayForm[RowBox[{SuperscriptBox[μ,ToString[k]<>","<>ToString[i-(j-k)]],SubsuperscriptBox[σ,indMorse[[k]],ToString[k]<>","<>ToString[i-(j-k)]]}]],{k,colMorse[[n]],j}];
Print["Numerator:  ",
TableForm[matnumerador,TableHeadings->{colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],musigmahead}]];
If[Abs[matnumerador[[All,j-colMorse[[n]]+1]]].Table[1,{poscolumnpart[[n]]}]==0,,
 Print["Denominator:  "];
If[j==colMorse[[n]], Print["There is no first term denominator"];,matdenominador1=Transpose[Table[musigma[Max[i-1-(j-k),0]][[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1],k]],{k,colMorse[[n]],j-1}]];
musigmahead=Table[DisplayForm[RowBox[{SuperscriptBox[μ,ToString[k]<>","<>ToString[i-(j-k)]],SubsuperscriptBox[σ,indMorse[[k]],ToString[k]<>","<>ToString[i-(j-k)]]}]],{k,colMorse[[n]],j-1}];
 Print["1st term of denominator: ",TableForm[matdenominador1,TableHeadings->{colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],musigmahead}]];
];
If[n==nMorse,Print["There is no second term denominator"];aux[[j]]=DisplayForm[RowBox[{"ℤ[",matnumerador[[All,j-colMorse[[n]]+1]].colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],"]"}]], 
If[(j-1+i<colMorse[[n]]+poscolumnpart[[n]])||(indMorse[[colMorse[[n]]]]+1≠indMorse[[colMorse[[n+1]]]]) , Print["There is no second term denominator"];aux[[j]]=DisplayForm[RowBox[{"ℤ[",matnumerador[[All,j-colMorse[[n]]+1]].colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],"]"}]],matdenominador2=Transpose[Table[(bordo[Max[k-j-1,0]].DiagonalMatrix[mu[Max[k-j-1,0]]])[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1],k]],{k,colMorse[[n+1]],Min[j-1+i,colMorse[[n+1]]+poscolumnpart[[n+1]]-1]}]];
musigmahead=Table[DisplayForm[RowBox[{SuperscriptBox["∂μ",ToString[k]<>","<>ToString[k-j]],SubsuperscriptBox[σ,indMorse[[k]],ToString[k]<>","<>ToString[k-j]]}]],{k,colMorse[[n+1]],Min[j-1+i,colMorse[[n+1]]+poscolumnpart[[n+1]]-1]}];
 Print["2nd term of denominador:  ",TableForm[matdenominador2,TableHeadings->{colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],musigmahead}]];
m=Last[Flatten[Position[Abs[matnumerador[[All,j-colMorse[[n]]+1]]],m_/;m≥1]]];
If[Abs[matdenominador2[[m,All]]].Table[1,{Min[j-1+i,colMorse[[n+1]]+poscolumnpart[[n+1]]-1]-colMorse[[n+1]]+1}]==0,aux[[j]]=DisplayForm[RowBox[{"ℤ[",matnumerador[[All,j-colMorse[[n]]+1]].colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],"]"}]];,
matdenominador2=Transpose[Partition[Flatten[Append[Flatten[Transpose[matdenominador2]],Flatten[Transpose[matnumerador]]]],poscolumnpart[[n]]]];
m=CalculoEspacoQuociente[matdenominador2];
If[m≤1,,aux[[j]]=DisplayForm[RowBox[{SubscriptBox["ℤ",ToString[m]],"[",matnumerador[[All,j-colMorse[[n]]+1]].colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]],"]"}]]];]]]]];
,{j,colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1}];
,{n,1,nMorse}];
erp[i]=aux;
Do[If[(aux[[j]]==="{0}") || (aux[[j-i]]==="{0}"),,drp[[j]]=currentmatrix[[j-i,j]]],{j,i+1,d}];
Print[Table[DisplayForm[RowBox[{DisplayForm[SubsuperscriptBox["E",j-1,i]],"=",erp[i][[j]]}]],{j,1,d}]]; Print[Table[DisplayForm[RowBox[{DisplayForm[SubsuperscriptBox["d",j-1,i]],"=",drp[[j]]}]],{j,1,d}]];];
(* Computation of (Subscript[E, j]^i) e (Subscript[d, j]^i) --- END *)
(* Matrix Delta update *)
currentmatrix=Inverse[changeofbasismatrix[i]].matrix.changeofbasismatrix[i];
Print["---------",i,"---------",i,"---------",i,"---------"];,
{i,1,iter}];
(* END of DO in i (diagonal index) *)
Print[Superscript[Δ,#]&[iter+1]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,ToString[indMorse[[d+1-row]]],ToString[d+1-row]<>","<>ToString[iter+1]],"=",(colhead.changeofbasismatrix[iter])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,ToString[indMorse[[col]]],ToString[col]<>","<>ToString[iter+1]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->{600,500}]];
Clear[ul,lr,pp,pmd,colhead];
iter
]

Sweeping algorithm over (Subscript[ℤ, 2]) without quotient space computation

(* global variables: ul, lr, pp, pmd, colhead, changeofbasismatrix, indMorse, colMorse *)
BinarySpectralSequenceConstruction=.;
BinarySpectralSequenceConstruction[matrix_List,columnpartition_List]:=
Module[{d=Dimensions[matrix][[1]],Parar,firstcolumn,firstrow,lastrow,lastcolumn,n,i,j,k,currentmatrix,iter,marcasprimarias,marcasmudancabase,marcas,ind,nula,pivotsprimarios},
Clear[ul,lr,pp,pmd,colhead,changeofbasismatrix];
nula=Table[0,{d},{d}];
Parar=False;
currentmatrix=matrix;
nb=QualityControl[currentmatrix,columnpartition,z2->True];
(* iter = numero de iteracoes = numero de diagonais auxiliares a serem varridas *)
If[nb==0,Parar=True,iter=Max[Array[lr[#][[2]]-ul[#][[1]]&,nb]]];
If[Parar,Return["erro"]];
marcas={};
Print["iter = ",iter];
pmd=Table[{},{iter}];
pivotsprimarios=Table[{},{iter}];
(* pp[[j]] = 0, se coluna j ainda não contém pivot primário
	= i, se col. j contém pivot primário na linha i *)
marcasprimarias={LightYellow,EdgeForm[Thick]};
pp=Table[0,{d}];
Print["i = 0"];
Print[Superscript[Δ,#]&[0]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[{Text[colhead[[d+1-col]]//DisplayForm,{-1,col-.2},{1,1}],Text[colhead[[col]]//DisplayForm,{col-.5,d+1},{0,0}]},{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
Print["---------",0,"---------",0,"---------",0,"---------"];
changeofbasismatrix[0]=IdentityMatrix[d];
Do[(* em i *)
changeofbasismatrix[i]=changeofbasismatrix[i-1];
marcasmudancabase={LightBlue,EdgeForm[Directive[Thick,Dotted]]};
Do[(* em n *)
{firstrow,firstcolumn}=ul[n];
{lastrow,lastcolumn}=lr[n];
Do[(* em j *)
(* Varrendo elemento (k,j) na diagonal auxiliar i, portanto k = j-i *)
k=j-i;
If[currentmatrix[[k,j]]≠0,If[pp[[j]]>k,(*Print["pp[[j]]>k, nao marca, existe pivot prim abaixo"]*),If[MemberQ[pp[[1;;j]],k],(*Print["Existe pivot primario a esquerda, faz mudanca de base em {",k,",",j,"}"];*)
ind=Flatten[Position[pp[[1;;j]],k]][[1]];
pmd[[i]]=Append[pmd[[i]],{k,j}];
marcasmudancabase=Append[marcasmudancabase,Rectangle[{j-1,d-k},{j,d-k+1}]];
changeofbasismatrix[i]=changeofbasismatrix[i].ReplacePart[IdentityMatrix[d],1,{ind,j}];
(* Print["changeofbasismatrix[",i,"] = ",mudabase[i]//MatrixForm];*),
pp[[j]]=k;
pivotsprimarios[[i]]=Append[pivotsprimarios[[i]],{k,j}];marcas=Append[marcas,{k,j}];
marcasprimarias=Append[marcasprimarias,Rectangle[{j-1,d-k},{j,d-k+1}]];(*Print["Marca como pivot primario, pp[[",j,"]] = ",k]*)]]],
{j,Max[{firstcolumn,i+firstrow}],Min[{lastcolumn,i+lastrow}]}];,
{n,1,nb}];
Print["i = ",i];
Print[Superscript[Δ,#]&[i]," = ",Show[Graphics[{{GrayLevel[.7],Thickness[0.01],Line[{{i,d},{d,i}}]},{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,marcasmudancabase,
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,ToString[indMorse[[d+1-row]]],ToString[d+1-row]<>","<>ToString[i]],"=",(colhead.changeofbasismatrix[i])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,ToString[indMorse[[col]]],ToString[col]<>","<>ToString[i]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
currentmatrix=Mod[Inverse[changeofbasismatrix[i]].matrix.changeofbasismatrix[i],2];
Print["---------",i,"---------",i,"---------",i,"---------"];,
{i,1,iter}];
Print[Superscript[Δ,#]&[iter+1]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,ToString[indMorse[[d+1-row]]],ToString[d+1-row]<>","<>ToString[iter+1]],"=",(colhead.changeofbasismatrix[iter])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,ToString[indMorse[[col]]],ToString[col]<>","<>ToString[iter+1]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
Clear[ul,lr,pp,pmd,colhead];
iter
]

Sweeping algorithm over (Subscript[ℤ, 2]) without quotient space computation, with (Subscript[H, k]) computation

leading[row_]:=If[Norm[row]>0,Position[row,1][[1,1]]];
kernelbasis[matrix_List]:=Module[{reduced,basis={},ncols=Dimensions[matrix][[2]],zeros,sup},
zeros=Table[0,{ncols}];
If[Norm[Flatten[matrix]]>0,reduced=RowReduce[matrix,Modulus->2];
leaders=DeleteCases[Map[leading,reduced],Null];
indices=Complement[Range[1,ncols],leaders];
Do[sup=Join[leaders[[Flatten[Position[reduced[[All,indices[[i]]]],1]]]],{indices[[i]]}];
basis=Append[basis,ReplacePart[zeros,Thread[sup->Table[1,{Length[sup]}]]]],{i,1,Length[indices]}]];
basis
]

drawedge[{x_,y_}]:=Spline[{{y,0},{.2 x+0.8 y,.2(x-y)},{.5(y+x),.3(x-y)},{.8 x+0.2 y,.2(x-y)},{x,0}},Bezier];

drawgraph[matrix_List,partition_List,n_Integer]:=Module[{A=matrix,dd=Dimensions[matrix][[1]],unsd,nn,nnb,numunsbloco,inicio,fim,cores,edges},
(* n = number of distinct Morse indices *)
(* at most n-1 sets of arcs between nodes with different indices *)
nnb=n-1;
unsd=Table[1,{dd}];
numunsbloco=Table[0,{nnb}];
Do[inicio=Sum[partition[[j]],{j,1,nn}]+1;fim=Sum[partition[[j]],{j,1,nn+1}];numunsbloco[[nn]] = (Table[1,{fim+1-inicio}].A[[Range[inicio,fim],All]]).unsd,{nn,1,nnb }];
cores={};
Do[cores=Append[cores,Array[Hue[nn/nnb]&,numunsbloco[[nn]]]],{nn,1,nnb}];
cores=Flatten[cores];
edges=Map[(-#+{dd+1,dd+1})&,Position[A,1]];
Print[Show[Graphics[{{EdgeForm[Thin],White,Table[Disk[{i,-.3},.3],{i,1,dd}],Black,Table[Text[dd+1-i,{i,-.3}],{i,1,dd}]},MapThread[{#1,#2}&,{cores,Map[drawedge,edges]}]}],ImageSize->600]]
];

(* global variables: ul, lr, pp, pmd, colhead, changeofbasismatrix, indMorse, colMorse *)
BinarySpectralSequenceConstructionWithHk=.;
BinarySpectralSequenceConstructionWithHk[matrix_List,columnpartition_List]:=
Module[{d=Dimensions[matrix][[1]],Parar,firstcolumn,firstrow,lastrow,lastcolumn,n,i,j,k,currentmatrix,iter,marcasprimarias,marcasmudancabase,marcas,ind,nula,pivotsprimarios,imageBasis,kBasis,hkBasis,ihk,poscolumnpart,nMorse},
Clear[ul,lr,pp,pmd,colhead,changeofbasismatrix];
nula=Table[0,{d},{d}];
Parar=False;
currentmatrix=matrix;
nb=QualityControl[currentmatrix,columnpartition,z2->True];
(* poscolumnpart = vector with positive components of columnpartition *)
poscolumnpart=Delete[columnpartition,Position[columnpartition,0]];
nMorse=Dimensions[colMorse][[1]];
(* iter = numero de iteracoes = numero de diagonais auxiliares a serem varridas *)
If[nb==0,Parar=True,iter=Max[Array[lr[#][[2]]-ul[#][[1]]&,nb]]];
If[Parar,Return["erro"]];
marcas={};
Print["iter = ",iter];
pmd=Table[{},{iter}];
pivotsprimarios=Table[{},{iter}];
(* pp[[j]] = 0, se coluna j ainda não contém pivot primário
	= i, se col. j contém pivot primário na linha i *)
marcasprimarias={LightYellow,EdgeForm[Thick]};
pp=Table[0,{d}];
Print["i = 0"];
Print[Superscript[Δ,#]&[0]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[{Text[colhead[[d+1-col]]//DisplayForm,{-1,col-.2},{1,1}],Text[colhead[[col]]//DisplayForm,{col-.5,d+1},{0,0}]},{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
Print["---------",0,"---------",0,"---------",0,"---------"];
changeofbasismatrix[0]=IdentityMatrix[d];
Do[(* em i *)
changeofbasismatrix[i]=changeofbasismatrix[i-1];
marcasmudancabase={LightBlue,EdgeForm[Directive[Thick,Dotted]]};
Do[(* em n *)
{firstrow,firstcolumn}=ul[n];
{lastrow,lastcolumn}=lr[n];
Do[(* em j *)
(* Varrendo elemento (k,j) na diagonal auxiliar i, portanto k = j-i *)
k=j-i;
If[currentmatrix[[k,j]]≠0,If[pp[[j]]>k,(*Print["pp[[j]]>k, nao marca, existe pivot prim abaixo"]*),If[MemberQ[pp[[1;;j]],k],(*Print["Existe pivot primario a esquerda, faz mudanca de base em {",k,",",j,"}"];*)
ind=Flatten[Position[pp[[1;;j]],k]][[1]];
pmd[[i]]=Append[pmd[[i]],{k,j}];
marcasmudancabase=Append[marcasmudancabase,Rectangle[{j-1,d-k},{j,d-k+1}]];
changeofbasismatrix[i]=changeofbasismatrix[i].ReplacePart[IdentityMatrix[d],1,{ind,j}];
(* Print["changeofbasismatrix[",i,"] = ",mudabase[i]//MatrixForm];*),
pp[[j]]=k;
pivotsprimarios[[i]]=Append[pivotsprimarios[[i]],{k,j}];marcas=Append[marcas,{k,j}];
marcasprimarias=Append[marcasprimarias,Rectangle[{j-1,d-k},{j,d-k+1}]];(*Print["Marca como pivot primario, pp[[",j,"]] = ",k]*)]]],
{j,Max[{firstcolumn,i+firstrow}],Min[{lastcolumn,i+lastrow}]}];,
{n,1,nb}];
Print["i = ",i];
Print[Superscript[Δ,#]&[i]," = ",Show[Graphics[{{GrayLevel[.7],Thickness[0.01],Line[{{i,d},{d,i}}]},{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,marcasmudancabase,
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,ToString[indMorse[[d+1-row]]],ToString[d+1-row]<>","<>ToString[i]],"=",(colhead.changeofbasismatrix[i])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,ToString[indMorse[[col]]],ToString[col]<>","<>ToString[i]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
(* Computation of (Subscript[H, k])(C)=((Ker (Subscript[Δ, k]))/(Im (Subscript[Δ, k+1]))) --- BEGIN *)
Do[(* Do in n from 1 to nMorse, the number of distinct Morse indices *)
(* Calculating the basis of the kernel of Delta_n *)
If[Norm[Flatten[currentmatrix[[All,Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]]]]==0,
(* Columns of Delta^i with nth Morse index are null *)
kBasis=IdentityMatrix[poscolumnpart[[n]]];,
(* Columns of Delta^i with nth Morse index are not null, which means the nth Morse index corresponds to a block b < n *)
indbloco=Flatten[Position[Table[ul[j][[2]],{j,1,nb}],colMorse[[n]]]][[1]];
kBasis=kernelbasis[currentmatrix[[ul[indbloco][[1]];;lr[indbloco][[1]],ul[indbloco][[2]];;lr[indbloco][[2]]]]];];
(* Calculating the basis of the image of Delta_{n+1} *)
If[n==nMorse,
(* n is the last Morse index *)
Print[SubscriptBox[H,ToString[indMorse[[colMorse[[n]]]]]]//DisplayForm," = ",Table[RowBox[{Table[SubsuperscriptBox[σ,ToString[indMorse[[colMorse[[n]]]]],ToString[col]<>","<>ToString[i]],{col,colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1}].kBasis[[ihk]]//DisplayForm,"=",colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]].kBasis[[ihk]]}]//DisplayForm,{ihk,1,Length[kBasis]}]];,
If[indMorse[[colMorse[[n]]]]+1≠indMorse[[colMorse[[n+1]]]],
(* There is no Morse index = indMorse[[colMorse[[n]]]]+1 *)
If[Length[kBasis]==0,Print[DisplayForm[SubscriptBox[H,ToString[indMorse[[colMorse[[n]]]]]]]," = {0}"],Print[SubscriptBox[H,ToString[indMorse[[colMorse[[n]]]]]]//DisplayForm," = ",Table[RowBox[{Table[SubsuperscriptBox[σ,ToString[indMorse[[colMorse[[n]]]]],ToString[col]<>","<>ToString[i]],{col,colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1}].kBasis[[ihk]]//DisplayForm,"=",colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]].kBasis[[ihk]]}]//DisplayForm,{ihk,1,Length[kBasis]}]];],
(* There IS a Morse index = indMorse[[colMorse[[n]]]]+1 *)
indbloco=Flatten[Position[Table[ul[j][[2]],{j,1,nb}],colMorse[[n+1]]]][[1]];
imageBasis=RowReduce[Transpose[currentmatrix[[ul[indbloco][[1]];;lr[indbloco][[1]],ul[indbloco][[2]];;lr[indbloco][[2]]]]],Modulus->2];
(* Begin of calculation of hkbasis =  basis of H_colMorse[[n]]*)
If[Norm[Flatten[imageBasis]]==0,hkBasis=kBasis,
While[Norm[imageBasis[[-1,All]]]==0,imageBasis=Delete[imageBasis,-1]];
hkBasis={};
Do[imageBasis=RowReduce[Append[imageBasis,kBasis[[j]]],Modulus->2];
If[Norm[imageBasis[[-1,All]]]>0,hkBasis=Append[hkBasis,kBasis[[j]]],imageBasis=Delete[imageBasis,-1]],{j,1,Length[kBasis]}]];
(* End of calculation of hkbasis *)
If[Length[hkBasis]==0,Print[DisplayForm[SubscriptBox[H,ToString[indMorse[[colMorse[[n]]]]]]]," = {0}"],Print[SubscriptBox[H,ToString[indMorse[[colMorse[[n]]]]]]//DisplayForm," = ",Table[RowBox[{Table[SubsuperscriptBox[σ,ToString[indMorse[[colMorse[[n]]]]],ToString[col]<>","<>ToString[i]],{col,Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]}].hkBasis[[ihk]]//DisplayForm,"=",colhead[[Range[colMorse[[n]],colMorse[[n]]+poscolumnpart[[n]]-1]]].hkBasis[[ihk]]}]//DisplayForm,{ihk,1,Length[hkBasis]}]]]]];
,{n,1,nMorse}];
(* Computation of (Subscript[H, k])(C)=((Ker (Subscript[Δ, k]))/(Im (Subscript[Δ, k+1]))) --- END *)
drawgraph[currentmatrix,poscolumnpart,nMorse];
currentmatrix=Mod[Inverse[changeofbasismatrix[i]].matrix.changeofbasismatrix[i],2];
Print["---------",i,"---------",i,"---------",i,"---------"];
,
{i,1,iter}];
Print[Superscript[Δ,#]&[iter+1]," = ",Show[Graphics[{{GrayLevel[.8],Table[Rectangle[{ul[i][[2]]-1,d-lr[i][[1]]},{lr[i][[2]],d+1-ul[i][[1]]}],{i,1,nb}]},Table[Line[{{0,row},{d,row}}],{row,0,d}],Table[Line[{{col,0},{col,d}}],{col,0,d}],marcasprimarias,
Black,Table[Text[RowBox[{SubsuperscriptBox[σ,ToString[indMorse[[d+1-row]]],ToString[d+1-row]<>","<>ToString[iter+1]],"=",(colhead.changeofbasismatrix[iter])[[d+1-row]]}]//DisplayForm,{-1,row-.2},{1,1}],{row,1,d}],Table[Text[SubsuperscriptBox[σ,ToString[indMorse[[col]]],ToString[col]<>","<>ToString[iter+1]]//DisplayForm,{col-.5,d+1},{0,0}],{col,1,d}],Table[Table[Text[Style[currentmatrix[[row,col]],Medium,Red,Bold],{col-.5,d+1-row-.5},{0,0}],{row,ul[k][[1]],lr[k][[1]]},{col,ul[k][[2]],lr[k][[2]]}],{k,1,nb}],Thickness[.01]}],PlotRange->All,ImageSize->600]];
Clear[ul,lr,pp,pmd,colhead];
iter
]