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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1
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 Paso_Pattern* p=A->mainBlock->pattern;
70 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