zp182931zyy 发表于 2011-4-4 21:33:40

问题求助

在有限元编程中,涉及到带宽的问题,本人的数据由ANSYS生成,但是结点编号没有规律,如果不进行带宽优化,则求解较慢,甚至可能得不到结果,不知道哪位编过带宽优化这样的程序,麻烦给点意见,谢谢!

eccentric 发表于 2011-4-5 12:56:15

CM算法或RCM算法。

eccentric 发表于 2011-4-5 12:58:52

#pragma once
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
//图的邻接表类及相关计算:
struct chain
{
        int data;
        chain* next;
};

//RCM算法类:
class RCM
{
public:
        RCM(const int Node);                                                        //构造函数:接受单元信息矩阵。
        int* renumRCM(int* Elm,int* ElmIndex,int ElmNum,int* oldlist);       
                                                                                                        //对点列oldlist RCM编号,返回编号数组。
        ~RCM(){delete[] adj_list;delete[] adj_point;};
        void outputtest(int pointA);                                        //输出检测:输出PointA的邻接点。

private:
        int degree(int pointA);                                                        //计算度。
        void Addin(int pointA,int pointB);                                //将点对(a,b)加入邻接表。
        int* pseudodia(int &pointA,int &dia,int &num);        //查询由pointA出发的伪直径,返回pointA的层级表。
        void state(int* Elm,int* ElmIndex,int ElmNum);        //统计当前信息
private:
        chain* adj_list;               
        chain* adj_point;                //邻接表。
        int Nodenum;                        //节点数。
        int* flag;                                //节点标记。

};
int RCM::degree(int pointA)
{        chain *p=adj_point.next;
        int deg=0;
        while(p!=adj_point.next){
                deg++;
                p=p->next;
        }
        return deg;
}
int* RCM::pseudodia(int &pointA,int &dia,int &num)
{        //从任意点开始建立层级结构:(广度优先遍历)
        int q;                        //标记矩阵:标记访问次序
        int front=0;
        int rear=1;
        dia=0;
        for(int s=0;s<Nodenum;s++)
                flag=0;
        flag=1;
        q=pointA;
        q=0;rear++;
        while(front!=rear){
                if(q==0){
                        dia++;
                        q=0;rear++;
                        front++;       
                        if(q==0) break;
                }
                int p1=q;front++;
                chain *v=adj_point.next;
                while(v!=adj_point.next){
                        if(flag==0){
                                flag=1;
                                q=v->data;rear++;
                        }
                        v=v->next;
                }
        }
        //统计最后一层的最小度:

        dia=dia-1;
        int min_degree=degree(q);
        rear-=2;
        int i=rear;
        int tempA=q;
        while(q!=0){
                if(min_degree>degree(q)){
                        min_degree=degree(q);
                        tempA=q;
                }
                i--;
        }
        pointA=tempA;
        int *rootPA=new int;
        num=rear;
        for(i=0;i<rear;i++)
                rootPA=q;
        return rootPA;
}
void RCM::outputtest(int pointA)
{        chain *p=adj_point.next;
        cout<<"Output adjacent point of "<<pointA<<":"<<endl;
        while(p!=adj_point.next){
                cout<<p->data<<" "<<ends;
                p=p->next;
        }
        cout<<endl;
}
RCM::RCM(const int Node)
{        //生成一个Node长度的空链表;Adj_point指向相应元素:
        int i=0;
        adj_list=new chain;
        adj_list->data=0;
        adj_list->next=NULL;
        adj_point=new chain;
        (adj_point).data=1;
        (adj_point).next=adj_list;
        chain* point=adj_list;
        for(i=1;i<Node;i++){
                chain* temp=new chain;
                point->next=temp;
                temp->data=0;
                temp->next=NULL;
                point=temp;
                adj_point.data=i+1;
                adj_point.next=temp;
        }

        adj_point.next=NULL;               
        //初始化单元数&标记向量。
        Nodenum=Node;
        flag=new int();
        for(i=0;i<Node;i++)
                flag=0;
}
void RCM::Addin(int pointA,int pointB)
{        pointA-=1;
        pointB-=1;
        int flag=0;
        chain *p=adj_point.next;
        chain *p1=NULL;
        while(p!=adj_point.next){
                if(p->data==pointB+1){
                        flag=1;
                        break;
                }
                p1=p;
                p=p->next;
        }
        if(flag==0){
        //p1和adj_point间插入pointB:
                chain* newpoint=new chain;
                newpoint->data=pointB+1;
                newpoint->next=adj_point.next;
                p1->next=newpoint;
                //相应在pointB的邻点集中加入pointA:
                chain *p2=adj_point.next;
                chain *p3=NULL;
                while(p2!=adj_point.next){
                        p3=p2;
                        p2=p2->next;
                }
                chain* newpoint1=new chain;
                newpoint1->data=pointA+1;
                newpoint1->next=adj_point.next;
                p3->next=newpoint1;                       
        }

}

