利用稀疏矩阵的程序

#define R1
#define G2
#define VS3
#define CS4
#define VCCS5
#define VCVS6
#define CCCS7
#define CCVS8
#define OPAMP9
typedef struct tagBranch{
int number; //支路号
char type; //类型号
int nfrom; //始节点
int nto; //终节点
float value; //参数值
int ncfrom; //控制始节点/控制终节点；G/R支路的“设置为电压定义支路”标记int ncto; //控制终节点；G/R支路的“控制用电压定义支路”标记
} BRANCH;
BRANCH Branch0; *BRANCH;
Branch = new BRANCH[BranchNumber+1];
int BranchNumber;
int NodeNumber;
int vBranchNumber;
int NN;
double val[NN];
int row[NN];
int col[NN];
int up[NN];
int down[NN];
int left[NN];
int right[NN];
int rp[NN];
int cp[NN];
int del[NN];
int nza;
int n;
typedef struct tagLINKED_LIST{
int nza; //非零元素个数
int n; //矩阵阶数
int *row; //元素行号
int *col; //元素列号
int *up; //元素的上邻元素号
int *down; //元素的下邻元素号
int *left; //元素的左邻元素号
int *right; //元素的右邻元素号
int *rp; //行链指针
int *cp; //列链指针
char *del; //元素删除标志；1=删除
} LINKED_LIST;
typedef struct tagTRIANGULAR_TABLE{
int nza; //非零元素个数
int n; //矩阵阶数
double *val; //矩阵及右端向量/中间解向量元素值
int *roco; //对角元素行列号，U元素列号，L元素行号int *urp; //U阵行指针
int *lcp; //L阵列指针
int lun; //由符号LU分解确定的需修正的元素个数
int fen; //由符号前消确定的需修正的向量元素个数int *lup; //由符号LU分解确定的需修正的元素号
int *fep; //前消需修正的元素号
} TRIANGULAR_TABLE;
typedef struct tagTRIANGULAR_TABLE_C{
int nza;
int n;
int lun;
int fen;
double *vre;
double *vim;
int *roco;
int *urp;
int *lcp;
int *lup;
int *fep;
} tagTRIANGULAR_TABLE_C;
#include "XISHU.h"
//创建双重链接表
void new_list(int rank_max,int nonzero,LINKED_LIST *list) {
row=new int[nonzero];
col=new int[nonzero];
up=new int[nonzero];
down=new int[nonzero];
left=new int[nonzero];
right=new int[nonzero];
rp=new int[rank_max];
cp=new int[rank_max];
del=new int [nonzero];
nza=0; //创建时尚无元素
int i;
for(i=0;i<rank_max;i++){ //初始化
rp[i]=0;
cp[i]=0;
}
for(i=0;i<nonzero+1;i++){
row[i]=0;
col[i]=0;
up[i]=0;
down[i]=0;
left[i]=0;
right[i]=0;
del[i]=0;
}
}
//-----------------------------------------------------------------
char insert_row(int nza,int x,int y,LINKED_LIST *list)
{
int k,kl,z;
k=(*list).rp[x];
if(k==0){ //第x行还未有非零元素，a[x][y]为第一个元素(*list)->rp[x]=nza;
(*list)->left[nza]=0;
(*list)->right[nza]=0;
return(0);
} //成功插入
for(;;){ //搜索第x行非零元素
z=(*list).col[k]; //a[x][y]为非零元素
if(y<z){ //nza在k左
kl=(*list).left[k];
if(kl==0){ //k为原行首，nza在其左，为新行首
(*list)->left[nza]=0;
(*list)->right[nza]=k;
(*list)->left[k]=nza;
(*list)->rp[x]=nza;
}
else{ //nza在kl和k之间
(*list)->left[nza]=kl;
(*list)->right[nza]=k;
(*list)->left[k]=nza;
(*list)->right[k]=nza;
}
return(0);
}
else if(y>z){ //nza在k右
kl=(*list).right[k];
if(kl==0){ //k为原行尾，nza在其中，为新行尾
(*list)->right[nza]=0;
(*list)->right[k]=nza;
(*list)->left[nza]=k;
return(0);
}
k=kl; //右边还有非零元素，继续向右搜索
continue;
}
else{ //y=z,a[x][y]已存在
return(1); //已有，未插入
}
}
}
//----------------------------------------------------------------
char insert_col(int nza,int x,int y,LINKED_LIST *list)
{
int k,kl,z;
k=(*list).cp[y];
if(k==0){ //第y列还未有非零元素，a[x][y]为第一个元素(*list)->cp[y]=nza;
(*list)->up[nza]=0;
(*list)->down[nza]=0;
return(0);
} //成功插入
for(;;){ //搜索第y列非零元素
z=(*list).row[k]; //a[x][y]为非零元素
if(x<z){ //nza在k左
kl=(*list).up[k];
if(kl==0){ //k为原列首，nza在其上，为新列首
(*list)->up[nza]=0;
(*list)->down[nza]=k;
(*list)->up[k]=nza;
(*list)->cp[x]=nza;
}
else{ //nza在kl和k之间
(*list)->up[nza]=kl;
(*list)->down[nza]=k;
(*list)->up[k]=nza;
(*list)->down[k]=nza;
}
return(0);
}
else if(x>z){ //nza在k下
kl=(*list).down[k];
if(kl==0){ //k为原列尾，nza在其下，为新列尾
(*list)->down[nza]=0;
(*list)->down[k]=nza;
(*list)->up[nza]=k;
return(0);
}
k=kl; //下面还有非零元素，继续向下搜索
continue;
}
else{ //x=z,a[x][y]已存在
return(1); //已有，未插入
}
}
}
//----------------------------------------------------------------
//在双链接表list所描述的稀疏矩阵行链中，插入元素a[x][y]，返回其标
//号；如果未插入，返回0
int insert_ele(int x,int y,LINKED_LIST *list)
{
int k,kl,z,nza;
nza=(*list).nza;
nza++;
if(insert_row(nza,x,y,list)==1) return(0); //插入行链，如果已有，则不插入insert_col(nza,x,y,list); //插入列链
(*list).del[nza]=0;
(*list).row[nza]=x;
(*list).col[nza]=y;
(*list).nza=nza;
return(nza);
}
//----------------------------------------------------------------
//创建三角形表结构变量
void new_table(int rank_max,int nonzero,TRIANGULAR_TABLE *table)
{ //生成三角形表结构table,最大阶数为rank_max,非零元素数为nonzero int i,non;
(*table).val=new double[nonzero];
(*table)->roco=new int[nonzero];
(*table)->urp=new int[rank_max];
(*table)->lcp=new int[rank_max];
(*table)->lup=new int[nonzero];
non=rank_max*5; //一般前消需要修正的元素数不超过阶数的5倍
(*table)->fep=new int[non];
for(i=0;i<nonzero+1;i++) (*table)->val[i]=0;
}
//----------------------------------------------------------------
//双链表-三角形表转换
void list_tablel(LINKED_LIST *list,int *prow,int *pcol,int *old_newr
,int *old_newc,TRIANGULAR_TABLE *table)
{
int n; //矩阵阶数
int i,ii,jj,kr,kc,x,y,nn;
int nza;
int si,sj,sm,sim,st,sst,ss[500];
n=(*list).n;
(*table)->n=n;
(*table)->nza=(*list).nza;
for(i=1;i<n+1;i++) (*table)->roco[i]=i; //对角元素优先nza=n+1; //已建立n个非零元素
for(i=1;i<n+1;i++){
(*table)->urp[i]=nza;
ii=prow[i]; //取原行号ii
kr=(*list).rp[ii];
sm=0;
while(kr!0){ //扫描第i/ii行，将新列号保存在ss[]
y=(*list).col[kr];
y=old_newc[y]; //将旧列号转换为新列号
if(y>i){ //U行元素
ss[sm]=y;
sm++; //U行非零元个数
}
kr=(*list).right[kr];
}
if(sm>0){
for(si=0,;si<sm-1;si++){ //将U行非零元列号记入rocol[] sim=si;
sst=ss[si];
for(sj=si+1;sj<sm;sj++){ //按新列号排序
if(sst>ss[sj]){
sim=sj;
sst=ss[sj];
}
}
(*table)->roco[nza]=sst;
nza++;
ss[sim]=ss[si];
}
(*table)->roco[nza]=ss[sm-1];
nza++;
}
(*table)->lcp[i]=nza;
jj=pcol[i]; //取原列号jj
kc=(*list).cp[jj];
sm=0;
while(kc!0){ //扫描第i/jj列
x=(*list).row[kc];
x=old_newr[x]; //将旧行号转换为新行号
if(x>i){ //L列元素
ss[sm]=x;
sm++;
}
kc=(*list).down[kc];
}
if(sm>0){
for(si=0,;si<sm-1;si++){
sim=si;
sst=ss[si];
for(sj=si+1;sj<sm;sj++){
if(sst>ss[sj]){
sim=sj;
sst=ss[sj];
}
}
(*table)->roco[nza]=sst;
nza++;
ss[sim]=ss[si];
}
(*table)->roco[nza]=ss[sm-1];
nza++;
}
}
for(i=1;i<nza+1;i++) (*table)->val[i]=0;
}
//----------------------------------------------------------------
//符号LU分解，生成table.lup[]
void sym_lu(TRIANGULAR_TABLE *table)
{
int n; //矩阵阶数
int k,kk,kkk,r,rr,i,c,cc,j,m,t;
int lu; //LU分解修正元素号
n=(*table).n;
lu=0;
for(k=1;k<n;k++){ //符号分解第k步，主元为k号元素
r=(*table).lcp[k]; //主列指针
c=(*table).urp[k]; //主行指针
rr=(*table).urp[k+1]; //主列非零元素号从r到rr-1
cc=(*table).lcp[k]; //主行非零元素号从c到cc-1
for(kk=r;kk<rr;kk++){ //扫描主列非零元素
i=(*table).roco[kk]; //i为主列非零元素行号
for(kkk=c;kkk<cc;kkk++){ //扫描主行非零元素
j=(*table).roco[kkk]; //j为主行非零元素列号
if(i==j) m=j; //对角元素
else if(i<j){ //U行元素
m=(*table).urp[i];
while(j!=(*table)->roco [m]) m++; //扫描第iL行
}
else{ //L列元素
m=(*table).lcp[j];
while(i!=(*table)->roco[m]) m++; //扫描第jU列
}
lu++;
(*table)->lup[lu]=m;
}
}
}
(*table)->lun=lu; //LU分解需要修正的元素个数
}
//----------------------------------------------------------------
//数值LU分解
char num_lu(TRIANGULAR_TABLE *table)
{
int n; //矩阵阶数
int k,kk,kkk,r,rr,c,cc,t;
int lu; //LU分解修正元素号
double x,val;
n=(*table).n;
lu=0;
for(k=1;k<n;k++){ //符号分解第k步，主元为k号元素
x=(*table).val[k];
if(x==0) return(1);
x=1.0/x;
c=(*table).urp[k]; //主行指针
cc=(*table).lcp[k]; //k-U行非零元素号从c到cc-1
for(kk=c;kk<cc;kk++)
(*table)->val[kk]=(*table)->val[kk]*x; //归一化U主行r=(*table).lcp[k]; //主列指针
rr=(*table).urp[k+1]; //k-L列非零元素号从r到rr-1
for(kk=r;kk<rr;kk++){ //扫描k-L主列非零元素
for(kkk=c;kkk<cc;kkk++){ //扫描k-U主行非零元素
lu++;
t=(*table).lup[lu]; //修正元素号
(*table)->val[t]=(*table)->val[t]-
(*table)->val[kk]*(*table)->val[kkk];
}
}
}
return(0);
}
//----------------------------------------------------------------
//右端向量填元处理
void vector_fillin(TRIANGULAR_TABLE *table)
{
char f;
int n; //矩阵阶数
int b; //右端元素向量号
int nza; //非零元素个数
int k,kk,l,ll,i,r,p;
n=(*table).n;
nza=(*table).nza;
b=(*table).lcp[n]; //右端元素从lup[n]开始存放
while (b<nza+1) { //处理第b号右端向量非零元素
k=(*table).roco[b]; //向量行号
if(k==n){ //最后一个向量元素不会产生填元
(*table)->nza=nza;
return;
}
l=(*table).lcp[k]; //第k-L列
ll=(*table).urp[k+1]; //第k-L列非零元素号从l到ll-1
for(kk=l;kk<ll;kk++){ //扫描k-L主列,kk为其非零元素号
f=0;
i=(*table)->roco[kk]; //kk号右端向量非零元素行号
r=b+1; //b号的下一个右端向量非零元素号
while (r<nza+1) { //搜索b以下的右端向量非零元素r
if((*table)->roco[r]<i) r++;
else if((*table)->roco[r]==i){
f=1;
break; //已有，此kk非填元
}
else if((*table)->roco[r]>i){ //b和kk生成r之上的i行向量填元
nza++;
p=nza;
while (p>r) { //r及其后的向量元素后移一位
(*table)->roco[p]=(*table)->roco[p-1];
p--;
}
(*table)->roco[r]=i; //插入填元
f=1;
break;
}
}
if(f==0){
nza++;
(*table)->roco[nza]=i;
}
}
b++;
}
(*table)->nza=nza;
}
//----------------------------------------------------------------
//符号前消
void sym_fe(TRIANGULAR_TABLE *table)
{
int n; //矩阵阶数
int b; //右端元素向量号
int nza; //非零元素个数
int fe;
int k,kk,l,ll,i,r;
n=(*table).n;
nza=(*table).nza;
b=(*table).lcp[n]; //右端元素从lup[n]开始存放
fe=0; //fep[]标号
while (b<nza) { //处理第b号右端向量非零元素
k=(*table).roco[b]; //向量行号
l=(*table).lcp[k]; //第k-L列
ll=(*table).urp[k+1]; //第k-L列非零元素号从l到ll-1
r=b+1;
for(kk=l;kk<ll;kk++){ //扫描k-L主列,kk为其非零元素号i=(*table)->roco[kk]; //kk号右端向量非零元素行号
while (r<nza+1) { //搜索b以下的右端向量非零元素r
if((*table)->roco[r]<i){
fe++;
(*table)->fep[fe]=r;
r++;
break;
}
r++;
}
}
b++;
}
(*table)->fen=fe;
}
//----------------------------------------------------------------
//数值前消
char num_fe(TRIANGULAR_TABLE *table)
{
int n; //矩阵阶数
int b; //右端元素向量号
int nza; //非零元素个数
int fe;
int k,kk,kkk,l,ll;
double x;
n=(*table).n;
nza=(*table).nza;
b=(*table).lcp[n]; //右端元素从lup[n]开始存放
fe=0; //fep[]标号
while (b<nza) { //处理第b号右端向量非零元素
k=(*table).roco[b]; //向量行号
x=(*table).val[k];
if(fabs(x)<1.0e-30) return(1);
else x=1.0/x;
x=(*table).val[b]*x;
(*table).val[b]=x;
l=(*table).lcp[k]; //第k-L列
ll=(*table).urp[k+1]; //第k-L列非零元素号从l到ll-1
for(kk=l;kk<ll;kk++){ //扫描k-L主列,kk为其非零元素号fe++;
kkk=(*table).fep[fe];
(*table)->val[kkk]=(*table)->val[kkk]-
(*table)->val[kk]*x;
}
b++;
}
k=(*table).roco[nza]; //向量行号
x=(*table).val[k];
if(fabs(x)<1.0e-30) return(1);
else x=1.0/x;
(*table)->val[nza]=(*table)->val[nza]*x;
return(0);
}
//----------------------------------------------------------------
//后代程序
void num_bs(TRIANGULAR_TABLE *table,double *x)
{
int n; //矩阵阶数
int nza,nn;
int i,m,j,l,ll;
double xx;
n=(*table).n;
nza=(*table).nza;
l=(*table).lcp[n];
nn=n+1;
for(i=1;i<nn;i++) x[i]=0;
ll=nza+1;
for(m=l;m<ll;m++){ //将中间解向量值传送到解向量
i=(*table).roco[m];
x[i]=(*table)->val[m];
}
nn=n-1;
for(i=nn;i>0;i--){ //从后往前代入
xx=x[i];
l=(*table).urp[i];
ll=(*table).lcp[i];
for(m=l;m<ll;m++){ //扫描i-U行非零元素
j=(*table).roco[m];
xx=xx-(*table).val[m]*x[j];
}
x[i]=xx;
}
}
#include "linear.h"
#include "XISHU.h"
//输入网络结构数据，构成Branch[]数组
//----------------------------------------------------------------
//-----------------------------------------------------------------
//实现Branch[]数组按支路号排序
void sortBranch(BRANCH *Branch, int BranchNumber)//根据支路号对Branch数组排序{
int len,i,j,k,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;}
}
memcpy(&Branch0,Branch+i,len);
memcpy(Branch+i,Branch+m,len);
memcpy(Branch+m,&Branch0,len);
}
}
//----------------------------------------------------------------
//-----------------------------------------------------------------
//提取电压定义支路
void vFoud()//确定电压定义支路的条数bv，并顺序编号；
//同时建立一般支路编号bNUM和电压定义支路编//号vbNUM之间的影射关系{
int vBranchNumber; //电压定义支路表
int *vbNUM_bNUM; //vbNUM=〉bNUM 映射表
int *bNUM_vbNUM; //bNUM=〉vbNUM 映射表
int i,j,k,t,tt;
bNUM_vbNUM = new int[BranchNumber+1];
vBranchNumber=0;//初始化
for(i=1;i<BranchNumber+1;i++)//建立bNUM_vbNUM映射表
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){
//if(Branch[k].ncto!=1)
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){
//if(Branch[k].ncto!=1)
vBranchNumber++;
bNUM_vbNUM[k]=vBranchNumber;
Branch[k].ncto=1;
}
}
else if((tt==G||tt==R)&Branch[i].ncfrom==1){//指定为电压定义支路的R/G vBranchNumber++;
bNUM_vbNUM[i]=vBranchNumber;
}
}
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
vbNUM_bNUM = new int[vBranchNumber+1];
for(i=1;i<BranchNumber+1;i++){//建立vbNUM_bNUM映射表
j=bNUM_vbNUM[i];
vbNUM_bNUM[j]=i;
}
cout<<"\nvBranchNumber="<<vBranchNumber;
}
//-------------------------------------------------------------------
//--------------------------------------------------------------------
// 建立改进节电方程稀疏矩阵的双链表结构
#define RANK100
#define NONZERO1000
LINKED_LIST mna_list;
TRIANGULAR mna_triangular_table;
int prow[RANK]; //各步主行号，新=>旧
int pcol[RANK]; //各步主列号，新=>旧
int old_newr[RANK]; //旧行号=>新行号
int old_newc[RANK]; //旧列号=>新列号
double *x; //按优化排序后顺序存放的解数组指针
new_list(RANK,NONZERO,&mna_list); //生成改进节电方程矩阵的双链表结构mna_list.n=NodeNumber+vBranchNumber; //赋方程阶数
BRANCH *br;
int i,j,k,kk,nb,t,tt,nf,nt,ncf,nct;
nb=NodeNumber+vBranchNumber+1; //方程阶数+1
x=new_double[nb]; //按优化顺序后顺序存放的解数组
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; //控制终节点
//----------------------------------------------------------------
if(t==R){ //电阻支路
if(ncf==0&nct==0){ //非电压定义支路
if(nf==0) insert_ele(nt,nt,&mna_list);
else if(nt==0) insert_ele(nf,nf,&mna_list);
else {
insert_ele(nf,nf,&mna_list);
insert_ele(nt,nt,&mna_list);
insert_ele(nf,nt,&mna_list);
insert_ele(nt,nf,&mna_list);
}
}
else{ //电压定义电阻支路
k=br->Number;
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,nf,&mna_list);
insert_ele(k,nt,&mna_list);
insert_ele(k,k,&mna_list);
}
}
//---------------------------------------------------------------else if(t==G){ //电导支路
if(ncf==0&nct==0){ //非电压定义支路
if(nf==0) insert_ele(nt,nt,&mna_list);
else if(nt==0) insert_ele(nf,nf,&mna_list);
else {
insert_ele(nf,nf,&mna_list);
insert_ele(nt,nt,&mna_list);
insert_ele(nf,nt,&mna_list);
insert_ele(nt,nf,&mna_list);
}
}
else{ //电压定义电阻支路
k=br->Number;
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,nf,&mna_list);
insert_ele(k,nt,&mna_list);
insert_ele(k,k,&mna_list);
}
}
//------------------------------------------------------------------else if(t==VCCS){ //电压控制电流源支路
insert_ele(nf,ncf,&mna_list);
insert_ele(nf,nct,&mna_list);
insert_ele(nt,ncf,&mna_list);
insert_ele(nt,nct,&mna_list);
}
//----------------------------------------------------------------else if(t==CCCS){ //电流控制电流源支路
k=ncf; //控制支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
}
//----------------------------------------------------------------
else if(t==VCVS){ //电压控制电压源支路
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,nf,&mna_list);
insert_ele(k,nt,&mna_list);
insert_ele(k,ncf,&mna_list);
insert_ele(k,nct,&mna_list);
}
//----------------------------------------------------------------
else if(t==CCVS){ //电流控制电压源支路
kk=br->ncf; //控制支路号
kk=bNUM_vbNUM[kk];
kk=kk+NodeNumber;
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,nf,&mna_list);
insert_ele(k,nt,&mna_list);
insert_ele(k,kk,&mna_list);
}
//----------------------------------------------------------------
else if(t==OPAMP){ //理想运算放大器
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,ncf,&mna_list);
insert_ele(k,nct,&mna_list);
}
//----------------------------------------------------------------
else if(t==VS){ //独立电压源
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(nf,k,&mna_list);
insert_ele(nt,k,&mna_list);
insert_ele(k,nf,&mna_list);
insert_ele(k,nt,&mna_list);
}
}
//----------------------------------------------------------------
//----------------------------------------------------------------
//符号解
//即优化排序,双链表,三角形表转换,符号LU分解，右端向量填元处理和//符号前消
//利用稀疏矩阵库
int n;
n=mna_list;
markowitzl(RANK,NONZERO,&mna_list,prow,pcol); //Markowitz优化排序
old_new(n,prow,pcol,old_newr,old_newc); //建立排序前=>后行/列号映射
new_table(RANK,NONZERO,&mna_triangular_table); //创建三角形表
list_tablel(&mna_list,prow,pcol,old_newr,old_newc,&mna_triangular_table); //双链表-三角形表转换sym_lu(&mna_triangular_table); //符号LU分解
//符号插入右端向量非零元素
for(i=1;i<BranchNumber+1;i++){ //扫描电路/支路数组
br=Branch+i;
t=br->type; //支路类型
nf=br->nfrom; //始节点
nt=br->nto; //终节点
//----------------------------------------------------------------
if(t==CS){ //独立电流源
insert_ele(nf,&mna_triangular_table,old_newr);
insert_ele(nt,&mna_triangular_table,old_newr);
}
//----------------------------------------------------------------
else if(t==VS){ ////独立电压源
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
insert_ele(k,&mna_triangular_table,old_newr);
}
}
vector_fillin(&mna_triangular_table); //确定和插入右端向量填元
sym_fe(&mna_triangular_table); //符号前消
//----------------------------------------------------------------
//-----------------------------------------------------------------
//数值解
//装配数值，进行数值求解
//利用稀疏矩阵库
//----------------------------------------------------------------
//装配系数数值
int ni,nj,nt,nn;
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){
ni=old_newr[nt];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
}
else if(nt==0){
ni=old_newr[nf];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
}
else{
ni=old_newr[nf];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nt];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nf];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
ni=old_newr[nt];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
}
}
else{ //电压定义电阻支路
k=br->Number;
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
ni=old_newr[k];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0*val;
}
}
//----------------------------------------------------------------else if(t==G){ //电导支路
if(ncf==0&nct==0){ //非电压定义支路
if(nf==0){
ni=old_newr[nt];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
}
else if(nt==0){
ni=old_newr[nf];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
}
else{
ni=old_newr[nf];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nt];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nf];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
ni=old_newr[nt];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
}
}
else{ //电压定义电阻支路
k=br->Number;
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0/val;
}
}
//----------------------------------------------------------------else if(t==BCCS){ //电压控制电流源支路
ni=old_newr[nf];
nj=old_newc[ncf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nf];
nj=old_newc[nct];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
ni=old_newr[nt];
nj=old_newc[ncf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
ni=old_newr[nt];
nj=old_newc[nct];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
}
//----------------------------------------------------------------else if(t==CCCS){ //电流控制电流源支路
k=ncf; //控制支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=val;
ni=old_newr[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0*val;
}
//----------------------------------------------------------------else if(t==VCVS){ //电压控制电压源支路
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[ncf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0*val;
ni=old_newr[k];
nj=old_newc[nct];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=val;
}
//----------------------------------------------------------------else if(t==CCVS){ //电流控制电压源支路
kk=br->ncf; //控制支路号
kk=bNUM_vbNUM[kk];
kk=kk+NodeNumber;
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[kk];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0*val;
}
//----------------------------------------------------------------else if(t==OPAMP){ //理想运算放大器
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[ncf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nj=old_newc[nct];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
}
//----------------------------------------------------------------else if(t==VS){ //独立电压源
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[nf];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
// aMatrix[nt][k]=-1.0;
ni=old_newr[nt];
nj=old_newc[k];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
ni=old_newr[k];
nj=old_newc[nf];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=1.0;
ni=old_newr[k];
nj=old_newc[nt];
nn=ele_num(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=-1.0;
}
}
//-----------------------------------------------------------------
//----------------------------------------------------------------
//数值LU分解
num_lu(&mna_triangular_table);
//----------------------------------------------------------------
//----------------------------------------------------------------
//装配右端向量数值
float val,vval;
BRANCH *br;
int i,j,k,kk,ni,t,tt,nf,nt,ncf,nct,nn;
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==CS){ //独立电流源
ni=old_newr[nf];
nn=ele_num2(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]+=val;
ni=old_newr[nt];
nn=ele_num2(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]-=val;
}
//----------------------------------------------------------------
else if(t==VS){ //独立电压源
k=br->Number; //支路号
k=bNUM_vbNUM[k]; //电压定义支路号
k=k+NodeNumber;
ni=old_newr[k];
nn=ele_num2(ni,nj,&mna_triangular_table);
mna_triangular_table.val[nn]=val;
}
}
//----------------------------------------------------------------
//-----------------------------------------------------------------
//数值前消
num_fe(&mna_triangular_table);
//后代
num_bs(&mna_triangular_table);。