/[escript]/trunk/paso/src/SystemMatrix_MIS.c
ViewVC logotype

Annotation of /trunk/paso/src/SystemMatrix_MIS.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2966 - (hide annotations)
Wed Mar 3 03:19:39 2010 UTC (9 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 6739 byte(s)
Fixing non-ansi compliant code

1 jfenwick 2965
2     #include "SystemMatrix_MIS.h"
3    
4    
5     #define MISIN 0
6     #define MISOUT 100
7     /* This macro is less complicated than it could have been because it is only used on available nodes */
8     #define ISLESS(x,y) (x<y)
9     #define ISAVAILABLE(x) ((x!=MISIN) && (x!=MISOUT))
10     #define IMAX(x,y) (x>y?x:y)
11    
12     #define MISSTRING(x) ((x==MISIN)?"IN":((x==MISOUT)?"OUT":"UNKNOWN"))
13    
14    
15     /* returns the nodes in this system matrix which connect to things outside the main block
16     the reference parameter count will be set to the number of nodes in the return value.
17     Cleanup of return value is the callers responsibility.
18    
19     This routine assumes that the System matrix is in default format with a CSR col_coupleBlock
20     */
21     index_t* Paso_SparseMatrix_getBorderNodes(Paso_SystemMatrix* A, index_t* count) {
22     const index_t MAXNEIGHBOURS=A->col_coupleBlock->len;
23     index_t* border=MEMALLOC(MAXNEIGHBOURS, index_t);
24     index_t len=0;
25     int i;
26     /* A node is in the border if it has a non-empty row in the coupling block */
27     #pragma omp parallel for schedule(static) private(i)
28     for (i=0;i<A->col_coupleBlock->pattern->numOutput;++i) {
29     if (A->col_coupleBlock->pattern->ptr[i]!=A->col_coupleBlock->pattern->ptr[i+1]) { /* non-empty row */
30     #pragma omp critical
31     {
32     border[len++]=i;
33     }
34     }
35     }
36     *count=len;
37     return border;
38     }
39    
40    
41     /* takes in a list of border nodes and weights and computes the MIS for all border ndoes
42     on all ranks. It does this one node at a time. Later this should use rank colouring
43     to do a number of ranks at once.
44     */
45     void Paso_SystemMatrix_CalcBorderMIS(Paso_SystemMatrix* A, index_t* border, index_t bordercount, double* weights, index_t n) {
46     index_t i=0;
47     index_t j=0;
48     index_t k=0;
49     int mpi_iam = 0; /* rank within world */
50     int mpi_num = 1; /* size of the world */
51     #ifdef PASO_MPI
52     double *remote_values=NULL;
53     #endif
54    
55     if (A->type!=MATRIX_FORMAT_DEFAULT) { /* We only support CSR matricies here */
56     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_CalcBorderMIS: Symmetric matrix patterns are not supported.");
57     }
58    
59     #ifdef PASO_MPI
60     MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
61     MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
62     #endif
63    
64     for (i=0;i<mpi_num;++i) {
65     if (i==mpi_iam) { /* mark my border in preparation for sending to other ranks */
66     for (j=0;j<bordercount;++j) {
67     index_t bnode=border[j];
68     if (ISAVAILABLE(weights[bnode])) {
69 jfenwick 2966 Paso_Pattern* p=A->mainBlock->pattern;
70 jfenwick 2965 weights[bnode]=MISIN;
71     /* Now walk the neighbours and mark them unavailable */
72     for (k=p->ptr[bnode];k<p->ptr[bnode+1];++k) { /* Walk along the row */
73     if (p->index[k]!=bnode) { /* ignore diagonal link to self */
74     weights[p->index[k]]=MISOUT; /* node can't be in the set */
75     }
76     } /* walk row */
77     } /* if avail */
78     } /* walk border */
79     } /* if not me */
80     #ifdef PASO_MPI
81     /* Now we set the weights using the coupler */
82     remote_values=NULL;
83     /* start exchange */
84     Paso_SystemMatrix_startCollect(A,weights);
85     /* finish exchange */
86     remote_values=Paso_SystemMatrix_finishCollect(A);
87     if (i!=mpi_iam) {
88     /* walk over nodes and see if they have a marked neighbour */
89     for (j=0;j<A->col_coupleBlock->pattern->numOutput;++j) {
90     for (k=A->col_coupleBlock->pattern->ptr[j];k<A->col_coupleBlock->pattern->ptr[j+1];++k) {
91     index_t con=A->col_coupleBlock->pattern->index[k];
92     if (remote_values[con]==MISIN) {
93     weights[j]=MISOUT;
94     break; /* done with this row */
95     }
96     }
97     }
98     }
99     #endif
100     } /* End loop over ranks */
101     }
102    
103    
104     /* used to generate pseudo random numbers: */
105    
106     static double Paso_Pattern_mis_seed=.4142135623730951;
107    
108     /* Return a list of nodes which belong the a maximal independent set.
109     Note: Only nodes local to this rank will be returned.
110    
111     Caller is responsible for cleaning up return value.
112     */
113     index_t Paso_SystemMatrix_getMIS(Paso_SystemMatrix* A, index_t** set) {
114     /* identify the nodes on my border */
115     index_t count=0;
116     index_t missize=0;
117     index_t* border=NULL;
118     index_t n=A->mainBlock->numRows;
119     char* inborder=NULL;
120     double* weights=NULL;
121     index_t* mis=NULL;
122     Paso_Pattern* pat=A->mainBlock->pattern;
123     int i,j,k, retry;
124     char done=1;
125     double seed=Paso_Pattern_mis_seed;
126     Paso_resetError();
127     border=Paso_SparseMatrix_getBorderNodes(A, &count);
128     if (!Paso_noError()) {
129     *set=NULL;
130     return 0;
131     }
132    
133     inborder=MEMALLOC(n,char);
134     weights=MEMALLOC(n,double);
135     mis=MEMALLOC(n,index_t);
136    
137     #pragma omp parallel for schedule(static) private(i)
138     for (i=0;i<n;++i) {
139     inborder[i]=0;
140     weights[i]=0.5;
141     }
142     /* This loop will use a different memory access pattern to the setup
143     ie a NUMA problem. However I'm gambling that the border is small and
144     that the memory won't be relocated.
145     */
146     #pragma omp parallel for schedule(static) private(i)
147     for (i=0;i<count;++i) {
148     inborder[border[i]]=1;
149     }
150    
151     /* compute set membership for border nodes */
152     Paso_SystemMatrix_CalcBorderMIS(A, border, count, weights, n);
153    
154     do {
155     done=1;
156     retry=0;
157     j=0;
158     /* Randomise weights of unmarked nodes */
159     #pragma omp parallel for schedule(static) private(i)
160     for (i=0;i<n;++i) {
161     if (ISAVAILABLE(weights[i])) {
162     weights[i]=fmod(seed*(i+1),1.)+0.1;
163     ++j;
164     }
165     }
166     k=0;
167     /* Now we evaluate nodes relative to their neighbours */
168     #pragma omp parallel for schedule(static) private(i)
169     for (i=0;i<n;++i) {
170     char ok=1;
171     for (j=pat->ptr[i];j<pat->ptr[i+1];++j) {
172     if (i!=pat->index[j] && !ISLESS(weights[i],weights[pat->index[j]])) {
173     ok=0;
174     break;
175     }
176     }
177     if (ok) {
178     weights[i]=MISIN;
179     k++;
180     retry=1; /* multiple threads writing same value should be fine */
181     }
182     }
183    
184     /* Go through and mark all the neighbours of nodes definitly in */
185     #pragma omp parallel for schedule(static) private(i)
186     for (i=0;i<n;++i) {
187     if (weights[i]==MISIN) {
188     for (j=pat->ptr[i];j<pat->ptr[i+1];++j) {
189     if (pat->index[j]!=i) {
190     weights[pat->index[j]]=MISOUT;
191     }
192     }
193     }
194     }
195    
196     /* Are we finished yet? */
197     #pragma omp parallel for schedule(static) private(i)
198     for (i=0;i<n;++i) {
199     if ((weights[i]!=MISIN) && (weights[i]!=MISOUT)) { /* neither in or neighbour of in */
200     done=0; /* multiple threads writing same value should be fine */
201     }
202     }
203     seed+=0.014159; /* Arbitrary change to seed */
204     } while (!done && retry);
205    
206     for (i=0;i<n;++i) {
207     if (weights[i]==MISIN) {
208     mis[missize++]=i;
209     }
210     }
211     MEMFREE(border);
212     MEMFREE(weights);
213     MEMFREE(inborder);
214     if (done==0) {
215     Paso_setError(NO_PROGRESS_ERROR,"Error in MIS - no progress.");
216     MEMFREE(mis);
217     *set=NULL;
218     return 0;
219     }
220     *set=mis;
221     return missize;
222     }

  ViewVC Help
Powered by ViewVC 1.1.26