利用稀疏矩阵的程序Word文档下载推荐.docx
《利用稀疏矩阵的程序Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《利用稀疏矩阵的程序Word文档下载推荐.docx（38页珍藏版）》请在冰豆网上搜索。
intdel[NN];
intnza;
intn;
typedefstructtagLINKED_LIST{
intnza;
//非零元素个数
intn;
//矩阵阶数
int*row;
//元素行号
int*col;
//元素列号
int*up;
//元素的上邻元素号
int*down;
//元素的下邻元素号
int*left;
//元素的左邻元素号
int*right;
//元素的右邻元素号
int*rp;
//行链指针
int*cp;
//列链指针
char*del;
//元素删除标志；
1=删除
}LINKED_LIST;
typedefstructtagTRIANGULAR_TABLE{
double*val;
//矩阵及右端向量/中间解向量元素值
int*roco;
//对角元素行列号，U元素列号，L元素行号
int*urp;
//U阵行指针
int*lcp;
//L阵列指针
intlun;
//由符号LU分解确定的需修正的元素个数
intfen;
//由符号前消确定的需修正的向量元素个数
int*lup;
//由符号LU分解确定的需修正的元素号
int*fep;
//前消需修正的元素号
}TRIANGULAR_TABLE;
typedefstructtagTRIANGULAR_TABLE_C{
double*vre;
double*vim;
}tagTRIANGULAR_TABLE_C;
#include"
XISHU.h"
//创建双重链接表
voidnew_list（intrank_max,intnonzero,LINKED_LIST*list）
{
row=newint[nonzero];
col=newint[nonzero];
up=newint[nonzero];
down=newint[nonzero];
left=newint[nonzero];
right=newint[nonzero];
rp=newint[rank_max];
cp=newint[rank_max];
del=newint[nonzero];
nza=0;
//创建时尚无元素
inti;
for（i=0;
i<
rank_max;
i++）{//初始化
rp[i]=0;
cp[i]=0;
}
nonzero+1;
i++）{
row[i]=0;
col[i]=0;
up[i]=0;
down[i]=0;
left[i]=0;
right[i]=0;
del[i]=0;
}
//-----------------------------------------------------------------
charinsert_row（intnza,intx,inty,LINKED_LIST*list）
intk,kl,z;
k=（*list）.rp[x];
if（k==0）{//第x行还未有非零元素，a[x][y]为第一个元素
（*list）->
rp[x]=nza;
left[nza]=0;
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）->
right[nza]=k;
left[k]=nza;
}
else{//nza在kl和k之间
left[nza]=kl;
right[k]=nza;
return（0）;
elseif（y>
z）{//nza在k右
kl=（*list）.right[k];
if（kl==0）{//k为原行尾，nza在其中，为新行尾
left[nza]=k;
return（0）;
k=kl;
//右边还有非零元素，继续向右搜索
continue;
else{//y=z,a[x][y]已存在
return
（1）;
//已有，未插入
//----------------------------------------------------------------
charinsert_col（intnza,intx,inty,LINKED_LIST*list）
k=（*list）.cp[y];
if（k==0）{//第y列还未有非零元素，a[x][y]为第一个元素
（*list）->
cp[y]=nza;
up[nza]=0;
down[nza]=0;
）{//搜索第y列非零元素
z=（*list）.row[k];
if（x<
kl=（*list）.up[k];
if（kl==0）{//k为原列首，nza在其上，为新列首
（*list）->
down[nza]=k;
up[k]=nza;
cp[x]=nza;
}
else{//nza在kl和k之间
up[nza]=kl;
down[k]=nza;
elseif（x>
z）{//nza在k下
kl=（*list）.down[k];
if（kl==0）{//k为原列尾，nza在其下，为新列尾
up[nza]=k;
return（0）;
k=kl;
//下面还有非零元素，继续向下搜索
continue;
else{//x=z,a[x][y]已存在
return
（1）;
}
//在双链接表list所描述的稀疏矩阵行链中，插入元素a[x][y]，返回其标
//号；
如果未插入，返回0
intinsert_ele（intx,inty,LINKED_LIST*list）
intk,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）;
//创建三角形表结构变量
voidnew_table（intrank_max,intnonzero,TRIANGULAR_TABLE*table）
{//生成三角形表结构table,最大阶数为rank_max,非零元素数为nonzero
inti,non;
（*table）.val=newdouble[nonzero];
（*table）->
roco=newint[nonzero];
（*table）->
urp=newint[rank_max];
lcp=newint[rank_max];
lup=newint[nonzero];
non=rank_max*5;
//一般前消需要修正的元素数不超过阶数的5倍
fep=newint[non];
for（i=0;
i++）（*table）->
val[i]=0;
//双链表-三角形表转换
voidlist_tablel（LINKED_LIST*list,int*prow,int*pcol,int*old_newr
int*old_newc,TRIANGULAR_TABLE*table）
inti,ii,jj,kr,kc,x,y,nn;
intsi,sj,sm,sim,st,sst,ss[500];
n=（*list）.n;
n=n;
nza=（*list）.nza;
for（i=1;
n+1;
roco[i]=i;
//对角元素优先
nza=n+1;
//已建立n个非零元素
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++;
lcp[i]=nza;
jj=pcol[i];
//取原列号jj
kc=（*list）.cp[jj];
while（kc!
0）{//扫描第i/jj列
x=（*list）.row[kc];
x=old_newr[x];
//将旧行号转换为新行号
if（x>
i）{//L列元素
ss[sm]=x;
kc=（*list）.down[kc];
si++）{
sj++）{
nza+1;
//符号LU分解，生成table.lup[]
voidsym_lu（TRIANGULAR_TABLE*table）
intk,kk,kkk,r,rr,i,c,cc,j,m,t;
intlu;
//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;
//对角元素
elseif（i<
j）{//U行元素
m=（*table）.urp[i];
while（j!
=（*table）->
roco[m]）m++;
//扫描第iL行
else{//L列元素
m=（*table）.lcp[j];
while（i!
roco[m]）m++;
//扫描第jU列
lu++;
（*table）->
lup[lu]=m;
}
lun=lu;
//LU分解需要修正的元素个数
//数值LU分解
charnum_lu（TRIANGULAR_TABLE*table）
intk,kk,kkk,r,rr,c,cc,t;
intlu;
doublex,val;
lu=0;
for（k=1;
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++）
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++）{//扫描k-L主列非零元素
for（kkk=c;
kkk++）{//扫描k-U主行非零元素
lu++;
t=（*table）.lup[lu];
//修正元素号
val[t]=（*table）->
val[t]-
val[kk]*（*table）->
val[kkk];
//右端向量填元处理
voidvector_fillin（TRIANGULAR_TABLE*table）
charf;
intb;
//右端元素向量号
intk,kk,l,ll,i,r,p;
nza=（*table）.nza;
b=（*table）.lcp[n];
//右端元素从lup[n]开始存放
while（b<
nza+1）{//处理第b号右端向量非零元素
k=（*table）.roco[b];
//向量行号
if（k==n）{//最后一个向量元素不会产生填元
nza=nza;
return;
l=（*table）.lcp[k];
//第k-L列
ll=（*table）.urp[k+1];
//第k-L列非零元素号从l到ll-1
for（kk=l;
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++;
elseif（（*table）->
roco[r]==i）{
f=1;
break;
//已有，此kk非填元
roco[r]>
i）{//b和kk生成r之上的i行向量填元
nza++;
p=nza;
while（p>
r）{//r及其后的向量元素后移一位
（*table）->
roco[p]=（*table）->
roco[p-1];
p--;
roco[r]=i;
//插入填元
if（f==0）{
roco[nza]=i;
b++;
//符号前消
voidsym_fe（TRIANGULAR_TABLE*table）
intfe;
intk,kk,l,ll,i,r;
fe=0;
//fep[]标号
nza）{//处理第b号右端向量非零元素
r=b+1;
ro
升级会员 