sparsem.cpp

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

    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 行 
            y=list.col[kr];
            y=old_newc[y];				 // 将旧列号转换为新列号 
            if(y>i){					 // U 行元素 
                ss[sm]=y;
                sm++;
            }
            kr=list.right[kr];
        }
        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++;
        }
        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;
};
//=========================================================================
void CSparseMatrix::list_table(LINKED_LIST &list, PRC *prow, PRC *pcol, 
    OLD_NEW *old_newr, OLD_NEW *old_newc, TRIANGULAR_TABLE_C &table)
// [ 功  能 ]   将稀疏矩阵的双链接表 list 格式转换成三角形表 table。 
{
    INT2 n;							// 矩阵阶数 
    INT2 i,ii,jj,kr,kc,x,y;
    INT2 nza;
    INT2 si,sj,sm,sim,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 行 
            y=list.col[kr];
            y=old_newc[y];				 // 将旧列号转换为新列号 
            if(y>i){					 // U 行元素 
                ss[sm]=y;
                sm++;
            }
            kr=list.right[kr];
        }
        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++;
        }
        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.vre[i]=0;
        table.vim[i]=0;
    }
}
// =========================================================================
INT2 CSparseMatrix::ele_num(INT2 x, INT2 y, TRIANGULAR_TABLE &table)
// 返回 a[x][y] 在三角形表 table 中的标号。若是零元素, 返回 0。
{
    INT2 n,l,ll,k,z;

    n=table.n;
    if(x==0 || y==0 || x>n || y>n) return(0);
    if(x==y) return(x);						// 对角元素 
    else if(x>y){							// L-列元素, 搜索 y-L-列 
        l=table.lcp[y];
        ll=table.urp[y+1];
        for(k=l;k<ll;k++){
            z=table.roco[k];
            if(z==x) return(k);
        }
    }
    else{									// x<y, U-行 
        l=table.urp[x];
        ll=table.lcp[x];
        for(k=l;k<ll;k++){
            z=table.roco[k];
            if(z==y) return(k);
        }
    }
    return(0);
};
//=========================================================================
INT2 CSparseMatrix::ele_num(INT2 x, INT2 y, TRIANGULAR_TABLE_C &table)
// 返回 a[x][y] 在三角形表 table 中的标号。若是零元素, 返回 0。
{
    INT2 n,l,ll,k,z;

    n=table.n;
    if(x==0 || y==0 || x>n || y>n) return(0);
    if(x==y) return(x); // 对角元素 
    else if(x>y){		// L-列元素, 搜索 y-L-列 
        l=table.lcp[y];
        ll=table.urp[y+1];
        for(k=l;k<ll;k++){
            z=table.roco[k];
            if(z==x) return(k);
        }
    }
    else{				// x<y, U-行 
        l=table.urp[x];
        ll=table.lcp[x];
        for(k=l;k<ll;k++){
            z=table.roco[k];
            if(z==y) return(k);
        }
    }
    return(0);
};
// ==========================================================================
INT2 CSparseMatrix::ele_num(INT2 x, TRIANGULAR_TABLE &table)
// 返回 b[x] 在三角形表 table 中的标号。若是零元素, 返回 0。
{
    INT2 n,l,ll,k,z;

    n=table.n;
    if(x==0 || x>n) return(0);
    l=table.lcp[n];
    ll=table.nza+1;
    for(k=l;k<ll;k++){
        z=table.roco[k];
        if(z==x) return(k);
    }
    return(0);
};
// =========================================================================
INT2 CSparseMatrix::ele_num(INT2 x, TRIANGULAR_TABLE_C &table)
// 返回 b[x] 在三角形表 table 中的标号。若是零元素, 返回 0。
{
    INT2 n,l,ll,k,z;

    n=table.n;
    if(x==0 || x>n) return(0);
    l=table.lcp[n];
    ll=table.nza+1;
    for(k=l;k<ll;k++){
        z=table.roco[k];
        if(z==x) return(k);
    }
    return(0);
};
// =========================================================================
void CSparseMatrix::sym_lu(TRIANGULAR_TABLE &table)
/* 根据稀疏矩阵的三角形表 table 进行符号 LU 分解, 确定每一步分解所需要修正
   的元素编号, 即生成 table.lup[]。*/
{
    INT2 n;							// 矩阵阶数 
    INT2 k,kk,kkk,r,rr,i,c,cc,j,m;
    INT2 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++; // 扫描第 i L-行 
                }
                else{							 // L 列元素 
                    m=table.lcp[j];
                    while(i!=table.roco[m]) m++; //扫描第 j 列 
                }
                lu++;
                table.lup[lu]=m;
            }
        }
    }
    table.lun=lu;								// LU 分解需要修正的元素个数
};
// =========================================================================
void CSparseMatrix::sym_lu(TRIANGULAR_TABLE_C &table)
/* 根据稀疏矩阵的三角形表 table 进行符号 LU 分解, 确定每一步分解所需要修正
   的元素编号, 即生成 table.lup[]。*/
{
    INT2 n;								// 矩阵阶数 
    INT2 k,kk,kkk,r,rr,i,c,cc,j,m;
    INT2 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++; // 扫描第 i L-行 
                    }
                else{					// L 列元素 
                    m=table.lcp[j];
                    while(i!=table.roco[m]) m++; // 扫描第 j 列 
                    }
                lu++;
                table.lup[lu]=m;
                }
            }
        }
    table.lun=lu;    // LU 分解需要修正的元素个数 
};
// =========================================================================
char CSparseMatrix::div_c(double are, double aim, double bre, double bim, double *cre,
    double *cim)
// 复数除法: (*cre,*cim) = (are,aim) / (bre,bim)
{
    double r;

    r=bre*bre+bim*bim;
    if(r<1.0e-30) return(1);
    *cre=(are*bre+aim*bim)/r;
    *cim=(aim*bre-are*bim)/r;
    return(0);
};
// ===================================================================================
void CSparseMatrix::mul_c(double are, double aim, double bre, double bim, double *cre,
    double *cim)
// 复数乘法: (*cre,*cim) = (are,aim) * (bre,bim)
{
    *cre=are*bre-aim*bim;
    *cim=aim*bre+are*bim;
};
// ===========================================================================
void CSparseMatrix::xy_rp(double x, double y, double *r, double *p)
// 复数直角表示 (x,y) => 极坐标表示 (r,p)
{
    *r=sqrt(x*x+y*y);
    if(y==0) *p=0;
    else if(x==0) *p=PI/2.0;
    else *p=atan2(y,x);
};
// ===========================================================================
char CSparseMatrix::num_lu(TRIANGULAR_TABLE &table)
// 根据稀疏矩阵的三角形表 table 进行数值 LU 分解。
{
    INT2 n;						// 矩阵阶数  
    INT2 k,kk,kkk,r,rr,c,cc,t;
    INT2 lu;					// LU 分解修正元素号 
    double x;

    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];
				if((table.roco[kk] ==  table.roco[kkk])||(t<=n))
				{
					if(fabs(table.val[t])==0)
					{
						 table.val[t] = 1;
					}
					else if(fabs(table.val[t])<=1.0e-12) 
					{
						table.val[t]=table.val[t]/fabs(table.val[t]);
					}
					else
					{