ism.m

自己编写的求可达矩阵

% ISM 解释型结构模型 根据邻接矩阵求可达矩阵，进行级划分的算法，2002-9-24，王新宇
% step 1 求可达矩阵
% A 是邻接矩阵，可以由用户输入
A = [1	0	0	0	1	0	0	1	1	0	1	1	1	0	0	0	1	1	1	1	1	0	1	1	1	1
1	1	0	0	1	1	0	1	1	0	1	1	1	0	1	0	1	1	1	1	1	0	1	1	1	1
0	0	1	0	0	0	0	1	1	0	0	1	1	1	0	0	1	1	1	1	1	0	1	1	1	1
0	0	0	1	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	0	0	1	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	1	0	0	1	1	0	1	0	0	0	0	0	0	0	1	1	1	0	0	1	0	0
0	0	0	0	1	0	0	1	1	0	1	0	0	0	0	0	0	0	1	1	1	0	0	1	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	1	1	1	1	0	0	0	1	1	1
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	0	0	0	1	1	1
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1
];


N=size(A,1);  % N 是矩阵的阶数，是所有要素的数目
r = [N:N];
r=A;
for(n=1:1:N)
  for(i=1:1:N)
    for(j=1:1:N)
     
       sum = 0.0;
           for(k=1:1:N)
             sum=sum+r(i,k)*A(k,j);%此处是采用普通的矩阵乘法，但是下面可以根据sum的值是否大于或等于1
                                   %来判断连通型，从而等价得到r(i,j)的实际值。
           end   
            if(sum >= 1)
                r(i,j)= 1;
              else
               r(i,j)=0;
            end
    end
   end
  end
  R=r; % R 是可达矩阵
  
  %step 2  级别划分
  
  L=[N:N];  %  L 是二维数组，存储级划分的结果，下面首先初始化为0
      for(i=1:1:N)
            for(j=1:1:N)
              L(i,j)=0;
          end
        end
        
for(p=1:1:N) % p 是存储层次级数 l 的变量，也控制了总循环的次数
  k=1;       % k 是记录每一级(p)内要素，在L内存储下标的变量,形式为L(p,k)
  
  %下面是利用二重循环，求解p级内的要素
  for(i=1:1:N)
      sign=0;  % 是标志变量，初始为0，如果对某要素考察后，其值仍为0，则表明该要素是顶点要素。否则，不是顶点要素
      sum=0; % sum 是计数器，用于判断当前考察的要素所在的矩阵行向量是否是值全为0的向量。因为本程序的算法是这样的，
      % 即如果已经发现某要素是某级的顶点，则在求下一级顶点时，需要去除上面各级的要素。在本程序是通过变通的方法实现
      % 相同的目的，即所有已经是顶点的要素的所在行和列的值全部重新置为0。这样就等价于删除了这些上级顶点。从而简化了
      % 程序的算法。
    for(j=1:1:N)
          sum=sum+r(i,j); 
         if(r(i,j) == 1 & (r(i,j) ~= r(j,i)))  
             % 算法的关键，R(i)=R(i)交A(i)，只需要判断要素i所在的行中，
              % 所有值为1的矩阵元素R(i,j), 其对称矩阵元素R(j,i)的值如果也为1，则说明R(i)=R(i)交A(i)成立。否则
              % 只要有一个元素不满足该条件（此时标志变量的值赋为 1 ），就说明该要素i不是顶点。进行下一要素i+1的考察。
             sign=1
             break
         end
    end
        if(sum~=0)  %sum不为0,说明该要素不是已经求出的上级的顶点，可能是新级别的顶点，如果sum为0，则说明是已经
                    % 求出的顶点
             if(sign==0)  % sign为0，说明是新顶点
                L(p,k)=i;  % 在L内记录新顶点，p为级别
                k=k+1;
             end
        end
        
   end

   % 已经对p级别的顶点计算完毕，下面的程序是把p级的顶点所在的行和列的矩阵元素全部置为0，以在后面的计算表示老顶点
          for(g=1:1:N)
            if(L(p,g)~=0)
                 for(h=1:1:N)
                   r(L(p,g),h)=0;   
                   r(h,L(p,g))=0;
                end
           end
         end
   
   % 进行下一次循环，找出p+1级别的顶点元素
   
  end

   clc  % 清命令窗口
        % 输出邻接矩阵，可达矩阵，级别划分矩阵。
   A
   R 
   L
 % 对于同级内部的强连通集的最大回路集合的划分，和不连通集的划分，比较容易。程序不再给出算法，自行判断。