fullmatrix_analysis.c.h

是稀疏矩阵库以及全矩阵和稀疏矩阵分析程序源代码

#define R     1
#define G     2
#define VS    3
#define CS    4
#define VCCS  5
#define VCVS  6
#define CCCS  7
#define CCVS  8
#define OPAMP 9
#define N     50

int BranchNumber;
int NodeNumber;
int vBranchNumber;
int *vbNUM_bNUM;
int *bNUM_vbNUM;
float aMatrix[N+1][N+2];
float a[N+1][N+2];
float p,absp,b,t1;
int i,j,k,m,t,n,tt,kk,nt,nb,nf,ncf,nct;

int x0,y0,y00,fg=0;

float B[10][50];
float x[N],Vb[N],Ib[N];

//----------------------------------------------------
//支路结构
typedef struct tagbranch
{

	int number;
	int type;
	int nfrom;
	int nto;
	float value;
	int ncfrom;
	int ncto;
} BRANCH;

BRANCH Branch0,*Branch;


void sortBranch(BRANCH *,int);
void vbBranch();
void allMatrix();
int GaussFE();
int GaussBS();


//------------------------------------------------------
//支路排序
void  sortBranch(BRANCH *Branch,int BranchNumber)
{
	int len,l,m;
	len=sizeof(Branch0);
	for(i=1;i<BranchNumber+1;i++)
	{
		k=Branch[i].number;
		m=i;
		for(j=i+1;j<BranchNumber+1;j++)
		{
			l=Branch[j].number;
			if(k>l){ k=l; m=j; }
		}
		memmove(&Branch0,Branch+i,len);
		memmove(Branch+i,Branch+m,len);
		memmove(Branch+m,&Branch0,len);
	}
}

