sparsem.cpp

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

void CSparseMatrix::del_c_f(INT2 m, LINKED_LIST &list, INT2 *fill)
/* [ 功  能 ]   在双链接表 list 所描述的稀疏矩阵中, 删去第 m 主列, 对其它候
 *              选主元产生的填元数 fill[] 的影响
 * =======================================================================*/
{
    INT2 mk,q,qk,j;

    mk=list.cp[m];						// mk 为第 m 列元素号 
    while(mk!=0){						// 搜索第 m 列非零元素 
        if(list.del[mk]==0) {			// 未删元素 
            q=list.row[mk];				// a[q][m] 为非零元素 
            qk=list.cp[q];				// qk 为第 q 列元素号 
            while(qk!=0){				// 搜索第 q 列非零元素 
                if(list.del[qk]==0) {   // 未删元素 
                    j=list.row[qk];		// a[j][q] 为非零元素 
                    /* 判断 a[j][m] 是否是零。若是零, 则候选主元 q 通过 
                       a[q][m] 和 a[j][q] 产生的填元 a[j][m] 将由于 m 列的
                       删除而不存在, fill[q] - 1。*/
                    if(ele_num(j,m,list)==0) fill[q]--; 
                }
                qk=list.down[qk];
            } 
        }
        mk=list.down[mk];
    }
};
// =====================================================================================
void CSparseMatrix::del_rc_p(INT2 mi, INT2 mj, LINKED_LIST &list, INT2 *nzr, INT2 *nzc)
/* 在双链接表 list 所描述的稀疏矩阵中, 删去第 mi 行 / mj 列, 对其它候选主元
   产生的列量/行量数组 nzc[]/nzr[] 的影响
   ====================================================================================*/
{
    INT2 mk,q;
   
	// 搜索第 mi 行非零元素 
    mk=list.rp[mi];					// mk 为第 mi 行元素号 
    while(mk!=0){					// 搜索第 mi 行非零元素 
        if(list.del[mk]==0) {		// 未删元素 
            q=list.col[mk];			// a[m][q] 为非零元素 
            if(q!=mj) nzc[q]--;
        }
        mk=list.right[mk];
    }
    // 搜索第 mj 列非零元素 
    mk=list.cp[mj];					// mk 为第 mj 列元素号 
    while(mk!=0){					// 搜索第 mj 行非零元素 
        if(list.del[mk]==0) {		// 未删元素 
            q=list.row[mk];			// a[q][m] 为非零元素 
            if(q!=mi) nzr[q]--;
        }
        mk=list.down[mk];
    }
};
//==========================================================================
INT2 CSparseMatrix::tinney_walker(LINKED_LIST &list, PRC *prc)
/* [ 功  能 ]   用 Tinney_Walker 算法排序双链接表 list 所描述的稀疏矩阵, 排
 *              序后的行/列顺序号放在数组 prc 中, 返回填元数。
 * =======================================================================*/
{
    INT2 n;						// 矩阵阶数 
    INT2 m;						// 第 k 步主行/列号 
    INT2 fm,fmm;				// m 主元时的填元数,及总填元数 
    INT2 i,j,k,q,x,y;
    INT2 nza,non;
    INT2 *fill;					// 各步填元数
    INT2 *fill_row,*fill_col;	// 存储生成的填元的行和列号 

    non=maxNonzero/2;
	ASSERT (non);
    fill_row=new INT2[non];
    fill_col=new INT2[non];
    fill=new INT2[maxRank+1];
    
	n=list.n;
    for(i=1;i<n+1;i++) prc[i]=i; // prc[] 原始排列 
    nza=list.nza;
    for(i=1;i<nza+1;i++) list.del[i]=0;
    fmm=0;
    for(i=1;i<n+1;i++) fill_in(i,list,fill,fill_row,fill_col,0);
	// 第一步,计算各候选主元对应的填元数, 记录在 fill[] 中 
    for(k=1;k<n;k++){					// k=1~n-1, 确定每一步的主元 m, 处理填元和主行/列 
        fillmin(n,k,prc,fill);			// 在 fill[] 中寻找最小填元数, 确定主元 m 
        m=prc[k];						// 第 k 步主元号 
        del_r_f(m,list,fill);			// 计入删去主行 m 对 fill[q] 的影响 
        del_c_f(m,list,fill);			// 计入删去主列 m 对 fill[q] 的影响 
		fm=fill[m];
        fmm=fmm+fm;
		if(fm!=0){
			fill_in(m,list,fill,fill_row,fill_col,2);
			// 主元 m 有填元需处理, fm 个填元的行/列标在 fill_row[]/fill_col[] 中 
			for(i=1;i<fm+1;i++){			// m 的填元对 fill[] 的影响 
				x=fill_row[i];
				y=fill_col[i];
				fillin_f(x,y,list,fill);
				// 增加填元 a[x][y] 后, 对候选主元 x 和 y 填元数的影响 
				// 其它候选主元 q 也以 a[x][y] 为主元时, fill[q]-1 
				for(j=k+1;j<n+1;j++){
					q=prc[j];
					if(q!=x){		// q=x 时, 候选主元 q 不会以 a[x][y] 为主元 
						if(q!=y){   // q=y 时, 候选主元 q 不会以 a[x][y] 为主元 
							if(ele_num(x,q,list)!=0 && ele_num(q,y,list)!=0){
								fill[q]--;
							}                            
						}
					}
				}
				insert_ele(x,y,list);
			}
		}
		del_row(m,list);    // 删去第 m 行 
		del_col(m,list);    // 删去第 m 列 
    }

    delete fill_row;
    delete fill_col;
    delete fill;

    return(fmm);
};
//=========================================================================
INT2 CSparseMatrix::markowitz(LINKED_LIST &list,PRC *prc)
/* [ 功  能 ]   用 Markowitz 算法排序双链接表 list 所描述的稀疏矩阵, 排
 *              序后的行/列顺序号放在数组 prc 中, 返回长运算数。
 * =======================================================================*/
{
    INT2 n;						// 矩阵阶数 
    INT2 m;						// 第 k 步主行/列号
    INT2 lopmin,lopmm;			// m 主元时的 LOP 数,及总长运算数 
    INT2 *nzr;					// 行量数组 
    INT2 *nzc;					// 列量数组 
    INT2 i,k,x,y;
    INT2 nza,fm;
    INT2 non;
    INT2 *fill_row,*fill_col;	// 存储生成的填元的行和列号 
    INT2 *fill,*lop;

    nzr=new INT2[maxRank+1];
    nzc=new INT2[maxRank+1];
    fill=new INT2[maxRank+1];
    lop=new INT2[maxRank+1];
    non=maxNonzero/2;
	ASSERT (non);
    fill_row=new INT2[non];
    fill_col=new INT2[non];

    n=list.n;
    for(i=1;i<n+1;i++) prc[i]=i;	// prc[] 原始排列 
    nza=list.nza;
    for(i=1;i<nza+1;i++) list.del[i]=0;
    nzerorc(list,nzr,nzc);			// 计算原始行量/列量 
    lopmm=0;
    for(k=1;k<n;k++){				// k=1~n-1, 确定每一步的主元 m, 处理填元和主行/列 
		lop_min(n,k,prc,nzr,nzc,&lopmin);
		m=prc[k];					// 第 k 步主元号
        lop[m]=lopmin;
        lopmm=lopmm+lopmin;
		del_rc_p(m,list,nzr,nzc);	// 主行/列删去对列/行量的影响 
        if(lopmin>1){
            fill_in(m,list,fill,fill_row,fill_col,3);   // 处理主元 m 的填入
            fm=fill[m];
			if(fm!=0){
				for(i=1;i<fm+1;i++){    // m 的填元 a[x][y] 对行/列量的影响
					x=fill_row[i];
                    y=fill_col[i];
                    nzr[x]++;
                    nzc[y]++;
                }
            }
        }
        del_row(m,list);				// 删去第 m 行
        del_col(m,list);				// 删去第 m 列 
    }

    delete nzr;
    delete nzc;
    delete fill_row;
    delete fill_col;
    delete fill;
    delete lop;

    return(lopmm);
};
// ==========================================================================
INT2 CSparseMatrix::markowitz(LINKED_LIST &list,PRC *prow,
    PRC *pcol)
