ClearAll["Global`*"]
(*specify size of the considered graphs*)

DeleteIso[gs_List]:=DeleteDuplicatesBy[gs,CanonicalGraph]

(*p number of variables; sn simple is 0, non-simple is 1; edgn specifies the edgenumber max is Binomial[p,2]; el exact edgenumber is 0, less is 1 *)

desiredGraphs[p0_,sn0_,edgn0_,el0_]:=
Module[{p=p0,sn=sn0,edgn=edgn0,el=el0,L,i,j,graphlist,ggraphlist,index},
(*First get all potential graphs as edgelist for edgn many edges exact or less depending on el*)
L=ConstantArray[0,{p,p}];
For[i=1,i≤p,i++,
For[j=1,j≤p,j++,
If[i ≠ j,
L[[i,j]]={i,j}
]
]
];
L=DeleteCases[Flatten[L,1],0];
If[el==0,
graphlist= Subsets[L,{edgn}];
,If[el==1,graphlist= Subsets[L,edgn];
graphlist=graphlist[2;;];
]
];
(*convert the edge lists to actual graph objects in Mathematica*)
ggraphlist=Graph[DirectedEdge@@@ #,VertexLabels->"Name"]&/@graphlist;
(*get indices of simple or nonsimple according to sn*)
index=List[];
For[i=1,i≤Length[ggraphlist],i++,
If[Length[Flatten[FindCycle[ggraphlist[[i]],{2}]]]>0,
index=AppendTo[index,i]];
];
If[sn==0,
ggraphlist=ggraphlist[[Complement[Table[x,{x,1,Length[graphlist]}],index]]];
,If[sn==1,ggraphlist=ggraphlist[[index]];]
];
(*Remove graphs that are the same up to relabelling*)
ggraphlist=DeleteIso[ggraphlist];
Return[ggraphlist];
]

p=5;
sn=1;
edgn=Binomial[p,2];
el=0;

(*apply function and directly display the graphs, be careful for larger numbers it might take forever to display all of them*)

graphlistnice=desiredGraphs[p,sn,edgn,el];

(*Converts the graphobject into edgelists for Macauly2 or for the subsequent genericidentifiability program*)

getgraphlist[graphlist0_]:=
Module[{graphlist=graphlist0,i,resgraphlist},
resgraphlist=List[];
For[i=1,i≤Length[graphlist],i++,
resgraphlist=AppendTo[resgraphlist,List@@@EdgeList[graphlist[[i]]]];
];
Return[resgraphlist]
]

graphlistlist=getgraphlist[graphlistnice]

(*If you have a reason to only consider connected graphs in the sense that the skeleton is connected*)

graphconv[graphlist0_,p0_]:=
Module[{graphlist=graphlist0,gc,p=p0},
If[ContainsAll[Flatten[graphlist],Table[x,{x,1,p}]], gc=ConnectedGraphQ[Graph[UndirectedEdge@@@ graphlist,VertexLabels->"Name"]],
gc=False;
];
Return[gc];
]

conind=graphconv[#,p]&/@graphlistlist;
graphlistlist=Pick[graphlistlist,conind];

(*Function for checking generic Identifiablity for a specific graph*)
genericidentifiability[graphlist0_,p0_]:=
Module[{graphlist=graphlist0,p=p0,M,M1,M2,S1,S2,L1,L2,A1,A2,L,M1t,M2t,graphc,index,zeroentries,nonzeroentries,nonzeroentriesl,i,j,k,A1l,M1l},


index=List[];
A1l=List[];
M1l=List[];
nonzeroentriesl=List[];
M1=ArrayReshape[RandomPrime[100,p^2],{p,p}];
M2=ArrayReshape[RandomPrime[100,p^2],{p,p}];

L=ConstantArray[0,{p,p}];
For[i=1,i≤p,i++,
For[j=1,j≤p,j++,
If[i ≠ j,
L[[i,j]]={i,j}
]
]
];
L=DeleteCases[Flatten[L,1],0];

For[i=1,i≤Length[graphlist],i++,
M1t=M1;
M2t=M2;
graphc=Complement[L,graphlist[[i]]];
For[j=1,j≤Length[graphc],j++,
M1t [[graphc[[j]][[1]],graphc[[j]][[2]]]]=0;
M2t [[graphc[[j]][[1]],graphc[[j]][[2]]]]=0;
];
zeroentries=List[];
For[k=1,k≤p,k++,
For[j=1,j≤p,j++,
If[M1t[[j,k]]==0,zeroentries=AppendTo[zeroentries,(k-1)*p+j]];
];
];
nonzeroentries=Complement[Table[x,{x,1,p^2}],zeroentries];
(*Print[nonzeroentries];*)
S1=LyapunovSolve[M1t,-2*IdentityMatrix[p]];
S2=LyapunovSolve[M2t,-2*IdentityMatrix[p]];
M=Array[Subscript[m,##]&,{p,p}];
L1=M.S1+S1.Transpose[M]+IdentityMatrix[p];
L2=M.S2+S2.Transpose[M]+IdentityMatrix[p];
A1=D[Flatten[L1],{DeleteCases[Flatten[Transpose[M]],0],1}];
A1=DeleteDuplicates[A1][[All,nonzeroentries]];

A1l=AppendTo[A1l,A1];
M1l=AppendTo[M1l,M1t];
nonzeroentriesl=AppendTo[nonzeroentriesl,nonzeroentries];

A2=D[Flatten[L2],{DeleteCases[Flatten[Transpose[M]],0],1}];
A2=DeleteDuplicates[A2][[All,nonzeroentries]];
If[MatrixRank[A1]==Dimensions[A1][[2]]||MatrixRank[A2]==Dimensions[A2][[2]],
index=AppendTo[index,i];
];
];
Return[List[M1l,A1l,index,nonzeroentriesl]];
]

result=genericidentifiability[graphlistlist,p];

Complement[Table[x,{x,1,Length[graphlistlist]}],result[[3]]]

getinfo[result0_,index0_]:=
Module[{result=result0,index=index0},
Print[MatrixForm[result[[1]][[index]]]];
Print[MatrixForm[result[[2]][[index]]]];
Print[result[[4]][[index]]];
]

getinfo[result,3]