Progetto MECOM, Programma Operativo Plurifondo 94/99
Misura 4.4 ”Ricerca scientifica e tecnologica”
Sviluppi ed applicazioni della meccanica computazionale
nella progettazione strutturale in campo civile ed industriale
.
Una implementazione in C++
di un algoritmo di rinumerazione per
la riduzione dell’ampiezza di banda
di matrici sparse
Antonio Bilotta
Report n. 193
Dicembre 1998
Dipartimento di Strutture
Università della Calabria 87030 Rende, Cosenza, Italy
tel.: +39 0984 494032 - fax: +39 0984 494045
e.mail: [email protected]
1
Una implementazione in C++
di un algoritmo di rinumerazione per
la riduzione dell’ampiezza di banda
di matrici sparse
A. Bilotta
Dipartimento Strutture
Università della Calabria, 87030 Rende (CS) – Italy
e.mail abilotta @ labmec.unical.it
Introduzione
Molti problemi ingegneristici richiedono la soluzione di un sistema lineare di equazioni
del tipo Ax = b, dove x è il vettore che raccoglie i parametri incogniti, b è il vettore
dei termini noti, A è una matrice N xN simmetrica definita positiva e, tipicamente,
con una struttura a banda. I due metodi più largamente adottati nei programmi
ad elementi finiti per la soluzione di tale sistema sono gli schemi profile o skyline e
frontale.
Lo schema profile si basa, nelle implementazioni solitamente adottate, su un
assemblaggio completo della matrice di rigidezza e quindi la sua efficacia computazionale dipende da due aspetti: (i) ottimizzazione delle risorse richieste per
la memorizzazione della matrice; (ii) riduzione dei tempi di calcolo necessari alla
soluzione del sistema. Lo schema frontale invece scorre le righe del sistema lineare
utilizzando al volo solo le equazioni effettivamente attive durante ogni singolo passo
del processo di decomposizione della matrice ed eliminazione delle incognite, perciò
la sua efficienza dipende essenzialmente dai tempi di calcolo richesti dalla soluzione
del sistema.
2
– Algoritmo di rinumerazione: implementazione in C++ –
3
Per quel che riguarda il primo aspetto, ovvero la memoria accupata dalla matrice di rigidezza, una sua ottimizzazione dipende dall’ordinamento delle equazioni
il quale, nel metodo ad elementi finiti, corrisponde alla numerazione adottata per
i nodi del modello discreto. Di conseguenza tale numerazione, vista la struttura
a banda della matrice di rigidezza, deve essere effettuata in modo da minimizzare
la somma delle ampiezze di banda delle singole righe. In maniera più formale è
possibile definire la quantità P denonimata profile in base alla seguente formula:
P =
N
bi
(1)
i=1
bi rappresenta l’ampiezza di banda della singola riga del sistema ovvero la differenza
di indici dell’ultimo coefficiente non nullo Aij , j > i, e del termine sulla diagonale
principale Aii . Risulta quindi fondamentale minizzare il valore assunto da P . D’altra parte anche i tempi di calcolo richiesti dalla soluzione del sistema, ad esempio
mediante una classica elinazione alla Gauss, dipendono dalle ampiezze di banda bi .
Infatti introducendo la quantità
12
N
1 2
bi
(2)
F =
N
i=1
(root–mean–square wavefront) è possibile mostrare come il numero di addizioni e
moltiplicazioni richieste per il processo di decomposizione della matrice sia di ordine
O(N F 2 ) [1]. Quindi il problema di una soluzione efficiente, rispetto alle risorse di
memoria richieste e ai tempi di calcolo impiegati, del sistema lineare Ax = b può
essere ridotto al problema della minimizzazione di P . Tale obbiettivo, per un modello
discreto ad elementi finiti, è raggiungibile mediante una numerazione adeguata dei
nodi che definiscono la struttura analizzata.
Nel seguito si intende descrivere l’implementazione in C++, noto linguaggio di
programmazione orientato agli oggetti, di un algoritmo di rinumerazione dei nodi
proposto da S. W. Sloan in [1]. Vengono quindi ripercorse rapidamente le caratteristiche salienti dell’algoritmo e, mediante dei risultati numerici, si evidenzia l’effettiva
convenienza di tale algoritmo rispetto al più noto procedimento di rinumerazione di
Cuthill e McKee [2]. Infine si riporta il listato del codice implementato il quale può
essere facilmente incluso in un tipico programma per l’analisi ad elementi finiti.
Notazione
Il problema della numerazione dei nodi di una discretizzazione FEM è equivalente
a quello della numerazione di un grafo non orientato. Conviene quindi rifersi alla
– Algoritmo di rinumerazione: implementazione in C++ –
4
teoria dei grafi ed introdurre un minimo di terminologia utile alla definizione del
problema.
Un grafo è formato da un insieme di membri detti nodi e da un insieme di coppie
di nodi distinti detti lati. Quando lungo i lati non viene stabilito un particolare
senso di percorrenza allora il grafo si dice non orientato. Si noti che loops, ossia
nodi connessi con se stessi, e lati multipli, ossia coppie di nodi connesse da più lati,
non sono previsti per le tipologie di grafi prese in esame. Per grado di un nodo si
intende il numero dei lati che convergono su quest’ultimo e può assumere un valore
intero compreso tra 1 e N . Due nodi si dicono connessi se è possibile stabilire una
sequenza di lati, percorso, che li congiunge. Se ciò è possibile per qualsiasi coppia di
nodi allora il grafo si dice connesso.
La distanza tra una coppia di nodi è definita dal numero di lati appartenenti al
percorso più breve che li congiunge, mentre il diametro di un grafo è la massima delle
distanze tra ogni coppia di nodi. Tutte le possibili coppie di nodi che definiscono
il diametro di un grafo sono denominati nodi periferici. Risulta utile introdurre,
per il tipo di algoritmo che verrà descritto, lo pseudo-diametro, ovvero una stima
del diametro del grafo, e, di conseguenza, i nodi pseudo-periferici, ovvero nodi che
definiscono uno pseudo-diametro.
Infine un concetto che risulta essere centrale per molti algoritmi di ordinamento: la struttura a livelli a partire da un nodo assegnato o radice. Si tratta di un
partizionamento dei nodi che assegna ciascun nodo ad uno dei livelli l1 , l2 , . . . , lh(r)
in base alla propria distanza dalla radice r. Tutti i nodi che si trovano alla stessa
distanza dalla radice appartengono allo stesso livello. La profondità di una struttura
a livelli, h(r), è data semplicemente dal numero di livelli che la compongono, mentre
la larghezza è fornita dal numero massimo di nodi che appartengono ad un singolo
livello. In modo più formale tale partizionamento può essere rappresentato come
L(r) = {l1 , l2 , . . . , lh(r) } dove
• l1 = r.
• Per i > 1 li è costituito dall’insieme dei nodi che si trovano ad una distanza
pari a i − 1 dalla radice.
L’algoritmo di Sloan
La strategia dell’algoritmo di rinumerazione di Sloan si basa su due passi ben distinti:
(i) selezione di due nodi periferici e della struttura a livelli ad essi associata la
quale, avendo come radice uno dei due nodi periferici, sarà adeguatamente stretta
– Algoritmo di rinumerazione: implementazione in C++ –
5
ed allungata; (ii) rinumerazione dei nodi utilizzando la struttura a livelli ottenuta
nel passo precedente.
Selezione dei nodi periferici
La numerazione vera e propria dei nodi viene preceduta da una fase in cui si effettua
una stima del diametro del grafo e quindi si procede alla selezione di due nodi
pseudo-periferici da utilizzare successivamente. Tale compito può essere portato a
termini nel modo seguente:
1 Ricerca, fra tutti i nodi del grafo, del nodo s avente grado minore.
2 Assumendo come radice s generazione della struttura a livelli L(s).
3 Ordinamento per grado ascendente dei nodi contenuti nel’ultimo livello lh(s) .
4 A partire da lh(s) , generazione di una lista di nodi q contenente per ciascun
grado solo un nodo.
5 Inizializzazione w(e) = ∞.
6 Per ogni nodo i ∈ q costruzione della struttura L(i) = {l1 , l2 , . . . , lh(i) }. Se
h(i) > h(s) e w(i) < w(e) allora s = i e ritorno al passo 3. Altrimenti se
w(i) < w(e) si pone e = i e w(e) = w(i).
7 Uscita con la determinazione dei nodi iniziale s e finale e.
Numerazione dei nodi
Dopo avere determinato una coppia di nodi pseudo-periferici si può procedere alla
rinumerazione dei nodi. La tecnica di rinumerazione fa uso delle seguenti definizioni:
un nodo che è stato rinumerato è detto postattivo; nodi adiacenti a nodi postattivi
sono invece attivi; qualsiasi nodo adiacente a nodi nodi attivi ma che non è postattivo
e neanche attivo viene definito preattivo; tutti gli altri nodi sono inattivi. Premesso
ciò l’algoritmo si articola nei seguenti passi:
1 Partenza dell’algoritmo con i nodi terminali s e e.
2 Costruzione della struttura a livelli prendendo come radice il nodo finale e e
quindi calcolo per ogni nodo i della distanza, δi , da e.
3 Per tutti i nodi viene assegnato lo stato inattivo e viene calcolata la priorità
iniziale pi mediante la formula:
pi = W1 δi − W 2(di + 1)
W1 e W2 sono dei valori interi e rappresentano dei pesi, di è il grado del nodo
i.
– Algoritmo di rinumerazione: implementazione in C++ –
6
4 Inizializzazione a zero del contatore l dei nodi che sono stati rinumerati e
della coda di priorità q. Assegnazione al nodo s dello stato di preattivo e suo
inserimento nella coda ponendo n = 1 e qn = s.
5 Fintantoché la coda di priorità è piena, n > 0, ripetizione dei passi dal 6 al 9.
6 Ricerca nella coda q del nodo avente la massima priorità: pi = maxj∈q {pj }.
Sia m l’indice del nodo i tale che qm = i.
7 Eliminazione del nodo i dall coda q ponendo qm = qn e decremento del contatore n = n − 1. Se il nodo non è preattivo si salta al passo successivo.
Altrimenti si esaminano tutti in nodi j adiacenti ad i incrementandone la priorità in base alla formula pj = pj + W2 . Inoltre se il nodo j è inattivo viene
inserito nella coda di priorità, n = n+1 e qn = j, ponendolo in stato preattivo.
8 Rinumerazione del nodo i. Ciò comporta l’incremento del contatore l, l = l+1,
e l’assegnazione ηi = l, dove η contiene le nuove etichette per i nodi. Passaggio
allo stato postattivo per il nodo i.
9 Esame di tutti i nodi j adiacenti al nodo appena rinumerato i. Se il generico
nodo j non è preattivo non si effettua alcuna operazione. Altrimenti si assegna
al nodo j lo stato di nodo attivo, si pone pj = pj + W2 , e si esamina ogni nodo
k adiacente a j. Se il nodo k non è postattivo se ne incrementa la priorità
pk = pk + W2 . Se il nodo k è inattivo viene inserito nella coda di priorità in
stato preattivo ponendo n = n + 1 e qn = k.
10 Uscita dall’algoritmo con la lista η deli nuovi numeri per i nodi del grafo.
Una precisazione finale riguarda i pesi W1 e W2 che riflettono l’importanza della
distanza dal nodo finale rispetto al grado corrente del nodo. Al fine di ottenere i
migliori risultati in [1] vengono suggeriti i seguenti valori: W1 = 1 e W2 = 2.
Implementazione
L’implementazione dell’algoritmo esaminato è stata effettuata mediante il noto linguaggio di programmazione orientato agli oggetti C++ [3, 4]. L’algoritmo di rinumerazione è stato quindi implementato realizzando una serie di classi dotate delle
strutture dati e dei metodi ritenuti più idonei per la risoluzione del problema in
esame. Verranno descritte quindi le classi adottate nell’implementazione a partire
dalla classe fondamentale Graph che formalizza il concetto di grafo non orientato.
– Algoritmo di rinumerazione: implementazione in C++ –
7
La classe Graph incapsula al suo interno tutti i nodi e i lati che definiscono
il grafo stesso e con i suoi metodi pubblici consente: l’aggiunta al grafo di nuovi
nodi e di nuovi lati; l’applicazione di uno degli algoritmi di rinumerazione proposti
rispettivamente da Sloan e da Cuthill e Mckee; il recupero delle nuove numerazioni
generate dagli algoritimi. Per la struttura dati fa uso delle classi GraphNode e
GraphEdge che raccolgono le informazioni relative a ciascun nodo e a ciascun lato.
La classe Graph utilizza la classe RootedLevelStruct che raccoglie le informazioni
relative alle strutture a livelli che vengono di volta in volta generate dall’algoritmo.
Tale compito viene svolto dal metodo privato Graph::GenerateRLS che, a partire
dal nodo radice assegnato, costruisce appunto la struttura a livelli. La classe RootedLevelStruct fa inoltre uso della astrazione intArray, che implementa semplicemente
un vettore di interi, per gestire comodamente i vari livelli nella forma di un array di
array.
In più punti del codice viene utilizzata la classe paramerizzata CArray <type,
type> che è una class di libreria che implementa l’array di oggetti generici e che è
dotata di una serie di metodi standard che ne consentono un facile uso dinamico.
Queste, in sintesi, le caratteristiche principali delle astrazioni adottate. Per
quanto riguarda gli altri dettagli la lettura del codice, che utilizza nomi abbastanza autoesplicativi, consente di riconoscere facilmente il compito svolto da ciascuna
variabile e di rintracciare i passi dell’algoritmo esposto nella sezione precedente.
L’algoritmo di Cuthill e McKee
Anche se nel presente lavoro l’interesse è focalizzato sull’algoritmo di Sloan, coviene
comunque ricordare brevemente l’organizzazione dell’approccio di rinumerazione
dovuto a Cuthill e McKee e del quale nel seguito viene riportata anche un’implementazione in C++.
L’algoritmo comincia selezionando una serie di nodi aventi valori ridotti del grado di connessione e che per tale motivi sono dei buoni candidati dai quali può partire
la numerazione. Per ciascun nodo selezionato i nodi rimanenti vengono rinumerati
seguendo la stuttura a livelli che parte da tale nodo (il nodo di partenza constituisce
da solo il primo livello). Quindi all’interno di ciascun livello i nodi vengono rinumerati per grado di connessione crescente. Il processo di rinumerazione di tutta la
struttura viene ripetuto per tutti i potenziali nodi di partenza precedentemente selezionati scegliendo, alla fine, la sequenza migliore sulla base della massima ampiezza
di banda prodotta.
– Algoritmo di rinumerazione: implementazione in C++ –
8
Risultati numerici
Al fine di mostrare l’efficacia computazionale dell’algoritimo di rinumerazione proposto da Sloan ma soprattutto per testare l’implementazione in C++ qui proposta,
sono stati effettuati alcuni test di confronto sulla base delle geometrie riportate in
figura 1. Quest’ultime si riferiscono a dominii bidimensionali discretizzati mediante
elementi triangolari aventi su ogni nodo tre gradi di libertà, come avviene per esempio nell’analisi di lastre inflesse (l’abbassamento del nodo e le due rotazioni rispetto
agli assi cartesiani contenuti nel piano medio della lastra). Le discretizzazioni sulle
quali sono stati applicati i due algoritmi di rinumerazione, per i risultati ottenuti
si veda la tabella 1, sono state generate per successive suddivisioni in quattro delle
triangolazioni minime necessarie alla definizione della geometria assegnata. Si ottiene in tale modo una mesh più fitta ma caratterizzata da una numerazione molto
casuale dei nodi ed assolutamente inadatta alla soluzione del sistema lineare, come è
anche dimostrato dai valori iniziali del profile, colonna P0 della tabella 1. Per quanto
riguarda i valori assunti da P è utile osservare come essi forniscano direttamente la
quantità memoria accupata dalla matrice rigidezza mediante la semplice formula:
mem = P ∗ sizeof (double)/1000 [Kb].
Dai risultati ottenuti risulta evidente la superiorità dell’algoritmo di Sloan rispetto al tradizionale algoritmo di Cuthill–McKee in termini soprattutto sia di profile
prodotto. Per quanto rigurda i tempi i due approcci mettono in luce una sensibile
differenza a favore del Cuthill–McKee ma solo nel caso in cui il numero di nodi
caratterizzati da un basso valore del grado di connessione è una piccola porzione
del numero totale dei nodi. Infatti per gli ultimi quattro test riportati in tabella
si ottengono dei tempi molto vicini in quanto si ha un numero elevato di nodi di
contorno.
– Algoritmo di rinumerazione: implementazione in C++ –
Mesh
9
Nvar
P0
PCM
PSloan
T empoCM [s]
T empoSloan [s]
S1
75
1761
1014
987
0.002
0.01
S2
243
14220
5382
5193
0.01
0.04
S3
867
144654
33414
32469
0.03
0.16
S4
3247
1732041
231174
226989
0.39
1.26
O1
144
5175
2556
2403
0.002
0.02
O2
480
45024
13947
12849
0.03
0.07
O3
1728
494028
86301
80829
0.10
0.44
O4
6528
6266175
590577
566421
1.51
4.56
H1
135
4716
1971
1764
0.004
0.01
H2
435
38220
11031
9222
0.01
0.07
H3
1534
401022
70263
59526
0.19
0.40
H4
5763
4955343
490263
426246
2.66
4.09
++1
171
7254
2727
2475
0.004
0.02
++2
555
60420
15717
13314
0.03
0.09
++3
1971
644058
101871
88668
0.35
0.60
++4
7395
8049756
720435
646752
5.09
6.51
Tabella 1: Valori di profile ottenuti prima e dopo l’applicazione degli algoritmi di
rinumerazione (I tempi di calcolo riportati sono stati ottenuti compilando il codice in
versione release mediante il compilatore Microsoft visual C++ ver. 6.0 ed eseguendo
i tests su un computer Pentium Pro a 200 MHz).
– Algoritmo di rinumerazione: implementazione in C++ –
Figura 1: Geometrie e mesh adottate.
10
– Algoritmo di rinumerazione: implementazione in C++ –
Il codice
//-----------------------------------------------------------// Graph.h
//
// Interface for the Graph class and other related classes.
//
//-----------------------------------------------------------class GraphNode
{
public:
GraphNode();
GraphNode(GraphNode& gn);
void operator=(GraphNode& gn);
~GraphNode() {}
int degree;
int distance;
int status;
int oldLabel;
double priority;
};
//-----------------------------------------------------------class GraphEdge
{
public:
GraphEdge() {}
GraphEdge(GraphEdge& e);
void operator=(GraphEdge& e);
~GraphEdge() {}
int i;
int j;
};
11
– Algoritmo di rinumerazione: implementazione in C++ –
//-----------------------------------------------------------class intArray
{
public:
intArray()
intArray(intArray& ia)
void operator=(intArray& ia)
int& operator[](int i)
operator const CArray<int,int>&()
int Add(int e)
int GetSize()
void RemoveAll()
void Copy(intArray& ia)
~intArray()
{}
{ array.Copy(ia.array); }
{ array.Copy(ia.array); }
{ return array[i]; }
{ return array; }
{ return array.Add(e); }
{ return array.GetSize(); }
{ array.RemoveAll(); }
{ array.Copy(ia.array); }
{ array.RemoveAll(); }
private:
CArray<int, int> array;
};
//-----------------------------------------------------------class RootedLevelStruct
{
public:
RootedLevelStruct();
RootedLevelStruct(int h);
~RootedLevelStruct();
void operator=(RootedLevelStruct& RLS);
void RemoveAll();
int height;
int width;
CArray<intArray, intArray> L;
CArray<int, int> q;
};
//------------------------------------------------------------
12
– Algoritmo di rinumerazione: implementazione in C++ –
13
class Graph
{
public:
Graph() {}
Graph(int nn, int ne);
~Graph();
BOOL Add(GraphNode& gn);
BOOL Add(GraphEdge& e);
void SetSize(int nn = 0, int ne = 0)
{ nodes.SetSize(nn);
edges.SetSize(ne);
}
BOOL ExistEdge(int i, int j);
void SloanLabelling(CProgressCtrl* prg);
void CutMckLabelling(CProgressCtrl* prg, int initial_band);
void GetNewLabels(CArray<int, int>* lbl);
private:
CArray<GraphNode, GraphNode> nodes;
CArray<GraphEdge, GraphEdge> edges;
CArray<int, int> labels;
int FindSmallestDegree();
BOOL GenerateRLS(int root, RootedLevelStruct* RLS, int xw = 1000000);
};
enum {INACTIVE, PREACTIVE, ACTIVE, POSTACTIVE};
//-----------------------------------------------------------// Graph.cpp
//
// Definition of the declared classes.
//
//-----------------------------------------------------------#include "Sloan.h"
– Algoritmo di rinumerazione: implementazione in C++ –
//-----------------------------------------------------------GraphNode::GraphNode()
{
degree = 0;
distance = 0;
status = INACTIVE;
oldLabel = 0;
priority = 0.0;
}
GraphNode::GraphNode(GraphNode& gn)
{
degree = gn.degree;
distance = gn.distance;
status = gn.status;
oldLabel = gn.oldLabel;
priority = gn.priority;
}
void GraphNode::operator=(GraphNode& gn)
{
degree = gn.degree;
distance = gn.distance;
status = gn.status;
oldLabel = gn.oldLabel;
priority = gn.priority;
}
GraphEdge::GraphEdge(GraphEdge& e)
{
i = e.i;
j = e.j;
}
void GraphEdge::operator=(GraphEdge& e)
{
i = e.i;
14
– Algoritmo di rinumerazione: implementazione in C++ –
j = e.j;
}
//-----------------------------------------------------------RootedLevelStruct::RootedLevelStruct()
{
height = width = 0;
}
RootedLevelStruct::RootedLevelStruct(int h)
{
width = 0;
height = h;
L.SetSize(h);
}
void RootedLevelStruct::operator=(RootedLevelStruct& RLS)
{
RemoveAll();
width = RLS.width;
height = RLS.height;
q.Copy(RLS.q);
L.SetSize(height);
for(int i=0; i<height; i++) L[i].Copy(RLS.L[i]);
}
void RootedLevelStruct::RemoveAll()
{ int i;
for(i=0; i<height; i++) L[i].RemoveAll();
L.RemoveAll();
q.RemoveAll();
height = width = 0;
}
RootedLevelStruct::~RootedLevelStruct()
{
RemoveAll();
15
– Algoritmo di rinumerazione: implementazione in C++ –
}
//-----------------------------------------------------------Graph::Graph(int nn, int ne)
{
nodes.SetSize(nn);
edges.SetSize(ne);
}
Graph::~Graph()
{
nodes.RemoveAll();
edges.RemoveAll();
labels.RemoveAll();
}
BOOL Graph::Add(GraphNode& gn)
{
return nodes.Add(gn);
}
BOOL Graph::Add(GraphEdge& e)
{
nodes[e.i].degree++;
nodes[e.j].degree++;
return edges.Add(e);
}
int Graph::FindSmallestDegree()
{ int nn = nodes.GetSize();
int md = 1000000;
if( !nn ) return -1;
for(int i=0, mi=0; i<nn; i++) if( nodes[i].degree < md )
{ mi = i;
md = nodes[i].degree;
}
return mi;
16
– Algoritmo di rinumerazione: implementazione in C++ –
17
}
BOOL Graph::GenerateRLS(int root, RootedLevelStruct* RLS, int xw)
{ intArray emptylevel, level0, level1;
int k, n, e, nl;
int added = 1, w1 = 0;
int nn = nodes.GetSize(), ne = edges.GetSize(), dst = 0;
for(n=0; n<nn; n++) nodes[n].distance = -1;
RLS->RemoveAll();
RLS->L.Add(emptylevel);
RLS->L[dst].Add(root);
RLS->height = 1;
RLS->width = 1;
nodes[root].distance = dst;
level0.Add(root);
while( added < n )
{ nl = level0.GetSize();
for(k=0; k<nl; k++) for(e=0; e<ne; e++)
{ if( edges[e].i == level0[k] && nodes[edges[e].j].distance == -1 )
{ level1.Add(edges[e].j);
nodes[edges[e].j].distance = dst+1;
w1++;
}
else if( edges[e].j == level0[k] &&
nodes[edges[e].i].distance == -1 )
{ level1.Add(edges[e].i);
nodes[edges[e].i].distance = dst+1;
w1++;
}
else ;
}
if( w1 > xw ) return FALSE;
if( nl = level1.GetSize() )
{ dst++;
RLS->height++;
RLS->width = max(RLS->width, w1);
RLS->L.Add(emptylevel);
level0.RemoveAll();
– Algoritmo di rinumerazione: implementazione in C++ –
18
for(k=0; k<nl; k++)
{ RLS->L[dst].Add(level1[k]);
level0.Add(level1[k]);
added++;
}
}
level1.RemoveAll();
w1 = 0;
}
RLS->q.Copy(RLS->L[dst]);
BOOL bubbled;
nl = RLS->q.GetSize();
do
{ bubbled = FALSE;
for(k=0; k<nl-1; k++)
{ if( nodes[RLS->q[k]].degree > nodes[RLS->q[k+1]].degree )
{ int box = RLS->q[k];
RLS->q[k] = RLS->q[k+1];
RLS->q[k+1] = box;
bubbled = TRUE;
}
else if( nodes[RLS->q[k]].degree == nodes[RLS->q[k+1]].degree )
{ RLS->q.RemoveAt(k--);
nl--;
}
else ;
}
} while( bubbled );
return TRUE;
}
BOOL Graph::ExistEdge(int i, int j)
{ BOOL exist = FALSE;
int e, ne = edges.GetSize();
for(e=0; e<ne; e++)
{ if( (edges[e].i == i && edges[e].j == j) ||
(edges[e].i == j && edges[e].j == i) )
– Algoritmo di rinumerazione: implementazione in C++ –
19
exist = TRUE;
}
return exist;
}
void Graph::GetNewLabels(CArray<int, int>* lbl)
{
lbl->RemoveAll();
lbl->Copy(labels);
}
//-----------------------------------------------------------// Sloan renumbering algorithm
//-----------------------------------------------------------void Graph::SloanLabelling(CProgressCtrl* prg)
{ int start, end, i;
RootedLevelStruct startRLS, endRLS, iRLS;
BOOL found_peripheral = FALSE;
start = FindSmallestDegree();
GenerateRLS(start, &startRLS);
endRLS.width = 1000000;
// Selection of pseudo-peripheral nodes
while( !found_peripheral )
{ prg->SetRange(1, startRLS.q.GetSize());
for(i=0; i<startRLS.q.GetSize(); i++)
{ prg->SetPos(i+1);
if( GenerateRLS(startRLS.q[i], &iRLS, endRLS.width) )
{ if( iRLS.height > startRLS.height && iRLS.width < endRLS.width )
{ start = startRLS.q[i];
startRLS = iRLS;
found_peripheral = FALSE;
break;
}
if( iRLS.width < endRLS.width )
{ end = startRLS.q[i];
endRLS = iRLS;
– Algoritmo di rinumerazione: implementazione in C++ –
found_peripheral = TRUE;
}
}
}
}
GenerateRLS(end, &endRLS);
// Node labelling
CArray<int, int> queue;
double maxp, pr;
double W1 = 1.0, W2 = 2.0;
int i_maxp, ni, nj, nk, lbl;
int ne = edges.GetSize(), nn = nodes.GetSize();
GraphNode *pn;
for(i=0, pn=nodes.GetData(); i<nn; i++, pn++)
{ pn->priority = W1*pn->distance - W2*(pn->degree+1);
pr = pn->priority;
pn->status = INACTIVE;
}
labels.RemoveAll();
labels.SetSize(nn);
lbl = 0;
nodes[start].status = PREACTIVE;
queue.Add(start);
prg->SetRange(1, nn);
while( queue.GetSize() )
{ maxp = -100000;
for(i=0; i<queue.GetSize(); i++)
{ if( maxp < nodes[queue[i]].priority )
{ maxp = nodes[queue[i]].priority;
i_maxp = i;
ni = queue[i];
}
else if( maxp == nodes[queue[i]].priority &&
nodes[queue[i]].degree < nodes[ni].degree )
{ i_maxp = i;
ni = queue[i];
}
20
– Algoritmo di rinumerazione: implementazione in C++ –
else;
}
queue.RemoveAt(i_maxp);
if( nodes[ni].status == PREACTIVE )
{ GraphEdge *pe;
for(i=0, pe=edges.GetData(); i<ne; i++, pe++)
{ if( pe->i == ni ) nj = pe->j;
else if( pe->j == ni ) nj = pe->i;
else nj = -1;
if( nj != -1 )
{ nodes[nj].priority += W2;
if( nodes[nj].status == INACTIVE )
{ queue.Add(nj);
nodes[nj].status = PREACTIVE;
}
}
}
}
prg->SetPos(lbl+1);
labels[lbl] = ni;
lbl++;
nodes[ni].status = POSTACTIVE;
GraphEdge *pe;
for(i=0, pe=edges.GetData(); i<ne; i++, pe++)
{ if( pe->i == ni ) nj = pe->j;
else if( pe->j == ni ) nj = pe->i;
else nj = -1;
if( nj != -1 && nodes[nj].status == PREACTIVE )
{ nodes[nj].status = ACTIVE;
nodes[nj].priority += W2;
GraphEdge *pek;
int k;
for(k=0, pek=edges.GetData(); k<ne; k++, pek++)
{ if( pek->i == nj ) nk = pek->j;
else if( pek->j == nj ) nk = pek->i;
else nk = -1;
if( nk != -1 && nodes[nk].status != POSTACTIVE )
nodes[nk].priority += W2;
21
– Algoritmo di rinumerazione: implementazione in C++ –
if( nk != -1 && nodes[nk].status == INACTIVE )
{ queue.Add(nk);
nodes[nk].status = PREACTIVE;
}
}
}
}
}
}
//-----------------------------------------------------------// Cuthill McKee renumbering algorithm
//-----------------------------------------------------------void Graph::CutMckLabelling(CProgressCtrl* prg, int initial_band)
{ int i, j, ni, nj, n, nn, ne, curr_band, band;
CArray<int, int> level0, level1;
GraphNode *pn;
GraphEdge *pe;
nn = nodes.GetSize();
ne = edges.GetSize();
labels.RemoveAll();
labels.SetSize(nn);
band = initial_band;
prg->SetRange(1, nn);
int aDegree = 0;
for(n=0, pn=nodes.GetData(); n<nn; n++, pn++)
aDegree += pn->degree;
aDegree /= nn;
aDegree--;
for(n=0; n<nn; n++)
{
prg->SetPos(n+1);
int lbl = 0;
if( nodes[n].degree < aDegree )
22
– Algoritmo di rinumerazione: implementazione in C++ –
23
for(i=0, pn=nodes.GetData(); i<nn; i++, pn++)
pn->status = INACTIVE;
curr_band = 0;
level0.RemoveAll();
level0.Add(n);
level1.RemoveAll();
while( level0.GetSize() )
{ int s1, s0 = level0.GetSize();
// rinumerazione del livello corrente e costruzione del successivo
for(i=0; i<s0; i++)
{ ni = level0[i];
nodes[lbl].oldLabel = ni;
nodes[ni].status = POSTACTIVE;
lbl++;
for(j=0, pe=edges.GetData(); j<ne; j++, pe++)
{ if( pe->i == ni ) nj = pe->j;
else if( pe->j == ni ) nj = pe->i;
else nj = -1;
if( nj != -1 && nodes[nj].status == INACTIVE )
{ level1.Add(nj);
nodes[nj].status = ACTIVE;
}
}
}
// controllo sulla convenienza della rinumerazione
curr_band = max(curr_band, s1 = level1.GetSize());
if( curr_band >= band ) break;
// ordinamento del livello successivo per grado di connessione crescente
BOOL bubbled;
do
{ bubbled = TRUE;
int *pl, *pl1;
for(i=0, pl=level1.GetData(), pl1=pl+1; i<s1-1; i++, pl++, pl1++)
if( nodes[*pl1].degree < nodes[*pl].degree )
{ int box = *pl1;
*pl1 = *pl;
*pl = box;
bubbled = FALSE;
– Algoritmo di rinumerazione: implementazione in C++ –
24
}
} while( !bubbled );
level0.RemoveAll();
level0.Copy(level1);
level1.RemoveAll();
}
// rinumerazione conveniente
if( curr_band < band )
{ for(i=0, pn=nodes.GetData(); i<nn; i++, pn++)
labels[i] = pn->oldLabel;
band = curr_band;
}
}
}
//------------------------------------------------------------
Ringraziamenti
Il presente lavoro è un contributo al programma di ricerca Progetto MECOM, Programma Operativo Plurifondo 94/99 - Misura 4.4 ”Ricerca scientifica e tecnologica”.
Riferimenti bibliografici
[1] S. W. Sloan, ‘A FORTRAN program for profile and wavefront reduction’ Int.
J. Num. Meth. Eng., 28, 2651–2679 (1989).
[2] E. Cuthill and J. McKee, ‘Reducing the bandwith of sparse symmetric matrices’
Proc. ACM 24th Nat. Conf. ACM, New York, 157–172 (1969).
[3] B. Stroustrup, Il linguaggio C++, Addison-Wesley, Milano (1990).
[4] G. Booch, Object-oriented anlysis and design, Benjamin/Cummings Publishing
Company, California (1994).