/*	用 Markowitz 算法排序双链接表 list 所描述的稀疏矩阵, 以下三角矩阵中的非
	零元作为候选主元, 排序后的行/列顺序号放在 prow[]/pcol[] 中, 返回长运算数
	========================================================================*/
{
    INT2 n;						// 矩阵阶数 
    INT2 mi,mj;					// 第 k 步主行/列号 
    INT2 lopmin,lopmm;			// m 主元时的 LOP 数及总长运算数 
    INT2 *nzr;					// 行量数组 
    INT2 *nzc;					// 列量数组 
    INT2 i,k,x,y;
    INT2 nza,fm;
    INT2 non;
    INT2 *fill_row,*fill_col;	// 存储生成的填元的行和列号 
    INT2 *fill,*lop;

    non=maxNonzero/2;
    nzr=new INT2[maxRank+1];
    nzc=new INT2[maxRank+1];
    fill=new INT2[maxRank+1];
    lop=new INT2[maxRank+1];
	ASSERT (non);
    fill_row=new INT2[non];
    fill_col=new INT2[non];
    n=list.n;
    for(i=1;i<n+1;i++){			// 原始排列 
        prow[i]=i; 
        pcol[i]=i;
        }
    nza=list.nza;
    for(i=1;i<nza+1;i++) list.del[i]=0;
    nzerorc(list,nzr,nzc);    // 计算原始行量/列量 
    lopmm=0;
    pcol[n]=n;
    for(k=1;k<n;k++){		  // k=1~n-1, 确定每一步的主元 mi/mj, 处理填元和主行/列 
		lop_min(k,list,prow,pcol,nzr,nzc,&lopmin);
		mi=prow[k];			  // 第 k 步主行号 
		mj=pcol[k];		      // 第 k 步主列号 
        pcol[n]=pcol[n]+k-mj;
        lop[k]=lopmin;
        lopmm=lopmm+lopmin;
		del_rc_p(mi,mj,list,nzr,nzc);  // 主行/列删去对列/行量的影响 
        if(lopmin>1){
			fill_in(k,mi,mj,list,fill,fill_row,fill_col);
			// 处理主元a[mi][mj] 的填入
            fm=fill[k];
			if(fm!=0){
				for(i=1;i<fm+1;i++){    //m 的填元 a[x][y] 对行/列量的影响
					x=fill_row[i];
                    y=fill_col[i];
                    nzr[x]++;
                    nzc[y]++;
                }
            }
        }
        del_row(mi,list);   // 删去第 m 行 
        del_col(mj,list);   // 删去第 m 列 
    }

    delete nzr;
    delete nzc;
    delete fill_row;
    delete fill_col;
    delete fill;
    delete lop;

    return(lopmm);
};
//==========================================================================
void CSparseMatrix::old_new(INT2 n, PRC *prc, OLD_NEW *old_newrc)
/* [ 功  能 ]   由优化后排序号, 即新=>旧映射 prc] 生成旧=>新映射 old_newrc[]
 * =======================================================================*/
{
    INT2 i,ii;

    for(i=1;i<n+1;i++){
        ii=prc[i];
        old_newrc[ii]=i;
    }
};
// ==========================================================================
void CSparseMatrix::old_new(INT2 n, PRC *prow, PRC *pcol, OLD_NEW *old_newr, 
    OLD_NEW *old_newc)
/* 由优化后排序号, 即新=>旧映射 prow[] / pcol[] 生成旧=>新映
   射 old_newr[] / old_newc[]
   =========================================================================*/
{
    INT2 i,ii,jj;

    for(i=1;i<n+1;i++){
        ii=prow[i];
        jj=pcol[i];
        old_newr[ii]=i;
		old_newc[jj]=i;
    }
};
//==========================================================================
void CSparseMatrix::list_table(LINKED_LIST &list, PRC *prc, 
    OLD_NEW *old_newrc, TRIANGULAR_TABLE &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=prc[i];						// 取原行号 ii 
        kr=list.rp[ii];
        sm=0;
        while(kr!=0){					// 扫描第 i / ii 行 
            y=list.col[kr];
            y=old_newrc[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=prc[i];				// 取原列号 jj 
        kc=list.cp[jj];
        sm=0;
        while(kc!=0){			// 扫描第 i / jj 列 
            x=list.row[kc];
            x=old_newrc[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 &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;  // 对角元素优先