#include"iostream.h"
#include"stdlib.h"
#include "math.h"
#include "stdio.h"
#define L 8
//计算矩阵行列式的子函数
void computer_det(int In,int Fk[L][L],float *Fdet)
{
int i,j,i1,j1;
int t;
for(j=0;j<In;j++)
{
i=j;
while(Fk[i][j]==0)
i++;
if(i>=In)
*Fdet=0;
else
{
for(j1=0;j1<In;j1++)
{
t=Fk[i][j1];
Fk[i][j1]=Fk[j][j1];
Fk[j][j1]=t;
}
*Fdet=-1*float(Fk[0][0]);
for(i1=j+1;i1<In;i1++)
for(j1=In-1;j1>=0;j1--)
Fk[i1][j1]=Fk[i1][j1]-Fk[j][j1]*Fk[i1][j]/Fk[j][j];
}
}
for(j=1;j<In;j++)
*Fdet=*Fdet*float(Fk[j][j]);
}
//矩阵求逆的子函数
void reciprocal_matrix(int n,int A[L][L],int B[L][L])
{
int i,j,k,m;
int t;
for(i=0;i<n;i++)
A[i][i+n]=1;
for(j=0;j<n;j++)
{
for(i=j,m=j; ;i<n;i++)
if(fabs(A[i][j])>fabs(A[j][j]))
m=i; //找到主元所在行
for(k=0;k<2*n;k++)
{
t=A[j][k]; //主元所在行与对角元所在行交换
A[j][k]=A[m][k];
A[m][k]=t;
}
for(i=j+1;i<n;i++)
for(k=2*n-1;k>=j;k--)
A[i][k]=A[i][k]-A[j][k]*A[i][j]/A[j][j];//化为上三角阵
}
for(j=n-1;j>0;j--)
for(i=0;i<j;i++)
for(k=2*n-1;k>=0;k--)
A[i][k]=A[i][k]-A[j][k]*A[i][j]/A[j][j];//化为单位阵
for(i=0;i<n;i++)
for(k=2*n-1;k>=0;k--)
A[i][k]=A[i][k]/A[i][i]; //对角元化为1
for(i=0;i<n;i++)
for(j=0;j<n;j++)
B[i][j]=A[i][j+n];
}
//矩阵相乘的子函数
void multiply_matrix(int m,int n,int s,int P[L][L],int Q[L][L],int M[L][L])
{
int i,j,l;
for(i=0;i<m;i++)
for(j=0;j<s;j++)
{
for(l=0;l<n;l++)
M[i][j]=M[i][j]+P[i][l]*Q[l][j];
}
}
//矩阵转置并求反的子函数
void neg_matrix(int m,int n,int K0[L][L],int K1[L][L])
{
int i,j;
for(i=0;i<m;i++)
for(j=0;j<n;j++)
K1[j][i]=-1*K0[i][j];
}
//主函数
void main()
{
int i,j,j1,j2,JF,JT,IFM[L],ITO[L],Nbr_tr; //临时变量
float det;
int NAt[L][L]={0};
int Nnode; //节点数
int Nbr; //支路数
int Ntr; //树枝数
int Tree[L]={0}; //一组树支
int Link[L]={0}; //一组连支
int A[L][L]={0}; //关联矩阵
int At[L][L]={0};
int At1[L][L]={0};
int Al[L][L]={0};
int Bf[L][L]={0}; //基本回路矩阵
int Qf[L][L]={0}; //基本割集矩阵
FILE *fp;
//从"原始数据输入.txt"读入网络图原始数据,求得A矩阵，默认以最后一个节点作为参考节点
if((fp=fopen("原始数据输入.txt","r"))==NULL)
{
printf("Can not open the file.\n");
return;
}
fscanf(fp,"节点数Nnode=%d\n支路数Nbr=%d\n",&Nnode,&Nbr);
fscanf(fp,"网络图的各支路两端的节点（以下三列依次为支路编号，支路的起始节点号，支路的终止节点号）\n");
while(fscanf(fp,"%d\t%d\t%d\n",&i,&JF,&JT)!=EOF)
{
IFM[i]=JF;
ITO[i]=JT;
A[JF][i]=1;
A[JT][i]=-1;
}
fclose(fp);
Ntr=Nnode-1;
Nbr_tr=Nbr-Ntr;
//从A阵中找出一个树,并得到At,Al矩阵
for(i=0;i<Ntr;i++)
for(j=0;j<Ntr;j++)
At[i][j]=A[i][j];
computer_det(Ntr,At,&det);
for(i=0;i<Ntr;i++)
for(j=0;j<Ntr;j++)
At[i][j]=A[i][j];
for(i=0;i<Ntr;i++)
for(j=Ntr;j<Nbr;j++)
Al[i][j-Ntr]=A[i][j];
j=0;
j1=Ntr;
while(det==0)
{
for(i=0;i<Ntr;i++)
At[i][j]=A[i][j1];
computer_det(Ntr,At,&det);
for(i=0;i<Ntr;i++)
for(j2=0;j2<Nbr;j2++)
At[i][j2]=A[i][j2];
j++;
j1++;
}
if(j>0||j1>Ntr)
{
j--;
j1--;
for(i=0;i<j;i++)
{
Tree[i]=i;
Link[i]=i+Ntr;
}
for(i=j+1;i<Ntr;i++)
{
Tree[i]=i;
Link[i]=i+Ntr;
}
Tree[j]=j1;
Link[j1-Ntr]=j;
for(i=0;i<Ntr;i++)
{
Al[i][j1-Ntr]=At[i][j];
At[i][j]=A[i][j1];
}
}
else
{
for(i=0;i<Ntr;i++)
Tree[i]=i;
for(i=0;i<Ntr;i++)
Link[i]=i+Ntr;
}
//求Bf和Qf
for(i=0;i<Ntr;i++)
for(j=0;j<Ntr;j++)
At1[i][j]=At[i][j];
reciprocal_matrix(Ntr,At1,NAt);
multiply_matrix(Ntr,Ntr,Nbr_tr,NAt,Al,Qf);
neg_matrix(Ntr,Nbr_tr,Qf,Bf);
for(i=0;i<Ntr;i++)
Bf[i][i+Ntr]=1;
for(i=0;i<Ntr;i++)
{
for(j=0;j<Ntr;j++)
{
Qf[i][j+Ntr]=Qf[i][j];
Qf[i][j]=0;
}
Qf[i][i]=1;
}
//将结果保存到"运行结果.txt"中
if((fp=fopen("运行结果.txt","w+"))==NULL)
{
printf("Can not open the file.\n");
return;
}
fprintf(fp,"\n关联矩阵A:\n");
for(i=0;i<Ntr;i++)
{
for(j=0;j<Nbr;j++)
fprintf(fp,"%d\t",A[i][j]);
fprintf(fp,"\n");
}
fprintf(fp,"\n一组树支:\n");
for(i=0;i<Ntr;i++)
fprintf(fp,"%d,",Tree[i]);
fprintf(fp,"\n相应的连支:\n");
for(i=0;i<Ntr;i++)
fprintf(fp,"%d,",Link[i]);
fprintf(fp,"\n");
fprintf(fp,"\n相应的关联矩阵A(树支在前，连支在后):\n");
for(i=0;i<Ntr;i++)
{
for(j=0;j<Ntr;j++)
fprintf(fp,"%d\t",At[i][j]);
for(j=Ntr;j<Nbr;j++)
fprintf(fp,"%d\t",Al[i][j-Ntr]);
fprintf(fp,"\n");
}
fprintf(fp,"\n基本回路矩阵Bf(树支在前，连支在后):\n");
for(i=0;i<Ntr;i++)
{
for(j=0;j<Nbr;j++)
fprintf(fp,"%d\t",Bf[i][j]);
fprintf(fp,"\n");
}
fprintf(fp,"\n基本割集矩阵Qf(树支在前，连支在后):\n");
for(i=0;i<Ntr;i++)
{
for(j=0;j<Nbr;j++)
fprintf(fp,"%d\t",Qf[i][j]);
fprintf(fp,"\n");
}
fclose(fp);
}
