/[escript]/branches/doubleplusgood/paso/src/SystemMatrix_MIS.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SystemMatrix_MIS.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4275 - (show annotations)
Tue Mar 5 09:32:52 2013 UTC (6 years, 9 months ago) by jfenwick
File size: 6723 byte(s)
some experiments
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 caller's 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=new index_t[MAXNEIGHBOURS];
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
42 * border nodes on all ranks. It does this one node at a time. Later this
43 * should use rank colouring 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 ESYS_MPI
52 double *remote_values=NULL;
53 #endif
54
55 if (A->type!=MATRIX_FORMAT_DEFAULT) { /* We only support CSR matrices here */
56 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_CalcBorderMIS: Symmetric matrix patterns are not supported.");
57 }
58
59 #ifdef ESYS_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 ESYS_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 to 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 Esys_resetError();
127 border=Paso_SparseMatrix_getBorderNodes(A, &count);
128 if (!Esys_noError()) {
129 *set=NULL;
130 return 0;
131 }
132
133 inborder=new char[n];
134 weights=new double[n];
135 mis=new index_t[n];
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 i.e. 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 definitely 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 delete[] border;
212 delete[] weights;
213 delete[] inborder;
214 if (done==0) {
215 Esys_setError(NO_PROGRESS_ERROR,"Error in MIS - no progress.");
216 delete[] mis;
217 *set=NULL;
218 return 0;
219 }
220 *set=mis;
221 return missize;
222 }
223

  ViewVC Help
Powered by ViewVC 1.1.26