//-----------------------------------------------------
//扫描电压定义支路
void vbBranch()
{
	bNUM_vbNUM=new int[BranchNumber+1];
	vBranchNumber=0;
	for(int i=1;i<BranchNumber+1;i++)
		bNUM_vbNUM[i]=0;
	for(i=1;i<BranchNumber+1;i++)
	{ 
		t=Branch[i].type;
		if(t==VS||t==VCVS||t==OPAMP)
			｛
			vBranchNumber++;
		bNUM_vbNUM[i]=vBranchNumber;
	}
		 else if(t==CCVS)
		 {
			 vBranchNumber++;
			 bNUM_vbNUM[i]=vBranchNumber;
			 k=Branch[i].ncfrom;
			 tt=Branch[k].type;
			 if((tt==G||tt==R)&&Branch[k].ncfrom==0&&Branch[k].ncto==0)
			 { 
				 vBranchNumber++;
				 bNUM_vbNUM[k]=vBranchNumber;
				 Branch[k].ncto=1;
			 }
		 }
		 else if(t==CCCS)
		 { 
			 k=Branch[i].ncfrom;
			 tt=Branch[k].type;
			 if((tt==G||tt==R)&&Branch[k].ncfrom==0&&Branch[k].ncto==0)
			 {
				 vBranchNumber++;
				 bNUM_vbNUM[k]=vBranchNumber;
				 Branch[k].ncto=1;
			 }
		 }
		 else if((t==G||t==R)&Branch[i].ncfrom==1)
		 {
			 vBranchNumber++;
			 bNUM_vbNUM[k]=vBranchNumber;
		 }
}
vbNUM_bNUM=new int[vBranchNumber+1];
for(i=1;i<BranchNumber+1;i++)
{
	j=bNUM_vbNUM[i];
	vbNUM_bNUM[j]=i;
}
 }

 //------------------------------------------------------------------------------------
 //改进节点法建立全矩阵
 void allMatrix()
 {  float val;
 BRANCH *br;
 nb=NodeNumber+vBranchNumber+1;
 for(i=1;i<nb+1;i++)
	 for(j=1;j<nb+2;j++)
		 aMatrix[i][j]=0;
 for(i=1;i<BranchNumber+1;i++)
 {       
	 br=Branch+i;
	 t=br->type;
	 nf=br->nfrom;
	 nt=br->nto;
	 ncf=br->ncfrom;
	 nct=br->ncto;
	 val=br->value;
	 if(t==R)
	 {
		 if(ncf==0&&nct==0)
		 {
			 val=1.0/val;
			 if(nf==0) aMatrix[nt][nt]+=val;
			 else if(nt==0)aMatrix[nf][nf]+=val;
			 else {
				 aMatrix[nf][nf]+=val;
				 aMatrix[nt][nt]+=val;
				 aMatrix[nf][nt]-=val;
				 aMatrix[nt][nf]-=val;
			 }
		 }
		 else{
			 k=br->number;
			 k=bNUM_vbNUM[k];
			 k=k+NodeNumber;
			 aMatrix[nf][k]+=1.0;
			 aMatrix[nt][k]+=-1.0;
			 aMatrix[k][nf]+=1.0;
			 aMatrix[k][nt]-=1.0;
			 aMatrix[k][k]-=1.0*val;
		 }
	 }
	 else if(t==G)
	 {
		 if(ncf==0&&nct==0)
		 {  if(nf==0) aMatrix[nt][nt]+=val;
		 else if(nt==0)aMatrix[nf][nf]+=val;
		 else {
			 aMatrix[nf][nf]+=val;
			 aMatrix[nt][nt]+=val;
			 aMatrix[nf][nt]-=val;
			 aMatrix[nt][nf]-=val;
		 }
		 }
		 else{
			 k=br->number;
			 k=bNUM_vbNUM[k];
			 k=k+NodeNumber;
			 aMatrix[nf][k]+=1.0;
			 aMatrix[nt][k]-=1.0;
			 aMatrix[k][nf]+=1.0;
			 aMatrix[k][nt]-=1.0;
			 aMatrix[k][k]-=1.0/val;
		 }
	 }
	 else if(t==VCCS)
	 {  aMatrix[nf][ncf]+=val;
	 aMatrix[nf][nct]-=val;
	 aMatrix[nt][ncf]-=val;
	 aMatrix[nt][nct]+=val;
	 }
	 else if(t==CCCS)
	 {  k=ncf;
	 k=bNUM_vbNUM[k];
	 k=k+NodeNumber;
	 aMatrix[nf][k]=val;
	 aMatrix[nt][k]-=1.0*val;
	 }
	 else if(t==VCVS)
	 {  k=br->number;
	 k=bNUM_vbNUM[k];
	 k=k+NodeNumber;
	 aMatrix[nf][k]+=1.0;
	 aMatrix[nt][k]-=1.0;
	 aMatrix[k][nf]+=1.0;
	 aMatrix[k][nt]-=1.0;
	 aMatrix[k][ncf]-=1.0*val;
	 aMatrix[k][nct]+=val;
	 }
	 else if(t==CCVS)
	 { kk=br->ncfrom;
	 kk=bNUM_vbNUM[kk];
	 kk=kk+NodeNumber;
	 k=br->number;
	 k=bNUM_vbNUM[k];
	 k=k+NodeNumber;
	 aMatrix[nf][k]+=1.0;
	 aMatrix[nt][k]-=1.0;
	 aMatrix[k][nf]+=1.0;
	 aMatrix[k][nt]-=1.0;
	 aMatrix[k][kk]-=1.0*val;
	 }
	 else if(t==OPAMP)
	 { k=br->number;
	 k=bNUM_vbNUM[k];
	 k=k+NodeNumber;
	 aMatrix[nf][k]+=1.0;
	 aMatrix[nt][k]-=1.0;
	 aMatrix[k][ncf]+=1.0;
	 aMatrix[k][nct]-=1.0;
	 }
	 else if(t==CS)
	 { aMatrix[nf][nb]-=val;//+ or -?
	 aMatrix[nt][nb]+=val;
	 }
	 else if(t==VS)
	 { k=br->number;
	 k=bNUM_vbNUM[k];
	 k=k+NodeNumber;
	 aMatrix[nf][k]+=1.0;
	 aMatrix[nt][k]-=1.0;
	 aMatrix[k][nf]+=1.0;
	 aMatrix[k][nt]-=1.0;
	 aMatrix[k][nb]+=val;
	 }
 }
 }

 //-------------------------------------------------------------------------
 //高斯前消
 int GaussFE()
 {
	 for(i=1;i<nb;i++)
	 { 
		 for(j=1;j<nb+1;j++)
			 a[i][j]=aMatrix[i][j];
	 }
	 n=nb-1;
	 float x[N];
	 for(k=1;k<n;k++)
	 { p=a[k][k];
	 absp=fabs(p);
	 m=k;
	 for(i=k+1;i<n+1;i++)
	 {
		 if (fabs(a[i][k])>absp)
		 { absp=fabs(a[i][k]);
		 m=i;
		 }
	 }
	 if(absp<1e-21) return 0;
	 if(m!=k)
		 for(j=1;j<n+2;j++)
		 {t1=a[m][j];a[m][j]=a[k][j];a[k][j]=t1;}
		 b=a[k][k];
		 for(j=k+1;j<n+2;j++)
			 a[k][j]=a[k][j]/b;
		 for(i=k+1;i<n+1;i++)
		 { b=a[i][k];
		 for(j=k+1;j<n+2;j++)
			 a[i][j]=a[i][j]-b*a[k][j];
		 }
	 }
	 a[n][n+1]=a[n][n+1]/a[n][n];
	 return (0);
 }

 //------------------------------------------------------------------
 //高斯后代
 int GaussBS()
 {
	 for(i=n;i>0;i--)
	 { x[i]=a[i][n+1];
	 for(j=i+1;j<n+1;j++)
		 x[i]=x[i]-a[i][j]*x[j];
	 ltoa(x[i],string[i],10);
	 }
	 x[0]=0;
	 PrintResult(x);
	 return 0;
 }