1 /* This software was developed by Bruce Hendrickson and Robert Leland   *
2  * at Sandia National Laboratories under US Department of Energy        *
3  * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
4 
5 #include <stdio.h>
6 #include "params.h"
7 #include "defs.h"
8 #include "structs.h"
9 
10 int LIMIT_KL_EWGTS = FALSE;	/* limit range of edge weights in multilevel-KL? */
11 double EWGT_RATIO_MAX = .25;	/* if so, max allowed ewgt/nvtxs */
12 
coarsen_kl(graph,nvtxs,nedges,using_vwgts,using_ewgts,term_wgts,igeom,coords,vwgt_max,assignment,goal,architecture,hops,solver_flag,ndims,nsets,vmax,mediantype,mkconnected,eigtol,nstep,step,pbndy_list,weights,give_up)13 void      coarsen_kl(graph, nvtxs, nedges, using_vwgts, using_ewgts, term_wgts,
14 		     igeom, coords, vwgt_max, assignment, goal, architecture, hops,
15 		          solver_flag, ndims, nsets, vmax, mediantype, mkconnected,
16 		               eigtol, nstep, step, pbndy_list, weights, give_up)
17 /* Coarsen until nvtxs < vmax, compute and uncoarsen. */
18 struct vtx_data **graph;	/* array of vtx data for graph */
19 int       nvtxs;		/* number of vertices in graph */
20 int       nedges;		/* number of edges in graph */
21 int       using_vwgts;		/* are vertices weights being used? */
22 int       using_ewgts;		/* are edge weights being used? */
23 float    *term_wgts[];		/* weights for terminal propogation */
24 int       igeom;                /* dimension for geometric information */
25 float   **coords;               /* coordinates for vertices */
26 int       vwgt_max;		/* largest vertex weight */
27 short    *assignment;		/* processor each vertex gets assigned to */
28 double   *goal;			/* desired set sizes */
29 int       architecture;		/* 0 => hypercube, d => d-dimensional mesh */
30 short     (*hops)[MAXSETS];	/* cost of edge between sets */
31 int       solver_flag;		/* which eigensolver to use */
32 int       ndims;		/* number of eigenvectors to calculate */
33 int       nsets;		/* number of sets being divided into */
34 int       vmax;			/* largest subgraph to stop coarsening */
35 int       mediantype;		/* flag for different assignment strategies */
36 int       mkconnected;		/* make graph connected before eigensolver? */
37 double    eigtol;		/* tolerence in eigen calculation */
38 int       nstep;		/* number of coarsenings between RQI steps */
39 int       step;			/* current step number */
40 int     **pbndy_list;		/* pointer to returned boundary list */
41 double   *weights;		/* weights of vertices in each set */
42 int       give_up;		/* has coarsening bogged down? */
43 {
44     extern FILE *Output_File;	/* output file or null */
45     extern int DEBUG_TRACE;	/* trace the execution of the code */
46     extern int DEBUG_COARSEN;	/* debug flag for coarsening */
47     extern int DEBUG_CONNECTED;	/* debug flag for connectivity checking */
48     extern double COARSEN_RATIO_MIN;	/* vtx reduction demanded */
49     extern int COARSEN_VWGTS;	/* use vertex weights while coarsening? */
50     extern int COARSEN_EWGTS;	/* use edge weights while coarsening? */
51     extern int COARSE_KL_BOTTOM;/* force KL invocation at bottom level? */
52     extern int KL_ONLY_BNDY;	/* only initialize boundary vertices? */
53     extern double KL_IMBALANCE;	/* fractional imbalance allowed in KL */
54     extern int MAPPING_TYPE;	/* how to get from eigenvectors to partition */
55     struct connect_data *cdata;	/* data structure for enforcing connectivity */
56     struct vtx_data **cgraph;	/* array of vtx data for coarsened graph */
57     double   *yvecs[MAXDIMS + 1];	/* eigenvectors for subgraph */
58     double    evals[MAXDIMS + 1];	/* eigenvalues returned */
59     double    new_goal[MAXSETS];/* new goal if not using vertex weights */
60     double   *real_goal;	/* chooses between goal and new_goal */
61     double    goal_weight;	/* total weight of vertices in goal */
62     double   *vwsqrt = NULL;	/* square root of vertex weights */
63     double    maxdeg;		/* maximum weighted degree of a vertex */
64     double    temp_goal[2];     /* goal to simulate bisection while striping */
65     double   *fake_goal;        /* either goal or temp_goal */
66     float    *cterm_wgts[MAXSETS];	/* terminal weights for coarse graph */
67     float    *new_term_wgts[MAXSETS];	/* modified for Bui's method */
68     float   **real_term_wgts;	/* which of previous two to use */
69     float     ewgt_max;		/* largest edge weight in graph */
70     float    *twptr;		/* loops through term_wgts */
71     float    *twptr_save;	/* copy of twptr */
72     float    *ctwptr;		/* loops through cterm_wgts */
73     float   **ccoords;		/* coarse graph coordinates */
74     int      *active;		/* space for assign routine */
75     int      *v2cv;		/* mapping from fine to coarse vertices */
76     int      *cv2v;		/* mapping from coarse to fine vertices */
77     int      *bndy_list;	/* list of vertices on boundary */
78     int      *cbndy_list;	/* list of vertices of coarse graph on boundary */
79     int      *mflag;		/* flags indicating matching */
80     short    *cassignment;	/* set assignments for coarsened vertices */
81     int       clist_length;	/* length of coarse graph boundary vtx list */
82     int       list_length;	/* length of boundary vtx list */
83     int       vtx;		/* vertex in graph */
84     int       cvtx;		/* vertex in coarse graph */
85     int       cnvtxs;		/* number of vertices in coarsened graph */
86     int       cnedges;		/* number of edges in coarsened graph */
87     int       cvwgt_max;	/* largest vertex weight in coarsened graph */
88     int       nextstep;		/* next step in RQI test */
89     int       max_dev;		/* largest allowed deviation from balance */
90     int       set;		/* set a vertex is in */
91     int       i, j;		/* loop counters */
92     int       sfree();
93     double   *smalloc(), *srealloc(), find_maxdeg();
94     void      makevwsqrt(), make_connected(), print_connected(), eigensolve();
95     void      make_unconnected(), assign(), klspiff(), coarsen1(), free_graph();
96     void      compress_ewgts(), restore_ewgts(), count_weights();
97 
98     if (DEBUG_COARSEN > 0 || DEBUG_TRACE > 0) {
99 	printf("<Entering coarsen_kl, step=%d, nvtxs=%d, nedges=%d, vmax=%d>\n",
100 	       step, nvtxs, nedges, vmax);
101     }
102 
103     /* Is problem small enough to solve? */
104     if (nvtxs <= vmax || give_up) {
105 	if (using_vwgts) {
106 	    vwsqrt = (double *) smalloc((unsigned) (nvtxs + 1) * sizeof(double));
107 	    makevwsqrt(vwsqrt, graph, nvtxs);
108 	}
109 	else
110 	    vwsqrt = NULL;
111 
112 	/* Create space for subgraph yvecs. */
113 	for (i = 1; i <= ndims; i++) {
114 	    yvecs[i] = (double *) smalloc((unsigned) (nvtxs + 1) * sizeof(double));
115 	}
116 
117 	active = (int *) smalloc((unsigned) nvtxs * sizeof(int));
118 	if (mkconnected) {
119 	    make_connected(graph, nvtxs, &nedges, (short *) &(yvecs[1][0]),
120 			   active, &cdata, using_ewgts);
121 	    if (DEBUG_CONNECTED > 0) {
122 		printf("Enforcing connectivity on coarse graph\n");
123 		print_connected(cdata);
124 	    }
125 	}
126 
127 	/* If not coarsening ewgts, then need care with term_wgts. */
128 	if (!using_ewgts && term_wgts[1] != NULL && step != 0) {
129 	    twptr = (float *) smalloc((unsigned)
130 		(nvtxs + 1) * (nsets - 1) * sizeof(float));
131 	    twptr_save = twptr;
132 	    for (j = 1; j < nsets; j++) {
133 	        new_term_wgts[j] = twptr;
134 	        twptr += nvtxs + 1;
135 	    }
136 
137 	    for (j = 1; j < nsets; j++) {
138 	        twptr = term_wgts[j];
139 	        ctwptr = new_term_wgts[j];
140 	        for (i = 1; i <= nvtxs; i++) {
141 		    if (twptr[i] > .5) ctwptr[i] = 1;
142 		    else if (twptr[i] < -.5) ctwptr[i] = -1;
143 		    else ctwptr[i] = 0;
144 		}
145 	    }
146 	    real_term_wgts = new_term_wgts;
147 	}
148 	else {
149 	    real_term_wgts = term_wgts;
150 	    new_term_wgts[1] = NULL;
151 	}
152 
153 	if (!COARSEN_VWGTS && step != 0) {	/* Construct new goal */
154 	    goal_weight = 0;
155 	    for (i = 0; i < nsets; i++)
156 		goal_weight += goal[i];
157 	    for (i = 0; i < nsets; i++)
158 		new_goal[i] = goal[i] * (nvtxs / goal_weight);
159 	    real_goal = new_goal;
160 	}
161 	else
162 	    real_goal = goal;
163 
164 	if (ndims == 1 && nsets > 2) {		/* Striping. */
165 	    mediantype = 4;
166 	    temp_goal[0] = temp_goal[1] = 0;
167 	    for (i = 0; 2 * i + 1 < nsets; i++) {
168 		 temp_goal[0] += real_goal[i];
169 		 temp_goal[1] += real_goal[nsets - 1 - i];
170 	    }
171 	    i = nsets / 2;
172 	    if (2 * i != nsets) {
173 		temp_goal[0] += .5 * real_goal[i];
174 		temp_goal[1] += .5 * real_goal[i];
175 	    }
176 	    fake_goal = temp_goal;
177 	}
178 	else {
179 	    mediantype = MAPPING_TYPE;
180 	    fake_goal = real_goal;
181 	}
182 
183 	/* Partition coarse graph with spectral method */
184 	maxdeg = find_maxdeg(graph, nvtxs, using_ewgts, &ewgt_max);
185 	eigensolve(graph, nvtxs, nedges, maxdeg, vwgt_max, vwsqrt,
186 		   using_vwgts, using_ewgts, real_term_wgts, 0, (float **) NULL,
187 		   yvecs, evals, architecture, assignment, fake_goal,
188 		   solver_flag, FALSE, 0, ndims, mediantype, eigtol);
189 
190 	if (mkconnected)
191 	    make_unconnected(graph, &nedges, &cdata, using_ewgts);
192 
193 	assign(graph, yvecs, nvtxs, ndims, architecture, nsets, vwsqrt, assignment,
194 	       active, mediantype, real_goal, vwgt_max);
195 /*
196 simple_part(graph, nvtxs, assignment, nsets, 1, real_goal);
197 */
198 	sfree((char *) active);
199 
200 	count_weights(graph, nvtxs, assignment, nsets, weights, (vwgt_max != 1));
201 
202 	bndy_list = NULL;
203 	if (COARSE_KL_BOTTOM || !(step % nstep)) {
204 	    if (LIMIT_KL_EWGTS)
205 		compress_ewgts(graph, nvtxs, nedges, ewgt_max, using_ewgts);
206 
207 	    max_dev = (step == 0) ? vwgt_max : 5 * vwgt_max;
208 	    goal_weight = 0;
209 	    for (i = 0; i < nsets; i++) {
210 		goal_weight += real_goal[i];
211 	    }
212 	    goal_weight *= KL_IMBALANCE / nsets;
213 	    if (goal_weight > max_dev) {
214 		max_dev = goal_weight;
215 	    }
216 
217 	    if (KL_ONLY_BNDY) {
218 		/* On small graph, quickest to include everybody in list. */
219 		bndy_list = (int *) smalloc((unsigned) (nvtxs+1)*sizeof(int));
220 		for (i=0; i<nvtxs; i++) {
221 		    bndy_list[i] = i+1;
222 		}
223 		bndy_list[nvtxs] = 0;
224 	    }
225 	    klspiff(graph, nvtxs, assignment, nsets, hops, real_goal,
226 		    real_term_wgts, max_dev, maxdeg, using_ewgts,
227 		    &bndy_list, weights);
228 	    if (LIMIT_KL_EWGTS)
229 		restore_ewgts(graph, nvtxs);
230 	}
231 
232 	else if (KL_ONLY_BNDY) {	/* Find boundary vertices directly. */
233 	    bndy_list = (int *) smalloc((unsigned) (nvtxs+1)*sizeof(int));
234 	    list_length = 0;
235 	    for (i=1; i<=nvtxs; i++) {
236 		set = assignment[i];
237 		for (j=1; j<graph[i]->nedges; j++) {
238 		    if (assignment[graph[i]->edges[j]] != set) {
239 			bndy_list[list_length++] = i;
240 			break;
241 		    }
242 		}
243 	    }
244 	    bndy_list[list_length] = 0;
245 	    bndy_list = (int *) srealloc((char *) bndy_list,
246 			(unsigned) (list_length + 1) * sizeof(int));
247 	}
248 	*pbndy_list = bndy_list;
249 
250 	if (real_term_wgts != term_wgts && new_term_wgts[1] != NULL) {
251 	    sfree((char *) real_term_wgts[1]);
252 	}
253 	if (vwsqrt != NULL)
254 	    sfree((char *) vwsqrt);
255 	for (i = ndims; i > 0; i--)
256 	    sfree((char *) yvecs[i]);
257 	return;
258     }
259 
260     /* Otherwise I have to coarsen. */
261     if (coords != NULL) {
262         ccoords = (float **) smalloc((unsigned) igeom * sizeof(float *));
263     }
264     else {
265         ccoords = NULL;
266     }
267     coarsen1(graph, nvtxs, nedges, &cgraph, &cnvtxs, &cnedges, &v2cv,
268 	     igeom, coords, ccoords, using_ewgts);
269 
270     if (term_wgts[1] != NULL) {
271 	twptr = (float *) smalloc((unsigned)
272 	    (cnvtxs + 1) * (nsets - 1) * sizeof(float));
273 	twptr_save = twptr;
274 	for (i = (cnvtxs + 1) * (nsets - 1); i ; i--) {
275 	    *twptr++ = 0;
276 	}
277 	twptr = twptr_save;
278 	for (j = 1; j < nsets; j++) {
279 	    cterm_wgts[j] = twptr;
280 	    twptr += cnvtxs + 1;
281 	}
282 	for (j = 1; j < nsets; j++) {
283 	    ctwptr = cterm_wgts[j];
284 	    twptr = term_wgts[j];
285 	    for (i = 1; i < nvtxs; i++){
286 	        ctwptr[v2cv[i]] += twptr[i];
287 	    }
288 	}
289     }
290     else {
291 	cterm_wgts[1] = NULL;
292     }
293 
294     /* If coarsening isn't working very well, give up and partition. */
295     give_up = FALSE;
296     if (nvtxs * COARSEN_RATIO_MIN < cnvtxs && cnvtxs > vmax) {
297 	printf("WARNING: Coarsening not making enough progress, nvtxs = %d, cnvtxs = %d.\n",
298 	    nvtxs, cnvtxs);
299 	printf("         Recursive coarsening being stopped prematurely.\n");
300 	if (Output_File != NULL) {
301 	    fprintf(Output_File,
302 		"WARNING: Coarsening not making enough progress, nvtxs = %d, cnvtxs = %d.\n",
303 	    nvtxs, cnvtxs);
304 	    fprintf(Output_File,
305 		"         Recursive coarsening being stopped prematurely.\n");
306 	}
307 	give_up = TRUE;
308     }
309 
310     /* Now recurse on coarse subgraph. */
311     if (COARSEN_VWGTS) {
312 	cvwgt_max = 0;
313 	for (i = 1; i <= cnvtxs; i++) {
314 	    if (cgraph[i]->vwgt > cvwgt_max)
315 		cvwgt_max = cgraph[i]->vwgt;
316 	}
317     }
318     else
319 	cvwgt_max = 1;
320 
321     cassignment = (short *) smalloc((unsigned) (cnvtxs + 1) * sizeof(short));
322     nextstep = step + 1;
323     coarsen_kl(cgraph, cnvtxs, cnedges, COARSEN_VWGTS, COARSEN_EWGTS, cterm_wgts,
324 	       igeom, ccoords, cvwgt_max, cassignment, goal, architecture, hops,
325 	       solver_flag, ndims, nsets, vmax, mediantype, mkconnected,
326 	       eigtol, nstep, nextstep, &cbndy_list, weights, give_up);
327 
328     /* Interpolate assignment back to fine graph. */
329     for (i = 1; i <= nvtxs; i++)
330 	assignment[i] = cassignment[v2cv[i]];
331 
332     if (KL_ONLY_BNDY) {
333 	/* Construct boundary list from coarse boundary list. */
334 	cv2v = (int *) smalloc((unsigned) (cnvtxs + 1) * sizeof(int));
335         mflag = (int *) smalloc((unsigned) (nvtxs + 1) * sizeof(int));
336 	for (i=1; i<=cnvtxs; i++) {
337 	    cv2v[i] = 0;
338 	}
339 	for (i=1; i<=nvtxs; i++) {
340 	    mflag[i] = 0;
341 	    cvtx = v2cv[i];
342 	    if (cv2v[cvtx] == 0) {
343 		cv2v[cvtx] = i;
344 	    }
345 	    else {
346 		mflag[cv2v[cvtx]] = i;
347 	    }
348 	}
349 
350 	clist_length = 0;
351 	while (cbndy_list[clist_length] != 0) {
352 	    clist_length++;
353 	}
354 	bndy_list = (int *) smalloc((unsigned) (2 * clist_length + 1) * sizeof(int));
355 
356 	list_length = 0;
357 	for (i=0; i<clist_length; i++) {
358 	    vtx = cv2v[cbndy_list[i]];
359 	    bndy_list[list_length++] = vtx;
360 	    if (mflag[vtx] != 0) {
361 		bndy_list[list_length++] = mflag[vtx];
362 	    }
363 	}
364 	bndy_list[list_length] = 0;
365 	bndy_list = (int *) srealloc((char *) bndy_list,
366 			(unsigned) (list_length + 1) * sizeof(int));
367 
368         sfree((char *) mflag);
369 	sfree((char *) cv2v);
370 	sfree((char *) cbndy_list);
371     }
372     else {
373 	bndy_list = NULL;
374     }
375 
376     /* Free the space that was allocated. */
377     sfree((char *) cassignment);
378     if (cterm_wgts[1] != NULL)
379 	sfree((char *) cterm_wgts[1]);
380     free_graph(cgraph);
381     sfree((char *) v2cv);
382 
383     /* Smooth using KL every nstep steps. */
384     if (!(step % nstep)) {
385 	if (!COARSEN_VWGTS && step != 0) {	/* Construct new goal */
386 	    goal_weight = 0;
387 	    for (i = 0; i < nsets; i++)
388 		goal_weight += goal[i];
389 	    for (i = 0; i < nsets; i++)
390 		new_goal[i] = goal[i] * (nvtxs / goal_weight);
391 	    real_goal = new_goal;
392 	}
393 	else
394 	    real_goal = goal;
395 
396 	maxdeg = find_maxdeg(graph, nvtxs, using_ewgts, &ewgt_max);
397 	if (LIMIT_KL_EWGTS)
398 	    compress_ewgts(graph, nvtxs, nedges, ewgt_max,
399 			   using_ewgts);
400 
401 	/* If not coarsening ewgts, then need care with term_wgts. */
402 	if (!using_ewgts && term_wgts[1] != NULL && step != 0) {
403 	    twptr = (float *) smalloc((unsigned)
404 		(nvtxs + 1) * (nsets - 1) * sizeof(float));
405 	    twptr_save = twptr;
406 	    for (j = 1; j < nsets; j++) {
407 	        new_term_wgts[j] = twptr;
408 	        twptr += nvtxs + 1;
409 	    }
410 
411 	    for (j = 1; j < nsets; j++) {
412 	        twptr = term_wgts[j];
413 	        ctwptr = new_term_wgts[j];
414 	        for (i = 1; i <= nvtxs; i++) {
415 		    if (twptr[i] > .5) ctwptr[i] = 1;
416 		    else if (twptr[i] < -.5) ctwptr[i] = -1;
417 		    else ctwptr[i] = 0;
418 		}
419 	    }
420 	    real_term_wgts = new_term_wgts;
421 	}
422 	else {
423 	    real_term_wgts = term_wgts;
424 	    new_term_wgts[1] = NULL;
425 	}
426 
427 	max_dev = (step == 0) ? vwgt_max : 5 * vwgt_max;
428 	goal_weight = 0;
429 	for (i = 0; i < nsets; i++) {
430 	    goal_weight += real_goal[i];
431 	}
432 	goal_weight *= KL_IMBALANCE / nsets;
433 	if (goal_weight > max_dev) {
434 	    max_dev = goal_weight;
435 	}
436 
437 	if (!COARSEN_VWGTS) {
438 	    count_weights(graph, nvtxs, assignment, nsets, weights,
439 		(vwgt_max != 1));
440 	}
441 
442 	klspiff(graph, nvtxs, assignment, nsets, hops, real_goal,
443 		real_term_wgts, max_dev, maxdeg, using_ewgts,
444 		&bndy_list, weights);
445 
446 	if (real_term_wgts != term_wgts && new_term_wgts[1] != NULL) {
447 	    sfree((char *) real_term_wgts[1]);
448 	}
449 
450 	if (LIMIT_KL_EWGTS)
451 	    restore_ewgts(graph, nvtxs);
452     }
453     *pbndy_list = bndy_list;
454 
455     if (ccoords != NULL) {
456         for (i = 0; i < igeom; i++)
457             sfree((char *) ccoords[i]);
458         sfree((char *) ccoords);
459     }
460 
461     if (DEBUG_COARSEN > 0) {
462 	printf(" Leaving coarsen_kl, step=%d\n", step);
463     }
464 }
465