1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include "Multilevel.h"
15#include "PriorityQueue.h"
16#include "memory.h"
17#include "logic.h"
18#include "assert.h"
19#include "arith.h"
20
21
22Multilevel_control Multilevel_control_new(int scheme, int mode){
23 Multilevel_control ctrl;
24
25 ctrl = GNEW(struct Multilevel_control_struct);
26 ctrl->minsize = 4;
27 ctrl->min_coarsen_factor = 0.75;
28 ctrl->maxlevel = 1<<30;
29 ctrl->randomize = TRUE;
30 /* now set in spring_electrical_control_new(), as well as by command line argument -c
31 ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
32 ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET_RS;
33 ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE;
34 ctrl->coarsen_scheme = COARSEN_HYBRID;
35 ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST;
36 ctrl->coarsen_mode = COARSEN_MODE_FORCEFUL; or COARSEN_MODE_GENTLE;
37 */
38
39 ctrl->coarsen_scheme = scheme;
40 ctrl->coarsen_mode = mode;
41 return ctrl;
42}
43
44void Multilevel_control_delete(Multilevel_control ctrl){
45 FREE(ctrl);
46}
47
48static Multilevel Multilevel_init(SparseMatrix A, SparseMatrix D, real *node_weights){
49 Multilevel grid;
50 if (!A) return NULL;
51 assert(A->m == A->n);
52 grid = GNEW(struct Multilevel_struct);
53 grid->level = 0;
54 grid->n = A->n;
55 grid->A = A;
56 grid->D = D;
57 grid->P = NULL;
58 grid->R = NULL;
59 grid->node_weights = node_weights;
60 grid->next = NULL;
61 grid->prev = NULL;
62 grid->delete_top_level_A = FALSE;
63 return grid;
64}
65
66void Multilevel_delete(Multilevel grid){
67 if (!grid) return;
68 if (grid->A){
69 if (grid->level == 0) {
70 if (grid->delete_top_level_A) {
71 SparseMatrix_delete(grid->A);
72 if (grid->D) SparseMatrix_delete(grid->D);
73 }
74 } else {
75 SparseMatrix_delete(grid->A);
76 if (grid->D) SparseMatrix_delete(grid->D);
77 }
78 }
79 SparseMatrix_delete(grid->P);
80 SparseMatrix_delete(grid->R);
81 if (grid->node_weights && grid->level > 0) FREE(grid->node_weights);
82 Multilevel_delete(grid->next);
83 FREE(grid);
84}
85
86static void maximal_independent_vertex_set(SparseMatrix A, int randomize, int **vset, int *nvset, int *nzc){
87 int i, ii, j, *ia, *ja, m, n, *p = NULL;
88 assert(A);
89 assert(SparseMatrix_known_strucural_symmetric(A));
90 ia = A->ia;
91 ja = A->ja;
92 m = A->m;
93 n = A->n;
94 assert(n == m);
95 *vset = N_GNEW(m,int);
96 for (i = 0; i < m; i++) (*vset)[i] = MAX_IND_VTX_SET_U;
97 *nvset = 0;
98 *nzc = 0;
99
100 if (!randomize){
101 for (i = 0; i < m; i++){
102 if ((*vset)[i] == MAX_IND_VTX_SET_U){
103 (*vset)[i] = (*nvset)++;
104 for (j = ia[i]; j < ia[i+1]; j++){
105 if (i == ja[j]) continue;
106 (*vset)[ja[j]] = MAX_IND_VTX_SET_F;
107 (*nzc)++;
108 }
109 }
110 }
111 } else {
112 p = random_permutation(m);
113 for (ii = 0; ii < m; ii++){
114 i = p[ii];
115 if ((*vset)[i] == MAX_IND_VTX_SET_U){
116 (*vset)[i] = (*nvset)++;
117 for (j = ia[i]; j < ia[i+1]; j++){
118 if (i == ja[j]) continue;
119 (*vset)[ja[j]] = MAX_IND_VTX_SET_F;
120 (*nzc)++;
121 }
122 }
123 }
124 FREE(p);
125 }
126 (*nzc) += *nvset;
127}
128
129
130static void maximal_independent_vertex_set_RS(SparseMatrix A, int randomize, int **vset, int *nvset, int *nzc){
131 /* The Ruge-Stuben coarsening scheme. Initially all vertices are in the U set (with marker MAX_IND_VTX_SET_U),
132 with gain equal to their degree. Select vertex with highest gain into a C set (with
133 marker >= MAX_IND_VTX_SET_C), and their neighbors j in F set (with marker MAX_IND_VTX_SET_F). The neighbors of
134 j that are in the U set get their gains incremented by 1. So overall
135 gain[k] = |{neighbor of k in U set}|+2*|{neighbors of k in F set}|.
136 nzc is the number of entries in the restriction matrix
137 */
138 int i, jj, ii, *p = NULL, j, k, *ia, *ja, m, n, gain, removed, nf = 0;
139 PriorityQueue q;
140 assert(A);
141 assert(SparseMatrix_known_strucural_symmetric(A));
142
143 ia = A->ia;
144 ja = A->ja;
145 m = A->m;
146 n = A->n;
147 assert(n == m);
148 *vset = N_GNEW(m,int);
149 for (i = 0; i < m; i++) {
150 (*vset)[i] = MAX_IND_VTX_SET_U;
151 }
152 *nvset = 0;
153 *nzc = 0;
154
155 q = PriorityQueue_new(m, 2*(m-1));
156
157 if (!randomize){
158 for (i = 0; i < m; i++)
159 PriorityQueue_push(q, i, ia[i+1] - ia[i]);
160 } else {
161 p = random_permutation(m);
162 for (ii = 0; ii < m; ii++){
163 i = p[ii];
164 PriorityQueue_push(q, i, ia[i+1] - ia[i]);
165 }
166 FREE(p);
167 }
168
169 while (PriorityQueue_pop(q, &i, &gain)){
170 assert((*vset)[i] == MAX_IND_VTX_SET_U);
171 (*vset)[i] = (*nvset)++;
172 for (j = ia[i]; j < ia[i+1]; j++){
173 jj = ja[j];
174 assert((*vset)[jj] == MAX_IND_VTX_SET_U || (*vset)[jj] == MAX_IND_VTX_SET_F);
175 if (i == jj) continue;
176
177 if ((*vset)[jj] == MAX_IND_VTX_SET_U){
178 removed = PriorityQueue_remove(q, jj);
179 assert(removed);
180 (*vset)[jj] = MAX_IND_VTX_SET_F;
181 nf++;
182
183 for (k = ia[jj]; k < ia[jj+1]; k++){
184 if (jj == ja[k]) continue;
185 if ((*vset)[ja[k]] == MAX_IND_VTX_SET_U){
186 gain = PriorityQueue_get_gain(q, ja[k]);
187 assert(gain >= 0);
188 PriorityQueue_push(q, ja[k], gain + 1);
189 }
190 }
191 }
192 (*nzc)++;
193 }
194 }
195 (*nzc) += *nvset;
196 PriorityQueue_delete(q);
197
198}
199
200
201
202static void maximal_independent_edge_set(SparseMatrix A, int randomize, int **matching, int *nmatch){
203 int i, ii, j, *ia, *ja, m, n, *p = NULL;
204 assert(A);
205 assert(SparseMatrix_known_strucural_symmetric(A));
206 ia = A->ia;
207 ja = A->ja;
208 m = A->m;
209 n = A->n;
210 assert(n == m);
211 *matching = N_GNEW(m,int);
212 for (i = 0; i < m; i++) (*matching)[i] = i;
213 *nmatch = n;
214
215 if (!randomize){
216 for (i = 0; i < m; i++){
217 for (j = ia[i]; j < ia[i+1]; j++){
218 if (i == ja[j]) continue;
219 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
220 (*matching)[ja[j]] = i;
221 (*matching)[i] = ja[j];
222 (*nmatch)--;
223 }
224 }
225 }
226 } else {
227 p = random_permutation(m);
228 for (ii = 0; ii < m; ii++){
229 i = p[ii];
230 for (j = ia[i]; j < ia[i+1]; j++){
231 if (i == ja[j]) continue;
232 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
233 (*matching)[ja[j]] = i;
234 (*matching)[i] = ja[j];
235 (*nmatch)--;
236 }
237 }
238 }
239 FREE(p);
240 }
241}
242
243
244
245static void maximal_independent_edge_set_heavest_edge_pernode(SparseMatrix A, int randomize, int **matching, int *nmatch){
246 int i, ii, j, *ia, *ja, m, n, *p = NULL;
247 real *a, amax = 0;
248 int first = TRUE, jamax = 0;
249
250 assert(A);
251 assert(SparseMatrix_known_strucural_symmetric(A));
252 ia = A->ia;
253 ja = A->ja;
254 m = A->m;
255 n = A->n;
256 assert(n == m);
257 *matching = N_GNEW(m,int);
258 for (i = 0; i < m; i++) (*matching)[i] = i;
259 *nmatch = n;
260
261 assert(SparseMatrix_is_symmetric(A, FALSE));
262 assert(A->type == MATRIX_TYPE_REAL);
263
264 a = (real*) A->a;
265 if (!randomize){
266 for (i = 0; i < m; i++){
267 first = TRUE;
268 for (j = ia[i]; j < ia[i+1]; j++){
269 if (i == ja[j]) continue;
270 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
271 if (first) {
272 amax = a[j];
273 jamax = ja[j];
274 first = FALSE;
275 } else {
276 if (a[j] > amax){
277 amax = a[j];
278 jamax = ja[j];
279 }
280 }
281 }
282 }
283 if (!first){
284 (*matching)[jamax] = i;
285 (*matching)[i] = jamax;
286 (*nmatch)--;
287 }
288 }
289 } else {
290 p = random_permutation(m);
291 for (ii = 0; ii < m; ii++){
292 i = p[ii];
293 if ((*matching)[i] != i) continue;
294 first = TRUE;
295 for (j = ia[i]; j < ia[i+1]; j++){
296 if (i == ja[j]) continue;
297 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
298 if (first) {
299 amax = a[j];
300 jamax = ja[j];
301 first = FALSE;
302 } else {
303 if (a[j] > amax){
304 amax = a[j];
305 jamax = ja[j];
306 }
307 }
308 }
309 }
310 if (!first){
311 (*matching)[jamax] = i;
312 (*matching)[i] = jamax;
313 (*nmatch)--;
314 }
315 }
316 FREE(p);
317 }
318}
319
320
321
322
323
324#define node_degree(i) (ia[(i)+1] - ia[(i)])
325
326static void maximal_independent_edge_set_heavest_edge_pernode_leaves_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
327 int i, ii, j, *ia, *ja, m, n, *p = NULL, q;
328 real *a, amax = 0;
329 int first = TRUE, jamax = 0;
330 int *matched, nz, ncmax = 0, nz0, nzz,k ;
331 enum {UNMATCHED = -2, MATCHED = -1};
332
333 assert(A);
334 assert(SparseMatrix_known_strucural_symmetric(A));
335 ia = A->ia;
336 ja = A->ja;
337 m = A->m;
338 n = A->n;
339 assert(n == m);
340 *cluster = N_GNEW(m,int);
341 *clusterp = N_GNEW((m+1),int);
342 matched = N_GNEW(m,int);
343
344 for (i = 0; i < m; i++) matched[i] = i;
345
346 assert(SparseMatrix_is_symmetric(A, FALSE));
347 assert(A->type == MATRIX_TYPE_REAL);
348
349 *ncluster = 0;
350 (*clusterp)[0] = 0;
351 nz = 0;
352 a = (real*) A->a;
353 if (!randomize){
354 for (i = 0; i < m; i++){
355 if (matched[i] == MATCHED || node_degree(i) != 1) continue;
356 q = ja[ia[i]];
357 assert(matched[q] != MATCHED);
358 matched[q] = MATCHED;
359 (*cluster)[nz++] = q;
360 for (j = ia[q]; j < ia[q+1]; j++){
361 if (q == ja[j]) continue;
362 if (node_degree(ja[j]) == 1){
363 matched[ja[j]] = MATCHED;
364 (*cluster)[nz++] = ja[j];
365 }
366 }
367 ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
368 nz0 = (*clusterp)[*ncluster];
369 if (nz - nz0 <= MAX_CLUSTER_SIZE){
370 (*clusterp)[++(*ncluster)] = nz;
371 } else {
372 (*clusterp)[++(*ncluster)] = ++nz0;
373 nzz = nz0;
374 for (k = nz0; k < nz && nzz < nz; k++){
375 nzz += MAX_CLUSTER_SIZE - 1;
376 nzz = MIN(nz, nzz);
377 (*clusterp)[++(*ncluster)] = nzz;
378 }
379 }
380
381 }
382 #ifdef DEBUG_print
383 if (Verbose)
384 fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax);
385#endif
386 for (i = 0; i < m; i++){
387 first = TRUE;
388 if (matched[i] == MATCHED) continue;
389 for (j = ia[i]; j < ia[i+1]; j++){
390 if (i == ja[j]) continue;
391 if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
392 if (first) {
393 amax = a[j];
394 jamax = ja[j];
395 first = FALSE;
396 } else {
397 if (a[j] > amax){
398 amax = a[j];
399 jamax = ja[j];
400 }
401 }
402 }
403 }
404 if (!first){
405 matched[jamax] = MATCHED;
406 matched[i] = MATCHED;
407 (*cluster)[nz++] = i;
408 (*cluster)[nz++] = jamax;
409 (*clusterp)[++(*ncluster)] = nz;
410 }
411 }
412
413 /* dan yi dian, wu ban */
414 for (i = 0; i < m; i++){
415 if (matched[i] == i){
416 (*cluster)[nz++] = i;
417 (*clusterp)[++(*ncluster)] = nz;
418 }
419 }
420 assert(nz == n);
421
422 } else {
423 p = random_permutation(m);
424 for (ii = 0; ii < m; ii++){
425 i = p[ii];
426 if (matched[i] == MATCHED || node_degree(i) != 1) continue;
427 q = ja[ia[i]];
428 assert(matched[q] != MATCHED);
429 matched[q] = MATCHED;
430 (*cluster)[nz++] = q;
431 for (j = ia[q]; j < ia[q+1]; j++){
432 if (q == ja[j]) continue;
433 if (node_degree(ja[j]) == 1){
434 matched[ja[j]] = MATCHED;
435 (*cluster)[nz++] = ja[j];
436 }
437 }
438 ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
439 nz0 = (*clusterp)[*ncluster];
440 if (nz - nz0 <= MAX_CLUSTER_SIZE){
441 (*clusterp)[++(*ncluster)] = nz;
442 } else {
443 (*clusterp)[++(*ncluster)] = ++nz0;
444 nzz = nz0;
445 for (k = nz0; k < nz && nzz < nz; k++){
446 nzz += MAX_CLUSTER_SIZE - 1;
447 nzz = MIN(nz, nzz);
448 (*clusterp)[++(*ncluster)] = nzz;
449 }
450 }
451 }
452
453 #ifdef DEBUG_print
454 if (Verbose)
455 fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax);
456#endif
457 for (ii = 0; ii < m; ii++){
458 i = p[ii];
459 first = TRUE;
460 if (matched[i] == MATCHED) continue;
461 for (j = ia[i]; j < ia[i+1]; j++){
462 if (i == ja[j]) continue;
463 if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
464 if (first) {
465 amax = a[j];
466 jamax = ja[j];
467 first = FALSE;
468 } else {
469 if (a[j] > amax){
470 amax = a[j];
471 jamax = ja[j];
472 }
473 }
474 }
475 }
476 if (!first){
477 matched[jamax] = MATCHED;
478 matched[i] = MATCHED;
479 (*cluster)[nz++] = i;
480 (*cluster)[nz++] = jamax;
481 (*clusterp)[++(*ncluster)] = nz;
482 }
483 }
484
485 /* dan yi dian, wu ban */
486 for (i = 0; i < m; i++){
487 if (matched[i] == i){
488 (*cluster)[nz++] = i;
489 (*clusterp)[++(*ncluster)] = nz;
490 }
491 }
492
493 FREE(p);
494 }
495
496 FREE(matched);
497}
498
499
500
501static void maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
502 int i, ii, j, *ia, *ja, m, n, *p = NULL;
503 real *a, amax = 0;
504 int first = TRUE, jamax = 0;
505 int *matched, nz, nz0;
506 enum {UNMATCHED = -2, MATCHED = -1};
507 int nsuper, *super = NULL, *superp = NULL;
508
509 assert(A);
510 assert(SparseMatrix_known_strucural_symmetric(A));
511 ia = A->ia;
512 ja = A->ja;
513 m = A->m;
514 n = A->n;
515 assert(n == m);
516 *cluster = N_GNEW(m,int);
517 *clusterp = N_GNEW((m+1),int);
518 matched = N_GNEW(m,int);
519
520 for (i = 0; i < m; i++) matched[i] = i;
521
522 assert(SparseMatrix_is_symmetric(A, FALSE));
523 assert(A->type == MATRIX_TYPE_REAL);
524
525 SparseMatrix_decompose_to_supervariables(A, &nsuper, &super, &superp);
526
527 *ncluster = 0;
528 (*clusterp)[0] = 0;
529 nz = 0;
530 a = (real*) A->a;
531
532 for (i = 0; i < nsuper; i++){
533 if (superp[i+1] - superp[i] <= 1) continue;
534 nz0 = (*clusterp)[*ncluster];
535 for (j = superp[i]; j < superp[i+1]; j++){
536 matched[super[j]] = MATCHED;
537 (*cluster)[nz++] = super[j];
538 if (nz - nz0 >= MAX_CLUSTER_SIZE){
539 (*clusterp)[++(*ncluster)] = nz;
540 nz0 = nz;
541 }
542 }
543 if (nz > nz0) (*clusterp)[++(*ncluster)] = nz;
544 }
545
546 if (!randomize){
547 for (i = 0; i < m; i++){
548 first = TRUE;
549 if (matched[i] == MATCHED) continue;
550 for (j = ia[i]; j < ia[i+1]; j++){
551 if (i == ja[j]) continue;
552 if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
553 if (first) {
554 amax = a[j];
555 jamax = ja[j];
556 first = FALSE;
557 } else {
558 if (a[j] > amax){
559 amax = a[j];
560 jamax = ja[j];
561 }
562 }
563 }
564 }
565 if (!first){
566 matched[jamax] = MATCHED;
567 matched[i] = MATCHED;
568 (*cluster)[nz++] = i;
569 (*cluster)[nz++] = jamax;
570 (*clusterp)[++(*ncluster)] = nz;
571 }
572 }
573
574 /* dan yi dian, wu ban */
575 for (i = 0; i < m; i++){
576 if (matched[i] == i){
577 (*cluster)[nz++] = i;
578 (*clusterp)[++(*ncluster)] = nz;
579 }
580 }
581 assert(nz == n);
582
583 } else {
584 p = random_permutation(m);
585 for (ii = 0; ii < m; ii++){
586 i = p[ii];
587 first = TRUE;
588 if (matched[i] == MATCHED) continue;
589 for (j = ia[i]; j < ia[i+1]; j++){
590 if (i == ja[j]) continue;
591 if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
592 if (first) {
593 amax = a[j];
594 jamax = ja[j];
595 first = FALSE;
596 } else {
597 if (a[j] > amax){
598 amax = a[j];
599 jamax = ja[j];
600 }
601 }
602 }
603 }
604 if (!first){
605 matched[jamax] = MATCHED;
606 matched[i] = MATCHED;
607 (*cluster)[nz++] = i;
608 (*cluster)[nz++] = jamax;
609 (*clusterp)[++(*ncluster)] = nz;
610 }
611 }
612
613 /* dan yi dian, wu ban */
614 for (i = 0; i < m; i++){
615 if (matched[i] == i){
616 (*cluster)[nz++] = i;
617 (*clusterp)[++(*ncluster)] = nz;
618 }
619 }
620 FREE(p);
621
622 }
623
624 FREE(super);
625
626 FREE(superp);
627
628 FREE(matched);
629}
630
631static int scomp(const void *s1, const void *s2){
632 real *ss1, *ss2;
633 ss1 = (real*) s1;
634 ss2 = (real*) s2;
635
636 if ((ss1)[1] > (ss2)[1]){
637 return -1;
638 } else if ((ss1)[1] < (ss2)[1]){
639 return 1;
640 }
641 return 0;
642}
643
644static void maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(SparseMatrix A, int csize,
645 int randomize, int **cluster, int **clusterp, int *ncluster){
646 int i, ii, j, *ia, *ja, m, n, *p = NULL, q, iv;
647 real *a;
648 int *matched, nz, nz0, nzz,k, nv;
649 enum {UNMATCHED = -2, MATCHED = -1};
650 real *vlist;
651
652 assert(A);
653 assert(SparseMatrix_known_strucural_symmetric(A));
654 ia = A->ia;
655 ja = A->ja;
656 m = A->m;
657 n = A->n;
658 assert(n == m);
659 *cluster = N_GNEW(m,int);
660 *clusterp = N_GNEW((m+1),int);
661 matched = N_GNEW(m,int);
662 vlist = N_GNEW(2*m,real);
663
664 for (i = 0; i < m; i++) matched[i] = i;
665
666 assert(SparseMatrix_is_symmetric(A, FALSE));
667 assert(A->type == MATRIX_TYPE_REAL);
668
669 *ncluster = 0;
670 (*clusterp)[0] = 0;
671 nz = 0;
672 a = (real*) A->a;
673
674 p = random_permutation(m);
675 for (ii = 0; ii < m; ii++){
676 i = p[ii];
677 if (matched[i] == MATCHED || node_degree(i) != 1) continue;
678 q = ja[ia[i]];
679 assert(matched[q] != MATCHED);
680 matched[q] = MATCHED;
681 (*cluster)[nz++] = q;
682 for (j = ia[q]; j < ia[q+1]; j++){
683 if (q == ja[j]) continue;
684 if (node_degree(ja[j]) == 1){
685 matched[ja[j]] = MATCHED;
686 (*cluster)[nz++] = ja[j];
687 }
688 }
689 nz0 = (*clusterp)[*ncluster];
690 if (nz - nz0 <= MAX_CLUSTER_SIZE){
691 (*clusterp)[++(*ncluster)] = nz;
692 } else {
693 (*clusterp)[++(*ncluster)] = ++nz0;
694 nzz = nz0;
695 for (k = nz0; k < nz && nzz < nz; k++){
696 nzz += MAX_CLUSTER_SIZE - 1;
697 nzz = MIN(nz, nzz);
698 (*clusterp)[++(*ncluster)] = nzz;
699 }
700 }
701 }
702
703 for (ii = 0; ii < m; ii++){
704 i = p[ii];
705 if (matched[i] == MATCHED) continue;
706 nv = 0;
707 for (j = ia[i]; j < ia[i+1]; j++){
708 if (i == ja[j]) continue;
709 if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
710 vlist[2*nv] = ja[j];
711 vlist[2*nv+1] = a[j];
712 nv++;
713 }
714 }
715 if (nv > 0){
716 qsort(vlist, nv, sizeof(real)*2, scomp);
717 for (j = 0; j < MIN(csize - 1, nv); j++){
718 iv = (int) vlist[2*j];
719 matched[iv] = MATCHED;
720 (*cluster)[nz++] = iv;
721 }
722 matched[i] = MATCHED;
723 (*cluster)[nz++] = i;
724 (*clusterp)[++(*ncluster)] = nz;
725 }
726 }
727
728 /* dan yi dian, wu ban */
729 for (i = 0; i < m; i++){
730 if (matched[i] == i){
731 (*cluster)[nz++] = i;
732 (*clusterp)[++(*ncluster)] = nz;
733 }
734 }
735 FREE(p);
736
737
738 FREE(matched);
739}
740static void maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatrix A, int randomize, int **matching, int *nmatch){
741 int i, ii, j, *ia, *ja, m, n, *p = NULL;
742 real *a, amax = 0;
743 int first = TRUE, jamax = 0;
744
745 assert(A);
746 assert(SparseMatrix_known_strucural_symmetric(A));
747 ia = A->ia;
748 ja = A->ja;
749 m = A->m;
750 n = A->n;
751 assert(n == m);
752 *matching = N_GNEW(m,int);
753 for (i = 0; i < m; i++) (*matching)[i] = i;
754 *nmatch = n;
755
756 assert(SparseMatrix_is_symmetric(A, FALSE));
757 assert(A->type == MATRIX_TYPE_REAL);
758
759 a = (real*) A->a;
760 if (!randomize){
761 for (i = 0; i < m; i++){
762 first = TRUE;
763 for (j = ia[i]; j < ia[i+1]; j++){
764 if (i == ja[j]) continue;
765 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
766 if (first) {
767 amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
768 jamax = ja[j];
769 first = FALSE;
770 } else {
771 if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){
772 amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
773 jamax = ja[j];
774 }
775 }
776 }
777 }
778 if (!first){
779 (*matching)[jamax] = i;
780 (*matching)[i] = jamax;
781 (*nmatch)--;
782 }
783 }
784 } else {
785 p = random_permutation(m);
786 for (ii = 0; ii < m; ii++){
787 i = p[ii];
788 if ((*matching)[i] != i) continue;
789 first = TRUE;
790 for (j = ia[i]; j < ia[i+1]; j++){
791 if (i == ja[j]) continue;
792 if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
793 if (first) {
794 amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
795 jamax = ja[j];
796 first = FALSE;
797 } else {
798 if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){
799 amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
800 jamax = ja[j];
801 }
802 }
803 }
804 }
805 if (!first){
806 (*matching)[jamax] = i;
807 (*matching)[i] = jamax;
808 (*nmatch)--;
809 }
810 }
811 FREE(p);
812 }
813}
814
815SparseMatrix DistanceMatrix_restrict_cluster(int ncluster, int *clusterp, int *cluster, SparseMatrix P, SparseMatrix R, SparseMatrix D){
816#if 0
817 /* this construct a distance matrix of a coarse graph, for a coarsen give by merging all nodes in each cluster */
818 SparseMatrix cD = NULL;
819 int i, j, nzc;
820 int **irn, **jcn;
821 real **val;
822 int n = D->m;
823 int *assignment = NULL;
824 int nz;
825 int *id = D->ia, jd = D->ja;
826 int *mask = NULL;
827 int *nnodes, *mask;
828 real *d = NULL;
829
830
831 assert(D->m == D->n);
832 if (!D) return NULL;
833 if (D->a && D->type == MATRIX_TYPE_REAL) d = (real*) D->val;
834
835 irn = N_GNEW(ncluster,int*);
836 jcn = N_GNEW(ncluster,int*);
837 val = N_GNEW(ncluster,real*);
838 assignment = N_GNEW(n,int);
839 nz = N_GNEW(ncluster,int);
840 mask = N_GNEW(n,int);
841 nnodes = N_GNEW(ncluster,int);
842
843
844 /* find ncluster-subgrahs induced by the ncluster -clusters, find the diameter of each,
845 then use the radius as the distance from the supernode to the rest of the "world"
846 */
847 for (i = 0; i < ncluster; i++) nz[i] = 0;
848 for (i = 0; i < ncluster; i++){
849 for (j = clusterp[i]; j < clusterp[i+1]; j++){
850 assert(clusterp[i+1] > clusterp[i]);
851 assignment[cluster[j]] = i;
852 }
853 }
854
855 for (i = 0; i < n; i++){/* figure out how many entries per submatrix */
856 ic = asignment[i];
857 for (j = id[i]; j < id[i+1]; j++){
858 if (i != jd[j] && ic == assignment[jd[j]]) {
859 nz[ic]++;
860 }
861 }
862 }
863 for (i = 0; i < ncluster; i++) {
864 irn[i] = N_GNEW(nz[i],int);
865 jcn[i] = N_GNEW(nz[i],int);
866 val[i] = N_GNEW(nz[i],int);
867 val[i] = NULL;
868 }
869
870
871 for (i = 0; i < ncluster; i++) nz[i] = 0;/* get subgraphs */
872 for (i = 0; i < n; i++) mask[i] = -1;
873 for (i = 0; i < ncluster; i++) nnodes[i] = -1;
874 for (i = 0; i < n; i++){
875 ic = asignment[i];
876 ii = mask[i];
877 if (ii < 0){
878 mask[i] = ii = nnodes[ic];
879 nnodes[ic]++;
880 }
881 for (j = id[i]; j < id[i+1]; j++){
882 jc = assignment[jd[j]];
883 if (i != jd[j] && ic == jc) {
884 jj = mask[jd[j]];
885 if (jj < 0){
886 mask[jd[j]] = jj = nnodes[jc];
887 nnodes[jc]++;
888 }
889 irn[ic][nz[ic]] = ii;
890 jcn[ic][nz[ic]] = jj;
891 if (d) val[ic][nz[ic]] = d[j];
892 }
893 }
894 }
895
896 for (i = 0; i < ncluster; i++){/* form subgraphs */
897 SparseMatrix A;
898 A = SparseMatrix_from_coordinate_arrays(nz[nz[i]], nnodes[i], nnodes[i], irn[i], jcn[i], (void*) val[i], MATRIX_TYPE_REAL);
899
900 SparseMatrix_delete(A);
901 }
902
903
904 for (i = 0; i < ncluster; i++){
905 for (j = clusterp[i]; j < clusterp[i+1]; j++){
906 assert(clusterp[i+1] > clusterp[i]);
907 irn[nzc] = cluster[j];
908 jcn[nzc] = i;
909 val[nzc++] = 1.;
910 }
911 }
912 assert(nzc == n);
913 cD = SparseMatrix_multiply3(R, D, P);
914
915 SparseMatrix_set_symmetric(cD);
916 SparseMatrix_set_pattern_symmetric(cD);
917 cD = SparseMatrix_remove_diagonal(cD);
918
919 FREE(nz);
920 FREE(assignment);
921 for (i = 0; i < ncluster; i++){
922 FREE(irn[i]);
923 FREE(jcn[i]);
924 FREE(val[i]);
925 }
926 FREE(irn); FREE(jcn); FREE(val);
927 FREE(mask);
928 FREE(nnodes);
929
930 return cD;
931#endif
932 return NULL;
933}
934
935SparseMatrix DistanceMatrix_restrict_matching(int *matching, SparseMatrix D){
936 if (!D) return NULL;
937 assert(0);/* not yet implemented! */
938 return NULL;
939}
940
941SparseMatrix DistanceMatrix_restrict_filtering(int *mask, int is_C, int is_F, SparseMatrix D){
942 /* max independent vtx set based coarsening. Coarsen nodes has mask >= is_C. Fine nodes == is_F. */
943 if (!D) return NULL;
944 assert(0);/* not yet implemented! */
945 return NULL;
946}
947
948static void Multilevel_coarsen_internal(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD,
949 real *node_wgt, real **cnode_wgt,
950 SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){
951 int *matching = NULL, nmatch = 0, nc, nzc, n, i;
952 int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL;
953 real *val = NULL;
954 SparseMatrix B = NULL;
955 int *vset = NULL, nvset, ncov, j;
956 int *cluster=NULL, *clusterp=NULL, ncluster;
957
958 assert(A->m == A->n);
959 *cA = NULL;
960 *cD = NULL;
961 *P = NULL;
962 *R = NULL;
963 n = A->m;
964
965 *coarsen_scheme_used = ctrl->coarsen_scheme;
966
967 switch (ctrl->coarsen_scheme){
968 case COARSEN_HYBRID:
969#ifdef DEBUG_PRINT
970 if (Verbose)
971 fprintf(stderr, "hybrid scheme, try COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST first\n");
972#endif
973 *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST;
974 Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
975
976 if (!(*cA)) {
977#ifdef DEBUG_PRINT
978 if (Verbose)
979 fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST\n");
980#endif
981 *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST;
982 Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
983 }
984
985 if (!(*cA)) {
986#ifdef DEBUG_PRINT
987 if (Verbose)
988 fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST\n");
989#endif
990 *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
991 Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
992 }
993
994 if (!(*cA)) {
995#ifdef DEBUG_PRINT
996 if (Verbose)
997 fprintf(stderr, "switching to COARSEN_INDEPENDENT_VERTEX_SET\n");
998#endif
999 *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET;
1000 Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
1001 }
1002
1003
1004 if (!(*cA)) {
1005#ifdef DEBUG_PRINT
1006 if (Verbose)
1007 fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE\n");
1008#endif
1009 *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE;
1010 Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
1011 }
1012 ctrl->coarsen_scheme = COARSEN_HYBRID;
1013 break;
1014 case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST:
1015 case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST:
1016 case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST:
1017 if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST) {
1018 maximal_independent_edge_set_heavest_edge_pernode_leaves_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster);
1019 } else if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST) {
1020 maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster);
1021 } else {
1022 maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(A, 4, ctrl->randomize, &cluster, &clusterp, &ncluster);
1023 }
1024 assert(ncluster <= n);
1025 nc = ncluster;
1026 if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc == n || nc < ctrl->minsize) {
1027#ifdef DEBUG_PRINT
1028 if (Verbose)
1029 fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
1030#endif
1031 goto RETURN;
1032 }
1033 irn = N_GNEW(n,int);
1034 jcn = N_GNEW(n,int);
1035 val = N_GNEW(n,real);
1036 nzc = 0;
1037 for (i = 0; i < ncluster; i++){
1038 for (j = clusterp[i]; j < clusterp[i+1]; j++){
1039 assert(clusterp[i+1] > clusterp[i]);
1040 irn[nzc] = cluster[j];
1041 jcn[nzc] = i;
1042 val[nzc++] = 1.;
1043 }
1044 }
1045 assert(nzc == n);
1046 *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
1047 *R = SparseMatrix_transpose(*P);
1048
1049 *cD = DistanceMatrix_restrict_cluster(ncluster, clusterp, cluster, *P, *R, D);
1050
1051 *cA = SparseMatrix_multiply3(*R, A, *P);
1052
1053 /*
1054 B = SparseMatrix_multiply(*R, A);
1055 if (!B) goto RETURN;
1056 *cA = SparseMatrix_multiply(B, *P);
1057 */
1058 if (!*cA) goto RETURN;
1059
1060 SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
1061 *R = SparseMatrix_divide_row_by_degree(*R);
1062 SparseMatrix_set_symmetric(*cA);
1063 SparseMatrix_set_pattern_symmetric(*cA);
1064 *cA = SparseMatrix_remove_diagonal(*cA);
1065
1066
1067
1068 break;
1069 case COARSEN_INDEPENDENT_EDGE_SET:
1070 maximal_independent_edge_set(A, ctrl->randomize, &matching, &nmatch);
1071 case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE:
1072 if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE)
1073 maximal_independent_edge_set_heavest_edge_pernode(A, ctrl->randomize, &matching, &nmatch);
1074 case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED:
1075 if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED)
1076 maximal_independent_edge_set_heavest_edge_pernode_scaled(A, ctrl->randomize, &matching, &nmatch);
1077 nc = nmatch;
1078 if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc == n || nc < ctrl->minsize) {
1079#ifdef DEBUG_PRINT
1080 if (Verbose)
1081 fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
1082#endif
1083 goto RETURN;
1084 }
1085 irn = N_GNEW(n,int);
1086 jcn = N_GNEW(n,int);
1087 val = N_GNEW(n,real);
1088 nzc = 0; nc = 0;
1089 for (i = 0; i < n; i++){
1090 if (matching[i] >= 0){
1091 if (matching[i] == i){
1092 irn[nzc] = i;
1093 jcn[nzc] = nc;
1094 val[nzc++] = 1.;
1095 } else {
1096 irn[nzc] = i;
1097 jcn[nzc] = nc;
1098 val[nzc++] = 1;
1099 irn[nzc] = matching[i];
1100 jcn[nzc] = nc;
1101 val[nzc++] = 1;
1102 matching[matching[i]] = -1;
1103 }
1104 nc++;
1105 matching[i] = -1;
1106 }
1107 }
1108 assert(nc == nmatch);
1109 assert(nzc == n);
1110 *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
1111 *R = SparseMatrix_transpose(*P);
1112 *cA = SparseMatrix_multiply3(*R, A, *P);
1113 /*
1114 B = SparseMatrix_multiply(*R, A);
1115 if (!B) goto RETURN;
1116 *cA = SparseMatrix_multiply(B, *P);
1117 */
1118 if (!*cA) goto RETURN;
1119 SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
1120 *R = SparseMatrix_divide_row_by_degree(*R);
1121 SparseMatrix_set_symmetric(*cA);
1122 SparseMatrix_set_pattern_symmetric(*cA);
1123 *cA = SparseMatrix_remove_diagonal(*cA);
1124
1125
1126 *cD = DistanceMatrix_restrict_matching(matching, D);
1127 *cD=NULL;
1128
1129 break;
1130 case COARSEN_INDEPENDENT_VERTEX_SET:
1131 case COARSEN_INDEPENDENT_VERTEX_SET_RS:
1132 if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_VERTEX_SET){
1133 maximal_independent_vertex_set(A, ctrl->randomize, &vset, &nvset, &nzc);
1134 } else {
1135 maximal_independent_vertex_set_RS(A, ctrl->randomize, &vset, &nvset, &nzc);
1136 }
1137 ia = A->ia;
1138 ja = A->ja;
1139 nc = nvset;
1140 if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc == n || nc < ctrl->minsize) {
1141#ifdef DEBUG_PRINT
1142 if (Verbose)
1143 fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
1144#endif
1145 goto RETURN;
1146 }
1147 irn = N_GNEW(nzc,int);
1148 jcn = N_GNEW(nzc,int);
1149 val = N_GNEW(nzc,real);
1150 nzc = 0;
1151 for (i = 0; i < n; i++){
1152 if (vset[i] == MAX_IND_VTX_SET_F){
1153 ncov = 0;
1154 for (j = ia[i]; j < ia[i+1]; j++){
1155 if (vset[ja[j]] >= MAX_IND_VTX_SET_C){
1156 ncov++;
1157 }
1158 }
1159 assert(ncov > 0);
1160 for (j = ia[i]; j < ia[i+1]; j++){
1161 if (vset[ja[j]] >= MAX_IND_VTX_SET_C){
1162 irn[nzc] = i;
1163 jcn[nzc] = vset[ja[j]];
1164 val[nzc++] = 1./(double) ncov;
1165 }
1166 }
1167 } else {
1168 assert(vset[i] >= MAX_IND_VTX_SET_C);
1169 irn[nzc] = i;
1170 jcn[nzc] = vset[i];
1171 val[nzc++] = 1.;
1172 }
1173 }
1174
1175 *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
1176 *R = SparseMatrix_transpose(*P);
1177 *cA = SparseMatrix_multiply3(*R, A, *P);
1178 if (!*cA) goto RETURN;
1179 SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
1180 SparseMatrix_set_symmetric(*cA);
1181 SparseMatrix_set_pattern_symmetric(*cA);
1182 *cA = SparseMatrix_remove_diagonal(*cA);
1183
1184 *cD = DistanceMatrix_restrict_filtering(vset, MAX_IND_VTX_SET_C, MAX_IND_VTX_SET_F, D);
1185 break;
1186 default:
1187 goto RETURN;
1188 }
1189 RETURN:
1190 if (matching) FREE(matching);
1191 if (vset) FREE(vset);
1192 if (irn) FREE(irn);
1193 if (jcn) FREE(jcn);
1194 if (val) FREE(val);
1195 if (B) SparseMatrix_delete(B);
1196
1197 if(cluster) FREE(cluster);
1198 if(clusterp) FREE(clusterp);
1199}
1200
1201void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, real *node_wgt, real **cnode_wgt,
1202 SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){
1203 SparseMatrix cA0 = A, cD0 = NULL, P0 = NULL, R0 = NULL, M;
1204 real *cnode_wgt0 = NULL;
1205 int nc = 0, n;
1206
1207 *P = NULL; *R = NULL; *cA = NULL; *cnode_wgt = NULL, *cD = NULL;
1208
1209 n = A->n;
1210
1211 do {/* this loop force a sufficient reduction */
1212 node_wgt = cnode_wgt0;
1213 Multilevel_coarsen_internal(A, &cA0, D, &cD0, node_wgt, &cnode_wgt0, &P0, &R0, ctrl, coarsen_scheme_used);
1214 if (!cA0) return;
1215 nc = cA0->n;
1216#ifdef DEBUG_PRINT
1217 if (Verbose) fprintf(stderr,"nc=%d n = %d\n",nc,n);
1218#endif
1219 if (*P){
1220 assert(*R);
1221 M = SparseMatrix_multiply(*P, P0);
1222 SparseMatrix_delete(*P);
1223 SparseMatrix_delete(P0);
1224 *P = M;
1225 M = SparseMatrix_multiply(R0, *R);
1226 SparseMatrix_delete(*R);
1227 SparseMatrix_delete(R0);
1228 *R = M;
1229 } else {
1230 *P = P0;
1231 *R = R0;
1232 }
1233
1234 if (*cA) SparseMatrix_delete(*cA);
1235 *cA = cA0;
1236 if (*cD) SparseMatrix_delete(*cD);
1237 *cD = cD0;
1238
1239 if (*cnode_wgt) FREE(*cnode_wgt);
1240 *cnode_wgt = cnode_wgt0;
1241 A = cA0;
1242 D = cD0;
1243 node_wgt = cnode_wgt0;
1244 cnode_wgt0 = NULL;
1245 } while (nc > ctrl->min_coarsen_factor*n && ctrl->coarsen_mode == COARSEN_MODE_FORCEFUL);
1246
1247}
1248
1249void print_padding(int n){
1250 int i;
1251 for (i = 0; i < n; i++) fputs (" ", stderr);
1252}
1253static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl){
1254 Multilevel cgrid;
1255 int coarsen_scheme_used;
1256 real *cnode_weights = NULL;
1257 SparseMatrix P, R, A, cA, D, cD;
1258
1259#ifdef DEBUG_PRINT
1260 if (Verbose) {
1261 print_padding(grid->level);
1262 fprintf(stderr, "level -- %d, n = %d, nz = %d nz/n = %f\n", grid->level, grid->n, grid->A->nz, grid->A->nz/(double) grid->n);
1263 }
1264#endif
1265 A = grid->A;
1266 D = grid->D;
1267 if (grid->level >= ctrl->maxlevel - 1) {
1268#ifdef DEBUG_PRINT
1269 if (Verbose) {
1270 print_padding(grid->level);
1271 fprintf(stderr, " maxlevel reached, coarsening stops\n");
1272 }
1273#endif
1274 return grid;
1275 }
1276 Multilevel_coarsen(A, &cA, D, &cD, grid->node_weights, &cnode_weights, &P, &R, ctrl, &coarsen_scheme_used);
1277 if (!cA) return grid;
1278
1279 cgrid = Multilevel_init(cA, cD, cnode_weights);
1280 grid->next = cgrid;
1281 cgrid->coarsen_scheme_used = coarsen_scheme_used;
1282 cgrid->level = grid->level + 1;
1283 cgrid->n = cA->m;
1284 cgrid->A = cA;
1285 cgrid->D = cD;
1286 cgrid->P = P;
1287 grid->R = R;
1288 cgrid->prev = grid;
1289 cgrid = Multilevel_establish(cgrid, ctrl);
1290 return grid;
1291
1292}
1293
1294Multilevel Multilevel_new(SparseMatrix A0, SparseMatrix D0, real *node_weights, Multilevel_control ctrl){
1295 /* A: the weighting matrix. D: the distance matrix, could be NULL. If not null, the two matrices must have the same sparsity pattern */
1296 Multilevel grid;
1297 SparseMatrix A = A0, D = D0;
1298
1299 if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
1300 A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
1301 }
1302 if (D && (!SparseMatrix_is_symmetric(D, FALSE) || D->type != MATRIX_TYPE_REAL)){
1303 D = SparseMatrix_symmetrize_nodiag(D, FALSE);
1304 }
1305 grid = Multilevel_init(A, D, node_weights);
1306 grid = Multilevel_establish(grid, ctrl);
1307 if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */
1308 return grid;
1309}
1310
1311
1312Multilevel Multilevel_get_coarsest(Multilevel grid){
1313 while (grid->next){
1314 grid = grid->next;
1315 }
1316 return grid;
1317}
1318
1319