int* RCM::renumRCM(int* Elm,int* ElmIndex,int ElmNum,int* oldlist)
{        //遍历ElmIndex,生成邻接表:
        int i=0;
        cout<<"strat RCM"<<endl;
        cout<<"ElmNum: "<<ElmNum<<endl;       
        cout<<"NodeNum: "<<Nodenum<<endl;
        for(i=0;i<ElmNum;i++){
                int e=ElmIndex;
                int en=ElmIndex;
                for(int im=e;im<en-1;im++)
                        for(int in=im+1;in<en;in++)
                                Addin(Elm,Elm);               
        }
        adj_list=adj_list->next;
        chain *pp=adj_list;
        while(pp->next!=NULL){
                if(pp->next->data==0){
                        pp->next=pp->next->next;
                        pp=pp->next->next;
                }
                else pp=pp->next;
        }
        for(i=0;i<Nodenum;i++)
                adj_point.next=adj_point.next->next;
        //统计矩阵信息:(计算当前矩阵的列高和)
        state(Elm,ElmIndex,ElmNum);
        //对每个节点的邻点集按度排序:(直接排序法)
        for(i=0;i<Nodenum;i++){
                int deg=degree(adj_point.data);
                int* a=new int;
                chain* p=adj_point.next;
                for(int j=0;j<deg;j++){
                        if(p!=adj_point.next){
                                a=p->data;
                                p=p->next;
                        }
                }
                int temp;
                for(int k=0;k<deg-1;k++){
                        int pmin=k;
                        for(int n=k+1;n<deg;n++)
                                if(degree(a)>degree(a)){
                                        pmin=n;
                                }
                        if(pmin!=k){
                                temp=a;a=a;a=temp;
                        }
                }
                chain *pp=adj_point.next;
                int m=0;
                while(pp!=adj_point.next){
                        (pp->data)=a;
                        pp=pp->next;
                }
        }

        //从oldlist中选择起始点PointA,并生成连根层级结构。
        int pointA=5;
        int *rootPA=NULL;
        int numPA;
        int dia=0;
        int dia1=1;
        while(dia1>dia){
                dia=dia1;
                rootPA=pseudodia(pointA,dia1,numPA);
        /*        cout<<"PointA:"<<pointA<<endl;
                cout<<"dia1    :"<<dia1<<endl;
                cout<<"dia   :"<<dia<<endl;
                cout<<"numPA   :"<<(numPA-dia1)<<endl;
        */
        }
        //遍历oldlist,对每一点,根据邻点集编号。
        int* rcmIndex=new int;
        int j=Nodenum;
        for(i=0;i<numPA;i++){
                if(rootPA!=0)                       
                        rcmIndex-1]=j--;
        }
        //统计编号后信息:
        for(i=0;i<ElmIndex;i++){
                Elm=rcmIndex-1];
        }
        cout<<endl;
        cout<<endl;
        cout<<"After RCM re-numbering:"<<endl;
        state(Elm,ElmIndex,ElmNum);
        return rcmIndex;
}
void RCM::state(int* Elm,int* ElmIndex,int ElmNum)
{        int* cm=new int;
        int i;
        for(i=0;i<Nodenum;i++) cm=1;
        for(i=0;i<ElmNum;i++){
                int temp;
                int k=0;
                for(k=0;k<20;k++) temp=0;
                int        e=ElmIndex;
                int        en=ElmIndex;
                k=0;
                for(int j=e;j<en;j++) temp=Elm;
                int min_temp=temp;
                for(k=0;k<20;k++)
                        if((min_temp>temp)&&(temp>0))        min_temp=temp;
                for(k=0;k<20;k++)
                        if(temp>0)
                                cm-1]=((cm-1]<(temp-min_temp+1))?(temp-min_temp+1):cm-1]);
       
        }
        for(i=1;i<Nodenum;i++) cm+=cm;
        int EffectElm=2*cm;
        //for(i=0;i<Nodenum;i++) cout<<cm<<endl;
        cout<<"Number of Element: "<<Nodenum<<"."<<endl;
        cout<<"Element in effect: "<<EffectElm<<"("<<setprecision(3)<<(float)100*EffectElm/Nodenum/Nodenum<<"%)."<<endl;

}


无责任随便贴一个以前写的。刚看了一下发现自己都看不懂。

zp182931zyy 发表于 2011-4-7 18:45:18

这个难道没人会吗?
页: [1]
查看完整版本: 问题求助