mini_degree.cpp

最小度算法,为矩阵做最小度排序. 在有限元求解中较为常用.

//************************************************************
//  mini_degree.cpp
//  基于一维半带宽存贮的最小度算法
//************************************************************

#include<iostream>
#include <fstream>
#include <string>
#include <cmath>
using namespace std;

//自由度总数宏定义
#define degree 8
//const int degree = 4;
//一维总刚中元素总数
#define TotalGK 25

//定义指针类型
typedef int* IntArrayPtr;

//计算未消去前各顶点的相邻集与度
void Initi_cla(int n_degree,  int **NE, double *pGK, int *pDiag, int *nBandWidth, int *nAdj);

//生成重排序数组
void Reorder_ID(int n_degree, int **NE, int *nAdj, int *ID, int iNum);

//更新相关顶点的度
void Update_Degree(int n_degree, int **NE, int *nAdj, int iNum);



int main()
{	
	//
	cout<<"总自由度数:degree= "<<degree<<endl;
	//

	//顶点拓扑,一个顶点的相邻顶点集,不包括本身
	//相邻顶点可以定义小些而不用degree,待改
	
	//动态数组定义NE
	IntArrayPtr *NE= new IntArrayPtr[degree];
	for(int i=0; i<degree; i++)
		NE[i]= new int[degree];

	//主对角元定义
	int *pDiag= new int [degree];
	//主对角地址值
	int values_pDiag[degree]={0,2,5,9,13,17,20,24};
	pDiag=values_pDiag;

	cout<<endl<<"输入的pDiag="<<endl;
	for( i=0; i<degree; i++)
	cout<<"pDiag["<<i<<"]= "<<pDiag[i]<<endl;
	
	//定义总刚
	double *pGK=new double[TotalGK];
	//一维刚度矩阵赋值
	double values_pGK[TotalGK]={1,1,8,4,0,10,2,12,0,1,1,1,8,4,5,0,0,12,1,1,1,1,0,0,1};
    pGK=values_pGK;
	//
	cout<<endl<<"输入的pGK="<<endl;
	for(i=0; i<6; i++)
	cout<<"pGK["<<i<<"]= "<<pGK[i]<<endl;
	//

	//半带宽(行左边第一个非零元到对角元,含对角元在内)
	//半带为1时,直接做为主元??
	int *nBandWidth = new int[degree];
	//顶点度数
	int *nAdj = new int[degree];
	//最小度顶点所在顶点号
	int iNum=0;
	//重排序数组
	int *ID = new int[degree];

	//求顶点的相邻点及未消元前的顶点度计算
	Initi_cla(degree, NE, pGK, pDiag, nBandWidth, nAdj);

	//生成重排序数组ID
	Reorder_ID(degree, NE, nAdj, ID, iNum);
	
	//最终生成的ID数组
	for(int i0=0; i0<degree; i0++)
		cout<<"ID["<<i0<<"]="<<ID[i0]<<endl;

	return 0;
}





//*****************************************************************************************
//顶点的相邻集合初始化,求得的NE[i][]中前nAdj[i]个为相邻顶点
void Initi_cla(int n_degree, int **NE, double *pGK, int *pDiag, int *nBandWidth, int *nAdj)
{
	//相邻顶点拓扑集合初始化,均赋-1,i顶点的相邻集合中不包括i顶点自身
	//ii为自由度序号
	cout<<endl<<"NE的初始化:"<<endl;
	for(int ii=0; ii<n_degree; ii++)
	{
		//jj为ii的相邻顶点,一般不要设太大
		for(int jj=0; jj<n_degree; jj++)
		{
			NE[ii][jj]=-1;
			cout<<"NE["<<ii<<"]"<<"["<<jj<<"]="<<NE[ii][jj]<<endl;
		}
	}

	//各顶点的相邻顶点数(度数),初始化,赋0
	cout<<endl<<"nAdj的初始化: "<<endl;
	for(int t=0; t<n_degree; t++)
	{
		nAdj[t]=0;
		cout<<"nAdj["<<t<<"]="<<nAdj[t]<<endl;
	}

	//对半带宽初始化,初始化为1
	cout<<endl<<"nBandWidth的初始化: "<<endl;
	for(t=0; t<n_degree; t++)
	{
		nBandWidth[t]=1;
		cout<<"nBandWidth["<<t<<"]="<<nBandWidth[t]<<endl;
	}

	//半带宽计算,从第1行开始. 第零行为半带宽为1
	for (int i=1; i<n_degree; i++)
	{
		//半带宽计算,0行的为1,从第1行起计算,包含对角元
		nBandWidth[i]=pDiag[i]-pDiag[i-1];
	}
	
	//对各顶点/行循环,依次求NE,nAdj
	for (i=1; i<n_degree; i++)
	{
		//判断第0个顶点的相邻顶点集合
		//判断第0列各行(1~n-1)元素是否为0
		if(nBandWidth[i]==i+1)
		{
			NE[0][nAdj[0]]=i;
			cout<<"第0行NE[0]["<<i<<"]="<<NE[0][nAdj[0]]<<endl;
			//0顶点的相邻顶点数累加
			nAdj[0]++;
		}
		
		//求NE[i][j],i=1~n-1
		//对行半带宽内素循环,判断是否为0,求非零元素(度)数
		//如果带宽为1,则无需求NE与nAdj,为初始值
		if(nBandWidth[i]>1)
		{
			//对带宽内元素循环
			for(int j=1; j<nBandWidth[i]; j++)
			{
				//判断行半带宽内的元素是否为0
				if(pGK[pDiag[i]-j] !=0)
				{
					//第i个顶点的拓扑顶点
					NE[i][nAdj[i]]=i-j;
					cout<<"NE["<<i<<"]["<<nAdj[i]<<"]="<<NE[i][nAdj[i]]<<endl;
		
					//相邻顶点数累计
					nAdj[i]++;
					cout<<"nAdj["<<i<<"]="<<nAdj[i]<<endl;
				}
			}
		}
				

				//判断i列的i+1~n-1行中不为零的元素
				//当前行不是最后一行时
				if(i != n_degree-1)
				{
					for(int kk=i+1; kk<n_degree; kk++)
						if(nBandWidth[kk]>=(kk-i+1) && pGK[pDiag[kk]-(kk-i)] !=0)
						{
							NE[i][nAdj[i]]=kk;
							nAdj[i]++;
						}
				}	
	}
	
			//一些结果输出
			cout<<endl;
			cout<<"nBandWidth及Nadj[]的值="<<endl;		
			for(int i1=0; i1<n_degree; i1++)
			{
				cout<<"nBandWidth[ "<<i1<<"]="<<nBandWidth[i1]<<endl;
				cout<<"nAdj["<<i1<<"]="<<nAdj[i1]<<endl;
			}

			cout<<endl;
			cout<<"NE的未消去前值="<<endl;
			for(int i2=0; i2<n_degree; i2++)
			for(int j1=0; j1<n_degree; j1++)
			{
				cout<<"NE["<<i2<<"]"<<"["<<j1<<"]="<<NE[i2][j1]<<endl;
			}
}




//*******************************************************************************************
//Reorder_ID 函数,对NE及nAdj更新,生成重排序数组ID
void Reorder_ID(int n_degree, int **NE, int *nAdj, int *ID, int iNum)
{
	//ID数组初始化为0
	for(int i=0; i<n_degree; i++)
	{
		ID[i]=0;
	}

	//求ID[]重排序数组
	for(int ii=0; ii<n_degree; ii++)
	{
		int temp=0;
		temp=nAdj[0];
		iNum=0;
		
		//初始度[2 2 3 1]
		//
		for(int i11=0; i11<n_degree; i11++)
		{
			cout<<"nAdjjjjjjjjjjjjjjjjj= "<<nAdj[i11]<<endl;
		}
		//

		//对nAdj[],选出其中最小的一个
		for(int i=1; i<n_degree; i++)
		{
			if(temp>nAdj[i])
			{
				temp=nAdj[i];
				//最小的度所在位置
				iNum=i;
			}
		}
		//
		
		cout<<"第"<<ii<<"次的最小度所在顶点号= "<<iNum<<endl;
		
		//将已求的最小度所在序号赋给排序数组ID
		ID[ii]=iNum;
		//cout<<"ID["<<ii<<"]="<<ID[ii]<<endl;
		
		//更新与iNum顶点相邻的顶点的度nAdj
		Update_Degree(degree, NE, nAdj, iNum);
	}
}




//********************************************************************************************
//更新相关顶点的度
void Update_Degree(int n_degree, int **NE, int *nAdj, int iNum)
{
	//将iNum(最小度)顶点的相邻顶点数目赋给temp1
	int temp1=0, temp2=0;
	temp1=nAdj[iNum];

	//最小度顶点的拓扑关系全为负1
	for(int ij=0; ij<degree; ij++)
		NE[iNum][ij]=-1;

	//aNum[],求与最小度顶点相邻的顶点的集合
	int *aNum = new int [degree];
	//最小度相邻顶点号初始化
	for(int aa=0; aa<n_degree; aa++)
		aNum[aa]=0;

    //更新与最小度顶点相邻顶点的度,先将与最小度顶点相邻顶点的度减1
	for(int i=0; i<temp1; i++)
	{
		//求与最小度顶点相邻的顶点号
		aNum[i]=NE[iNum][i];
		//消去最小度顶点,将与最小度顶点相邻的顶点的度减去1
		nAdj[aNum[i]]=nAdj[aNum[i]]-1;
	}

	//对最小度顶点的所有相邻顶点循环,以更新它们的度及拓扑关系
	for(int kk=0; kk<temp1; kk++)
	{
		/////////////////////
		//对已消去的最小度顶点,在其相邻顶点拓扑集中顶点号变为负1
				for(int nn=0; nn<nAdj[aNum[kk]]; nn++)
				{
					if(NE[aNum[kk]][nn]==iNum) 
						NE[aNum[kk]][nn]=-1;
				}

				//////  
				//对最小度顶点的相邻顶点,将NE中顶点号重排,将-1放置最后
				for(int ss=0; ss<nAdj[aNum[kk]]+1; ss++)
				if(NE[aNum[kk]][ss]==-1)
				{
					NE[aNum[kk]][ss]=NE[aNum[kk]][nAdj[aNum[kk]]];
					NE[aNum[kk]][nAdj[aNum[kk]]]=-1;
				}
				/////////

		//对该最小度相邻顶点集合中各顶点的相邻顶点循环,不包括消去的最小度顶点
		for(int ll=0; ll<nAdj[aNum[kk]]; ll++)
		{			
			//如果最小度相邻顶点集合中多于一个顶点
			if(nAdj[iNum]>1)
			{
				//对其它相邻顶点
				for(int mm=ll+1; mm<nAdj[iNum]; mm++)
				{
					//最小度顶点的两个相邻顶点不相邻 && 
					if(NE[aNum[kk]][ll]!=aNum[mm] && kk != mm)
					{
						//顶点拓扑更新
						NE[aNum[kk]][nAdj[aNum[kk]+1]]=aNum[mm];
						//????
						NE[aNum[mm]][nAdj[aNum[mm]+1]]=aNum[kk];
						//相关度的更新
						nAdj[aNum[kk]]++;
						nAdj[aNum[mm]]++;
					}
				}
			}
		}
				
	}
		
	//将已消去的iNum顶点的度nAdj赋大值
	nAdj[iNum]=10000;
}