sparsem.cpp

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

					}
				}
            }
        }
    }
    return(0);
};
//==========================================================================
char CSparseMatrix::num_lu(TRIANGULAR_TABLE_C &table)
// [ 功  能 ]   根据复数稀疏矩阵的三角形表 table 进行数值 LU 分解。 
{
    char luid;
    INT2 n;						// 矩阵阶数  
    INT2 k,kk,kkk,r,rr,c,cc,t;
    INT2 lu;					// LU 分解修正元素号 
    double xre,xim,vre,vim; 

    n=table.n;
    lu=0;
    for(k=1;k<n;k++){			// 第  k 步分解, 主元为 k 号元素 
        xre=table.vre[k];
        xim=table.vim[k];
        luid=div_c(1.0,0,xre,xim,&xre,&xim);
        if(luid!=0) return(1);
        c=table.urp[k];			// 主行指针 
        cc=table.lcp[k];		// k-U 行非零元素号从 c 到 cc-1 
        for(kk=c;kk<cc;kk++){   // 归一化 U 主行 
            vre=table.vre[kk];
            vim=table.vim[kk];
            mul_c(vre,vim,xre,xim,&vre,&vim);   
            table.vre[kk]=vre;
            table.vim[kk]=vim;
        }
        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];		// 修正元素号 
                mul_c(table.vre[kk],table.vim[kk],table.vre[kkk],
                    table.vim[kkk],&vre,&vim);
                table.vre[t]=table.vre[t]-vre;
                table.vim[t]=table.vim[t]-vim;
            }
        }
    }
    return(0);
}
// ==========================================================================
void CSparseMatrix::vector_fillin(TRIANGULAR_TABLE &table)
// 根据已输入右端向量的稀疏矩阵的三角形表 table 确定和插入向量填元。
{
    char f;
    INT2 n;					// 矩阵阶数  
    INT2 b;					// 右端向量元素号 
    INT2 nza;				// 非零元素个数 
    INT2 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 CSparseMatrix::vector_fillin(TRIANGULAR_TABLE_C &table)
// 根据已输入右端向量的稀疏矩阵的三角形表 table 确定和插入向量填元。
{
    char f;
    INT2 n;			// 矩阵阶数  
    INT2 b;			// 右端向量元素号 
    INT2 nza;		// 非零元素个数 
    INT2 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 CSparseMatrix::sym_fe(TRIANGULAR_TABLE &table)
/* 根据已输入右端向量的稀疏矩阵的三角形表 table 进行符号前消,
   确定 table.fep[]。*/
{
    INT2 n;			// 矩阵阶数  
    INT2 b;			// 右端向量元素号 
    INT2 nza;		// 非零元素个数 
    INT2 fe;
    INT2 k,kk,l,ll,r,i;

    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;				// b 号的下一个右端向量非零元素号 
        for(kk=l;kk<ll;kk++){   // 扫描第 k-L 列, kk 为其非零元素号 
            i=table.roco[kk];   // kk 号右端向量非零元素行号
            while(r<nza+1){		// 搜索 b 以下的右端向量非零元素 
                if(table.roco[r]==i){
                    fe++;
                    table.fep[fe]=r;
                    r++;
                    break;
                }
                r++;
            }
        }
        b++;
    }
    table.fen=fe;
};
// ==========================================================================
void CSparseMatrix::sym_fe(TRIANGULAR_TABLE_C &table)
/* 根据已输入右端向量的稀疏矩阵的三角形表 table 进行符号前消,
   确定 table.fep[]。*/
{
    INT2 n;			// 矩阵阶数 
    INT2 b;			// 右端向量元素号 
    INT2 nza;		// 非零元素个数 
    INT2 fe;
    INT2 k,kk,l,ll,r,i;

    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;					// b 号的下一个右端向量非零元素号 
        for(kk=l;kk<ll;kk++){   // 扫描第 k-L 列, kk 为其非零元素号 
            i=table.roco[kk];   // kk 号右端向量非零元素行号 
            while(r<nza+1){		// 搜索 b 以下的右端向量非零元素 
                if(table.roco[r]==i){
                    fe++;
                    table.fep[fe]=r;
                    r++;
                    break;
                }
                r++;
            }
        }
        b++;
    }
    table.fen=fe;
};
// =========================================================================
char CSparseMatrix::num_fe(TRIANGULAR_TABLE &table)
{
    INT2 n;			// 矩阵阶数 
    INT2 b;			// 右端向量元素号 
    INT2 nza;		// 非零元素个数 
    INT2 fe;
    INT2 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);
}

//==========================================================================
BOOL CSparseMatrix::JudgeZeroT(TRIANGULAR_TABLE &table)                                       
// 判断矩阵对角元素是否为零
{
    INT2 n;			// 矩阵阶数 
    INT2 b;			// 右端向量元素号 
    INT2 nza;		// 非零元素个数 
    INT2 fe;
    INT2 k;
    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 FALSE;
        b++;
    }

	return TRUE;
}
//==========================================================================
char CSparseMatrix::num_fe(TRIANGULAR_TABLE_C &table)
/* [ 功  能 ]   根据复数稀疏矩阵的三角形表 table 进行数值前消。右端向量化为 
                中间解向量。 */
{
    char feid;
    INT2 n;				// 矩阵阶数  
    INT2 b;				// 右端向量元素号 
    INT2 nza;			// 非零元素个数 
    INT2 fe;
    INT2 k,kk,kkk,l,ll;
    double xre,xim,vre,vim;

    n=table.n;
    nza=table.nza;
    b=table.lcp[n];		// 右端向量从 lup[n] 开始存放 
    fe=0;				// fep[] 标号 

    while(b<nza){								// 处理第 b 号右端向量非零元素 
        k=table.roco[b];						// 向量行号 
        xre=table.vre[k];						// 对应主元实部 
        xim=table.vim[k];						// 对应主元虚部 
        feid=div_c(1.0,0,xre,xim,&xre,&xim);    // x=1.0/x 
        if(feid!=0) return(1);
        mul_c(table.vre[b],table.vim[b],xre,xim,&xre,&xim); // 归一化
        table.vre[b]=xre;
        table.vim[b]=xim;
        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];
            mul_c(table.vre[kk],table.vim[kk],xre,xim,&vre,&vim); 
            table.vre[kkk]=table.vre[kkk]-vre;
            table.vim[kkk]=table.vim[kkk]-vim;
        }
        b++;
    }
    // 归一化最后一个向量元素 
    k=table.roco[nza];
    xre=table.vre[k];
    xim=table.vim[k];
    feid=div_c(1.0,0,xre,xim,&xre,&xim);
    if(feid!=0) return(1);
    mul_c(table.vre[nza],table.vim[nza],xre,xim,&vre,&vim);
    table.vre[nza]=vre;
    table.vim[nza]=vim;
    return(0);
}
//==========================================================================
void CSparseMatrix::num_bs(TRIANGULAR_TABLE &table,double *x)
{
    INT2 n;					// 矩阵阶数  
    INT2 nza,nn;
    INT2 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;
    }
};
//==========================================================================
void CSparseMatrix::num_bs(TRIANGULAR_TABLE_C &table,double *xre,double *xim)
// [ 功  能 ]   根据数值前消后