# # Compute the "strongly connected" blocks of a square matrix # Input an n by n square matrix M. # Output a list of 0 <= r <= n non-zero square matrices [A1,A2,...,Ar]. # Let m = sum( dimension(Ai), i=1..r ). # The value of m may be less than n because we don't return # zero matrices. # The output blocks satisfy # (i) det(M) = 0 if m < n and product( det(Ai), i=1..r ) otherwise. # (ii) charpoly(M,x) = x^(n-m) * product( charpoly(Ai,x), i=1..r ) # # Example: # > A := Matrix( [[a,b,c,d,e], # > [f,g,h,i,j], # > [0,0,0,0,k], # > [0,0,0,v,w], # > [0,0,0,x,y]] ): # > SCCBlocks(A); # # [y x] [a b] # [[ ], [ ]] # [w v] [f g] # # Author: Simon Lo, June 2006. # SCCBlocks := proc(M::Matrix) local StrongConnect,i,j,Stack,top,Number,Lowlink,OnStack,u,n,m,k,SCC,t,A; StrongConnect := proc(v) local w; i := i+1; Number[v] := i; Lowlink[v] := i; top := top+1; Stack[top] := v; OnStack[v] := true; for w in A[v] do if Number[w] = 0 then StrongConnect(w); Lowlink[v] := min(Lowlink[v],Lowlink[w]) elif Number[w] < Number[v] then if OnStack[w] then Lowlink[v] := min(Lowlink[v],Number[w]) end if end if end do; if Lowlink[v] = Number[v] then t := top; while t>0 and Number[Stack[t]] >= Number[v] do OnStack[Stack[t]] := false; t := t-1; end do; SCC[k] := [seq(Stack[j],j=t+1..top)]; top := t; k := k+1 end if end proc: n,m := op(1,M); A:=[seq([seq(`if`(M[i,j]=0,NULL,j),j=1..n)],i=1..n)]; i := 0; k := 1; Stack := Array(1..n); top := 0; Number := Array(1..n); Lowlink := Array(1..n); OnStack := Array(1..n,fill=false); for u to n do if Number[u] = 0 then StrongConnect(u) end if end do; [seq(`if`(nops(SCC[i])=1 and M[op(SCC[i]),op(SCC[i])]=0, NULL, M[SCC[i],SCC[i]]),i=1..k-1)] end proc: