sparsem.cpp

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

        else if(x>z){		// nza 在 k 下
            k1=list.down[k];
            if(k1==0){		// k 为原列尾, nza 在其下, 为新列尾
                list.down[nza]=0;
                list.down[k]=nza;
                list.up[nza]=k;
                return(0);
            }
            k=k1;			// 下面还有非零元素, 继续向下搜索
            continue;
        }
        else{				// x=z, a[x][y] 已存在
            return(1);
        }
    }
};
//===========================================================================
INT2 CSparseMatrix::insert_ele(INT2 x, INT2 y, LINKED_LIST &list)
/* 在双链接表 list 所描述的稀疏矩阵行链中, 插入元素 a[x][y],返回其标号;
   如果未插入,返回 0*/
{
    INT2 nza;  

	nza=list.nza;
	//修改2003.3.29*****************
	INT2 n;
	n = list.n;
	if(x<1||x>n||y<1||y>n)
		return(0);
	//修改2003.3.29*****************
    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);
};
// ========================================================================
INT2 CSparseMatrix::ele_num(INT2 x, INT2 y, LINKED_LIST &list)
// 返回 a[x][y] 在双链接表 list 中的标号。若是零元素, 返回 0。
{
    INT2 k,z;

    k=list.rp[x];						// k 为第 x 行元素号 
    while(k!=0){						// 搜索第 x 行非零元素 
        if(list.del[k]==0) {			// 未删元素 
            z=list.col[k];
            if(z>y) break;  
            else if(z==y) return (k);   // a[x][y] 为非零元素 
        }
        k=list.right[k];
    }
    return (0);							// 未发现非零元素 a[x][y], 其为零元素 
};
//==========================================================================
void CSparseMatrix::del_row(INT2 m, LINKED_LIST &list)
// 在双链接表 list 所描述的稀疏矩阵中, 作上删去第 m 行的标记。
{
    INT2 mk;

    mk=list.rp[m];   // mk 为第 m 行元素号 
    while(mk!=0){
		list.del[mk]=1;
        mk=list.right[mk];
    }
};
//==========================================================================
void CSparseMatrix::del_col(INT2 m, LINKED_LIST &list)
// 在双链接表 list 所描述的稀疏矩阵中, 作上删去第 m 列的标记。
{
    INT2 mk;

    mk=list.cp[m];   // mk 为第 m 列元素号 
    while(mk!=0){
		list.del[mk]=1;
        mk=list.down[mk];
    }
};
//=========================================================================
void CSparseMatrix::nzerorc(LINKED_LIST &list, INT2 *nzr, INT2 *nzc)
/* 计算双链接表 list 所描述的稀疏矩阵各行/列量, 即非零元素个数,置于数
   组 nzr[]/nzc[]*/
{
    INT2 n;					// 矩阵阶数 
    INT2 i,k,kr,kc;

    n=list.n;
    for(i=1;i<n+1;i++){		// 计算原始的行量/列量 
        k=0;
        kr=list.rp[i];
        while(kr!=0){
            k++;
            kr=list.right[kr];
        }
        nzr[i]=k;
        k=0;
        kc=list.cp[i];
        while(kc!=0){
            k++;
            kc=list.down[kc];
        }
        nzc[i]=k;
    }
};
//=========================================================================
void CSparseMatrix::fillin_f(INT2 x, INT2 y, LINKED_LIST &list, INT2 *fill)
/* [ 功  能 ]   x 行/y 列增加填元 a[x][y] 后, 对候选主元 x/y 填元数的影响
 * =======================================================================*/
{
    INT2 i,j,kc,kr;  

    // 搜索 x 列非零元素 a[i][x] 
    kc=list.cp[x];					// kc 为第 x 列元素号
    while(kc!=0){					// 搜索第 x 列非零元素
        if(list.del[kc]==0) {		// 未删元素 
            i=list.row[kc];			// a[i][x] 为非零元素 
			if(i!=x && ele_num(i,y,list)==0) fill[x]++; // 原 a[i][y] 是零元素, 从而是填元 
		}
		kc=list.down[kc]; 
	}
    // 搜索 y 行非零元素 a[y][j] 
    kr=list.rp[y];					// kr 为第 y 行元素号
    while(kr!=0){					// 搜索第 y 行非零元素
        if(list.del[kr]==0) {		// 未删元素 
            j=list.col[kr];			// a[y][j] 为非零元素 
			if(j!=y && ele_num(x,j,list)==0) fill[y]++; // 原 a[x][j] 是零元素, 从而是填元 
        }
        kr=list.right[kr]; 
    }
};
//=================================================================================
void CSparseMatrix::fill_in(INT2 q, LINKED_LIST &list, INT2 *fill, INT2 *fill_row,
    INT2 *fill_col, char id)
/* [ 功  能 ]   处理双链接表 list 所描述的稀疏矩阵在 LU 分解过程中, 
 *              以 a[q][q] 为主元时产生的填元, 用于优化排序:	
 *              (1) 生成主元或候补主元 a[q][q] 产生的填元数 fill[q]
 *              (2) 将生成的填元行/列标存入 fill_row[]/fill_col[]
 *              (3) 确定非零元素总数 nza 的当前值, 及当前填元号
 *              (4) 插入填元, 修改双重链接表             
 *              id=1 时, 作 (1)。
 *              id=2 时, 作 (1),(2)。
 *              id=3 时, 作 (1),(2),(3),(4)。
 *=================================================================================*/
{
    INT2 kc,kr,x,y;  
    INT2 nza;

    nza=0;
    fill[q]=0;
    kc=list.cp[q];										// kc 为第 q 列元素号
    while(kc!=0){										// 搜索第 q 列非零元素
        if(list.del[kc]==0) {							// 未删元素 
            x=list.row[kc];								// a[x][q] 为非零元素 
            if(x!=q) {									// 非主元本身 
                kr=list.rp[q];							// kr 为第 q 行元素号
                while(kr!=0){							// 搜索第 q 行非零元素
                    if(list.del[kr]==0) {				// 未删元素 
                        y=list.col[kr];					// a[q][y] 为非零元素 
                        if(y!=q) {						// 非主元本身 
                            if(ele_num(x,y,list)==0){   // 原 a[x][y] 是零元素, 从而是填元 
                                nza++;
                                fill[q]++;
                                if(id>1){
                                    fill_row[nza]=x;
                                    fill_col[nza]=y;
                                    if(id>2) insert_ele(x,y,list);
                                }
                            }
                         }
                    }
                    kr=list.right[kr]; 
                }
            }
        }
        kc=list.down[kc]; 
    }
};
// ==================================================================================
void CSparseMatrix::fill_in(INT2 k, INT2 mi, INT2 mj, LINKED_LIST &list, INT2 *fill, 
    INT2 *fill_row, INT2 *fill_col)
/* 处理双链接表 list 所描述的稀疏矩阵在 LU 分解过程中, 第 k 步以 a[mi][mj] 为
   主元时产生的填元, 用于优化排序:	
   (1) 生成主元或候补主元 a[mi][mj] 产生的填元数 fill[k]
   (2) 将生成的填元行/列标存入 fill_row[]/fill_col[]
   (3) 插入填元, 修改双重链接表             
	==================================================================================
*/
{
    INT2 kc,kr,x,y;
    INT2 nza;

    nza=0;
    fill[k]=0;
    kc=list.cp[mj];										// kc 为第 mj 列元素号
    while(kc!=0){										// 搜索第 mj 列非零元素
        if(list.del[kc]==0) {							// 未删元素 
            x=list.row[kc];								// a[x][q] 为非零元素 
            if(x!=mi) {									// 非主元本身 
                kr=list.rp[mi];							// kr 为第 mi 行元素号
                while(kr!=0){							// 搜索第 mi 行非零元素
                    if(list.del[kr]==0) {				// 未删元素 
                        y=list.col[kr];					// a[q][y] 为非零元素 
                        if(y!=mj) {						// 非主元本身 
                            if(ele_num(x,y,list)==0){   // 原 a[x][y] 是零元素, 从而是填元 
                                nza++;
                                fill[k]++;
                                fill_row[nza]=x;
                                fill_col[nza]=y;
                                insert_ele(x,y,list);
                            }
                         }
                    }
                    kr=list.right[kr]; 
                }
            }
        }
        kc=list.down[kc]; 
    }
};
//==========================================================================
void CSparseMatrix::fillmin(INT2 n, INT2 k, PRC *prc, INT2 *fill)
/* [ 功  能 ]   根据填元数数组 fill[], 在 prci=k~n] 中找出填元数最小的候选主
 *              元号, 交换到 prck] 位置。
 *=========================================================================*/
{
    INT2 filmin;		// 最小填元数 
    INT2 i,j,im;

    // 在 fill[] 中寻找最小填元数, 确定主元 m 
    filmin=n*n;
    for(i=k;i<n+1;i++){
        j=prc[i];		// 取原始行/列号 
        if(fill[j]==0){
            im=i;
            break;
        }
        if(fill[j]<filmin){
            im=i;
            filmin=fill[j];
        }
    }            
    // 交换 prc[k] 和 prc[im] 
    i=prc[k];
    prc[k]=prc[im];
    prc[im]=i;
};
//========================================================================================
void CSparseMatrix::lop_min(INT2 n, INT2 k, PRC *prc, INT2 *nzr, INT2 *nzc, INT2 *lopmin)
/* [ 功  能 ]   根据行量/列量数组 nzr[]/nzc[] 计算  prci=k~n] 中各候选主元长
 *              运算量 LOP, 找出 LOP 数最小的候选主元号, 交换到 prck] 位置。
 * ======================================================================================*/
{
    INT2 i,j,im,lop;

    *lopmin=n*n;
    for(i=k;i<n+1;i++){
        j=prc[i];				// 取原始行/列号 
        lop=nzc[j]*(nzr[j]-1);  // 候选主元 j 对应的长运算数 
        if(lop==0){
            im=i;
            *lopmin=0;
            break;
        }
        else if(lop<*lopmin){
            im=i;
            *lopmin=lop;
        }
    }            
    // 交换 prck] 和 prcim] 
    i=prc[k];
    prc[k]=prc[im];
    prc[im]=i;
};
// ====================================================================================
void CSparseMatrix::lop_min(INT2 k, LINKED_LIST &list, PRC *prow, PRC *pcol, INT2 *nzr, 
    INT2 *nzc, INT2 *lopmin)
/* 根据行量/列量数组 nzr[]/nzc[] 计算  prow/pcol[i=k~n] 中各候选主元长运算
   量 LOP, 找出 LOP 数最小的候选主元号, 交换到 prow/pcol[k] 位置。
   ===================================================================================*/
{
    INT2 n;							// 矩阵阶数 
    INT2 i,ii,jj,im,jm,lop;
    INT2 nr,nc,kr;

    n=list.n;
    *lopmin=n*n;
    for(i=k;i<n+1;i++){
        ii=prow[i];					// 取原始行号 
        nr=nzr[ii]-1;
        kr=list.rp[ii];
        while(kr!=0){				// 搜索第 ii 行 
            if(list.del[kr]==0){	// 未删除元素 
                jj=list.col[kr]; 
                nc=nzc[jj];
                lop=nc*nr;			// 候选主元 a[ii][jj] 对应的长运算数 
                if(lop==0){
                    im=i;
                    jm=jj;
                    *lopmin=0;
                    break;
                }
                else if(lop<*lopmin){
                    im=i;
                    jm=jj;
                    *lopmin=lop;
                }
            }
            kr=list.right[kr];
        }
        if(*lopmin==0) break;
    }
    // 交换行/列 
    i=prow[k];
    prow[k]=prow[im];
    prow[im]=i;
    pcol[k]=jm;
};
//============================================================================
void CSparseMatrix::del_rc_p(INT2 m, LINKED_LIST &list, INT2 *nzr, INT2 *nzc)
/* [ 功  能 ]   在双链接表 list 所描述的稀疏矩阵中, 删去第 m 主行/列, 对其它
 *              候选主元产生的列量/行量数组 nzc[]/nzr[] 的影响
 * ==========================================================================*/
{
    INT2 mk,q;
    
	// 搜索第 m 行非零元素 
    mk=list.rp[m];					// mk 为第 m 行元素号 
    while(mk!=0){					// 搜索第 m 行非零元素 
        if(list.del[mk]==0) {		// 未删元素 
            q=list.col[mk];			// a[m][q] 为非零元素 
            if(q!=m) nzc[q]--;
        }
        mk=list.right[mk];
    }
    // 搜索第 m 列非零元素 
    mk=list.cp[m];					 // mk 为第 m 列元素号 
    while(mk!=0){					 // 搜索第 m 列非零元素 
        if(list.del[mk]==0) {		 // 未删元素 
            q=list.row[mk];			 // a[q][m] 为非零元素 
            if(q!=m) nzr[q]--;
        }
        mk=list.down[mk];
    }
};
//==========================================================================
void CSparseMatrix::del_r_f(INT2 m, LINKED_LIST &list, INT2 *fill)
/* [ 功  能 ]   在双链接表 list 所描述的稀疏矩阵中, 删去第 m 主行, 对其它候
 *              选主元产生的填元数 fill[] 的影响
 * =======================================================================*/
{
    INT2 mk,q,qk,j;

    mk=list.rp[m];						// mk 为第 m 行元素号 
    while(mk!=0){						// 搜索第 m 行非零元素 
        if(list.del[mk]==0) {			// 未删元素 
            q=list.col[mk];				// a[m][q] 为非零元素 
            qk=list.rp[q];				// qk 为第 q 行元素号 
            while(qk!=0){				// 搜索第 q 行非零元素 
                if(list.del[qk]==0) {   // 未删元素 
                    j=list.col[qk];		// a[q][j] 为非零元素 
                    /* 判断 a[m][j] 是否是零。若是零, 则候选主元 q 通过 
                       a[m][q] 和 a[q][j] 产生的填元 a[m][j] 将由于 m 行
                       的删除而不存在, fill[q] - 1。*/
                    if(ele_num(m,j,list)==0) fill[q]--; 
                }
                qk=list.right[qk];
            } 
        }
        mk=list.right[mk];
    }
};
//=========================================================================