1 /*********************************************************/
2 /* TAUCS */
3 /* Author: Doron Chen and Sivan Toledo */
4 /* File : taucs_vaidya.c */
5 /* Description: constructs Vaidya's preconditioners */
6 /*********************************************************/
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <math.h>
12 #include <assert.h>
13 #include "taucs.h"
14
15 /*long int random(void); omer*/
16
17 /*********************************************************/
18 /* */
19 /*********************************************************/
20
21 #ifdef TAUCS_CORE_DOUBLE
22
23 typedef unsigned char byte;
24 /*typedef int byte;*/
25
26 #define Do(i,n) for ((i)=0;(i)<(n);(i)++)
27
28 typedef struct {
29 int n;
30 int nent;
31 int max_size;
32 int *ivec1;
33 int *ivec2;
34 double *dvec;
35 } sym_matrix;
36
37 typedef struct {
38 int i;
39 int j;
40 double v;
41 } wedge; /* weighted edge */
42
43 typedef struct {
44 int n;
45 int nent;
46 int max_size;
47 wedge* edges;
48 } graph;
49
50 /************** UNION FIND ***********/
51
52 static char *label = NULL;
53 static int *p = NULL;
54 static int *rank = NULL;
55
56 static
unionfind_init(int size)57 int unionfind_init(int size)
58 {
59 int i;
60 p = (int *)taucs_malloc(size * sizeof(int));
61 rank = (int *)taucs_malloc(size * sizeof(int));
62 label = (char *)taucs_malloc(size * sizeof(char));
63 if (!p || !rank || !label) {
64 taucs_free(p);
65 taucs_free(rank);
66 taucs_free(label);
67 return -1;
68 }
69
70 Do(i,size)
71 {
72 p[i] = i;
73 rank[i] = 0;
74 label[i] = 0;
75 }
76
77 return 0;
78 }
79
80 static void
unionfind_free(void)81 unionfind_free (void)
82 {
83 taucs_free(p);
84 taucs_free(rank);
85 taucs_free(label);
86 }
87
88 static
Union(int a,int b,int x,int y,int l)89 int Union(int a,int b,int x,int y,int l) /* unite a's and b's trees, whose roots are x and y. returns the root
90 of the united tree */
91 {
92 if (rank[x] > rank[y])
93 {
94 p[y] = x;
95 label[y] = label[a] ^ label[b] ^ l;
96 return(x);
97 }
98 else
99 {
100 p[x] = y;
101 label[x] = label[a] ^ label[b] ^ l;
102 if (rank[x] == rank[y])
103 rank[y]++;
104 return(y);
105 }
106 }
107
108 static
find_set(int x)109 int find_set(int x)
110 {
111 int tmp;
112
113 if (x != p[x])
114 {
115 tmp = find_set(p[x]);
116 label[x] ^= label[p[x]];
117 p[x] = tmp;
118 }
119 return(p[x]);
120 }
121
122 /*********** HEAP OPERATIONS ************/
123
124 typedef struct hea {
125 int heap_size;
126 int alloc_size;
127 int *edges;
128 double *key;
129 } heap;
130
131 #define INF 100000000.0
132
133 #define Parent(i) ((((i)+1)/2) - 1)
134 #define Left(i) ((((i)+1)*2) - 1)
135 #define Right(i) ((((i)+1)*2 + 1) - 1)
136
137 static
exchange(heap A,int a,int b)138 void exchange(heap A,int a,int b)
139 {
140 int tmp1;
141 double tmp2;
142
143 tmp1 = A.edges[a];
144 A.edges[a] = A.edges[b];
145 A.edges[b] = tmp1;
146
147 tmp2 = A.key[a];
148 A.key[a] = A.key[b];
149 A.key[b] = tmp2;
150
151 }
152
153 static
Heapify(heap A,int i)154 void Heapify(heap A,int i)
155 {
156 int l,r,smallest;
157
158 l = Left(i);
159 r = Right(i);
160
161 if ((l < A.heap_size) && (A.key[l] < A.key[i]))
162 smallest = l;
163 else
164 smallest = i;
165
166 if ((r < A.heap_size) && (A.key[r] < A.key[smallest]))
167 smallest = r;
168
169 if (smallest != i)
170 {
171 exchange(A,i,smallest);
172 Heapify(A,smallest);
173 }
174 }
175
176 #if 0
177 static
178 int build_heap(int size,heap *h,graph *A)
179 {
180 int i,k=0;
181
182 h->heap_size = 0;
183 h->edges = (int *)taucs_malloc(size * sizeof(int));
184 h->key = (double *)taucs_malloc(size * sizeof(double));
185 if (!(h->edges) || !(h->key)) {
186 taucs_free(h->edges);
187 taucs_free(h->key);
188 return -1;
189 }
190
191 Do(i,size)
192 if ((A->edges)[i].i != (A->edges)[i].j)
193 {
194 h->edges[k] = i;
195 h->key[k] = fabs((A->edges)[i].v);
196 k++;
197 }
198
199 h->heap_size = k;
200 size = k;
201
202 /*
203 for(i=(size/2)-1;i>=0;i--)
204 Heapify(*h,i);
205 */
206
207 return 0;
208 }
209 #endif
210
211 static
free_heap(heap h)212 void free_heap(heap h)
213 {
214 taucs_free(h.edges);
215 taucs_free(h.key);
216 }
217
218
partition(heap h,int p,int r)219 static int partition(heap h, int p, int r)
220 {
221 int pivot;
222 double x;
223 int i,j;
224
225 if (r-p < 16) pivot = 0;
226 /* Sivan: chnaged random() to rand() to remove warning (random is not ansi C) */
227 else if ((r - p) < 8) pivot = p + (rand() % (r-p+1));
228 else {
229 int c[3]; /* candidates */
230 int t;
231 c[0] = p + (rand() % (r-p+1));
232 c[1] = p + (rand() % (r-p+1));
233 c[2] = p + (rand() % (r-p+1));
234
235 if (h.key[c[1]] < h.key[c[0]]) {t=c[0]; c[0]=c[1]; c[1]=t;}
236 if (h.key[c[2]] < h.key[c[1]]) {t=c[1]; c[1]=c[2]; c[2]=t;}
237 if (h.key[c[1]] < h.key[c[0]]) {t=c[0]; c[0]=c[1]; c[1]=t;}
238
239 pivot = c[1];
240 }
241
242 x = h.key[pivot];
243 /*
244 i = p-1;
245 j = r+1;
246
247 while (1) {
248 do {
249 j--;
250 } while ( h.key[j] > x );
251
252 do {
253 i++;
254 } while ( h.key[i] < x );
255
256 if (i < j)
257 exchange(h,i,j);
258 else
259 return j;
260 }
261 */
262 i = p-1;
263 j = r+1;
264
265 while (1) {
266 for (j--; h.key[j] > x; j--);
267 for (i++; h.key[i] < x; i++);
268
269 if (i < j)
270 exchange(h,i,j);
271 else
272 return j;
273 }
274 }
275
276 #if 0
277 static void insertion_sort(heap h, int p, int r)
278 {
279 int i,j;
280
281 for (j=p+1; j<r; j++) {
282 double key = h.key[j];
283 int edge = h.edges[j];
284
285 for (i=j-1; i>=p && h.key[i] > key; i--) {
286 h.key[i+1] = h.key[i];
287 h.edges[i+1] = h.edges[i];
288 }
289
290 h.key[i+1] = key;
291 h.edges[i+1] = edge;
292 }
293 }
294 #endif /* 0, we don't need insertion sort, heap sort */
295
296 static
heapify_offset(heap A,int p,int r,int i)297 void heapify_offset(heap A,int p,int r,int i)
298 {
299 int L,R,smallest;
300 int size = r-p+1;
301
302 L = Left(i);
303 R = Right(i);
304
305 if ((L < size) && (A.key[p+L] < A.key[p+i]))
306 smallest = L;
307 else
308 smallest = i;
309
310 if ((R < size) && (A.key[p+R] < A.key[p+smallest]))
311 smallest = R;
312
313 if (smallest != i)
314 {
315 exchange(A,p+i,p+smallest);
316 heapify_offset(A,p,r,smallest);
317 }
318 }
319
320
heapsort_sort(heap h,int p,int r)321 static void heapsort_sort(heap h, int p, int r)
322 {
323 int size = r-p+1;
324 int i;
325
326 for(i=(size/2)-1;i>=0;i--)
327 heapify_offset(h,p,r,i);
328
329 for(i=size-1;i>=1;i--)
330 {
331 exchange(h,0,i);
332 heapify_offset(h,0,i,i);
333 }
334 }
335
quick_sort(heap h,int p,int r)336 static void quick_sort(heap h, int p, int r)
337 {
338 int q;
339 if (p >= r) return;
340 if (r - p < 100) {
341 /*insertion_sort(h,p,r);*/
342 heapsort_sort(h,p,r);
343 return;
344 }
345 q = partition(h,p,r);
346 quick_sort(h,p,q);
347 quick_sort(h,q+1,r);
348 }
349
350 #if 0
351 static
352 int heap_sort(int size,heap *h,graph *A)
353 {
354 int i;
355 /*double *dvec;*/
356
357 if (build_heap(size,h,A) == -1)
358 return -1;
359 size = h->heap_size;
360
361 #define noQSORT
362 #ifdef QSORT
363 quick_sort(*h,0,size-1);
364 #else
365
366 for(i=(size/2)-1;i>=0;i--)
367 Heapify(*h,i);
368
369
370 for(i=size-1;i>=1;i--)
371 {
372 exchange(*h,0,i);
373 h->heap_size--;
374 Heapify(*h,0);
375 }
376 #endif
377
378 return(size); /* cannot be -1, so -1 is an error */
379 }
380 #endif /* 0, no heap_sort */
381
pqueue_fill(heap * h,graph * G)382 static int pqueue_fill(heap* h, graph* G)
383 {
384 int i,size;
385
386 size=0;
387 Do(i,G->nent) {
388 if ((G->edges)[i].i != (G->edges)[i].j) {
389 assert(size <= h->alloc_size);
390 h->edges[size] = i;
391 h->key[size] = fabs((G->edges)[i].v);
392 size++;
393 }
394 }
395
396 h->heap_size = size;
397
398 #define noQSORT
399 #ifdef QSORT
400 quick_sort(*h,0,size-1);
401 #else
402
403 for(i=(size/2)-1;i>=0;i--)
404 Heapify(*h,i);
405
406 for(i=size-1;i>=1;i--) {
407 exchange(*h,0,i);
408 h->heap_size--;
409 Heapify(*h,0);
410 }
411 #endif
412
413 h->heap_size = size;
414
415 return size;
416 }
417
pqueue_create(heap * h,int size)418 static int pqueue_create(heap* h, int size)
419 {
420 h->heap_size = 0;
421 h->alloc_size = size;
422 h->edges = (int *)taucs_malloc(size * sizeof(int));
423 h->key = (double *)taucs_malloc(size * sizeof(double));
424 if (!(h->edges) || !(h->key)) {
425 taucs_free(h->edges);
426 taucs_free(h->key);
427 return -1;
428 }
429 return 0;
430 }
431
432 /***************************************************/
433 #ifdef GRAPHSORT
434
435 #define Parent(i) ((((i)+1)/2) - 1)
436 #define Left(i) ((((i)+1)*2) - 1)
437 #define Right(i) ((((i)+1)*2 + 1) - 1)
438
439 static
new_exchange(wedge * e,int i,int j)440 void new_exchange(wedge* e,int i,int j)
441 {
442 wedge t;
443 t = e[i];
444 e[i] = e[j];
445 e[j] = t;
446 /*
447 int tmp1;
448 double tmp2;
449
450 tmp1 = A.edges[a];
451 A.edges[a] = A.edges[b];
452 A.edges[b] = tmp1;
453
454 tmp2 = A.key[a];
455 A.key[a] = A.key[b];
456 A.key[b] = tmp2;
457 */
458 }
459
460 static
heapify(wedge * e,int n,int i)461 void heapify(wedge* e,int n,int i)
462 {
463 int l,r,smallest;
464
465 l = Left(i);
466 r = Right(i);
467
468 if ((l < n) && (fabs(e[l].v) < fabs(e[i].v)))
469 smallest = l;
470 else
471 smallest = i;
472
473 if ((r < n) && (fabs(e[r].v) < fabs(e[smallest].v)))
474 smallest = r;
475
476 if (smallest != i)
477 {
478 new_exchange(e,i,smallest);
479 heapify(e,n,smallest);
480 }
481 }
482
483 static
new_heap_sort(wedge * e,int n)484 int new_heap_sort(wedge* e, int n)
485 {
486 int i;
487
488 for (i=(n/2)-1; i>=0; i--)
489 heapify(e,n,i);
490
491 for(i=n-1; i>=1; i--) {
492 new_exchange(e,0,i);
493 n--;
494 heapify(e,n,0);
495 }
496 }
497
498
wedge_compare(const void * ve1,const void * ve2)499 int wedge_compare(const void* ve1, const void* ve2)
500 {
501 wedge* e1 = (wedge*) ve1;
502 wedge* e2 = (wedge*) ve2;
503
504 double k1, k2;
505
506 /*k1 = fabs(e1->v);*/
507 /*k2 = fabs(e2->v);*/
508
509 k1 = fabs(e1->v);
510 k2 = fabs(e2->v);
511
512 if (k1 < k2) return -1;
513 if (k1 > k2) return 1;
514 return 0;
515 }
516
517
518
insertion_sort(wedge * e,int n)519 static insertion_sort(wedge* e, int n)
520 {
521 int i,j;
522
523 for (j=1; j<n; j++) {
524 double key = fabs(e[j].v);
525 wedge ej = e[j];
526 for (i=j-1; i>=0 && fabs(e[i].j) > key; i--) {
527 /*e[i+1] = e[i];*/
528
529 e[i+1].i = e[i].i;
530 e[i+1].j = e[i].j;
531 e[i+1].v = e[i].v;
532 }
533 e[i+1].i = ej.i;
534 e[i+1].j = ej.j;
535 e[i+1].v = ej.v;
536 /* e[i+1] = ej;*/
537 }
538 }
539
partition(wedge * e,int n)540 static int partition(wedge* e, int n)
541 {
542 int pivot = (rand() % n);
543 double x = fabs(e[pivot].v);
544 int i,j;
545
546 i = -1;
547 j = n;
548
549 while (1) {
550 do {
551 j--;
552 } while ( fabs(e[j].v) > x );
553
554 do {
555 i++;
556 } while ( fabs(e[i].v) < x );
557
558 if (i < j) {
559 /*
560 wedge t;
561 t = e[i];
562 e[i] = e[j];
563 e[j] = t;
564 */
565
566 int ti,tj; double tv;
567 ti = e[i].i;
568 tj = e[i].j;
569 tv = e[i].v;
570 e[i].i = e[j].i;
571 e[i].j = e[j].j;
572 e[i].v = e[j].v;
573 e[j].i = ti;
574 e[j].j = tj;
575 e[j].v = tv;
576 } else return j;
577 }
578 }
579
quick_sort(wedge * e,int n)580 static quick_sort(wedge* e, int n)
581 {
582 int q;
583 if (n <= 1) return;
584 if (n < 32) {
585 insertion_sort(e,n);
586 return;
587 }
588 q = partition(e,n);
589 quick_sort(e ,q+1);
590 quick_sort(e+q+1,n-q-1);
591 }
592
593 static int
graph_sort(graph * G)594 graph_sort(graph* G)
595 {
596 /*
597 qsort(G->edges,
598 G->nent,
599 sizeof(wedge),
600 wedge_compare);
601 */
602 quick_sort(G->edges,G->nent);
603 }
604
605 #endif /* GRAPHSORT */
606 /************ VAIDYA'S PRECONDITIONERS *************/
607
608 #define swap(a,b) {int TMP; TMP = a; a = b; b = TMP;}
609 #define EPSILON 0.00000001
610
611 /************ GRAPHS *************/
612
613 static
construct_graph(int size)614 graph* construct_graph(int size)
615 {
616 graph *out;
617
618 out = (graph *)taucs_malloc(sizeof(graph));
619 if (!out) return NULL;
620
621 out->edges = (wedge*) taucs_malloc(size*sizeof(wedge));
622 if (!(out->edges)) {
623 taucs_free(out);
624 return NULL;
625 }
626
627 out->max_size = size;
628
629 return out;
630 }
631
632 static
free_graph(graph * a)633 void free_graph(graph *a)
634 {
635 if(a)
636 {
637 taucs_free(a->edges);
638 taucs_free(a);
639 }
640 }
641
free_ccs_matrix(taucs_ccs_matrix * a)642 void free_ccs_matrix(taucs_ccs_matrix *a)
643 {
644 if (a)
645 {
646 taucs_free(a->rowind);
647 taucs_free(a->colptr);
648 taucs_free(a->values.d/*taucs_values*/);
649 taucs_free(a);
650 }
651 }
652
653 static
construct_ccs_matrix(int nent,int n)654 taucs_ccs_matrix* construct_ccs_matrix(int nent,int n)
655 {
656 taucs_ccs_matrix *out;
657
658 out = (taucs_ccs_matrix *)taucs_malloc(sizeof(taucs_ccs_matrix));
659 if (!out) return NULL;
660 out->colptr = (int *)taucs_malloc((n+1)*sizeof(int));
661 out->rowind = (int *)taucs_malloc(nent*sizeof(int));
662 out->values.d/*taucs_values*/ = (double *)taucs_malloc(nent*sizeof(double));
663 if (!(out->colptr) || !(out->rowind) || !(out->values.d/*taucs_values*/)) {
664 taucs_free(out->colptr);
665 taucs_free(out->rowind);
666 taucs_free(out->values.d/*taucs_values*/);
667 taucs_free(out);
668 return NULL;
669 }
670
671 out->n = n;
672 out->m = n;
673 out->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
674
675 return out;
676 }
677
678 #if 0
679 static
680 graph *ccs_matrix_to_graph(taucs_ccs_matrix *in)
681 {
682 graph *out;
683 int nent,n;
684 int j,ip;
685
686 nent = in->colptr[in->n];
687 out = construct_graph(nent);
688 if (!out) return NULL;
689
690 n = in->n;
691
692 out->n = n;
693 out->nent = nent;
694
695 for(j=0;j<n;j++) {
696 for(ip=in->colptr[j];ip<in->colptr[j+1];ip++) {
697 (out->edges)[ip].i = (in->rowind)[ip];
698 (out->edges)[ip].j = j;
699 (out->edges)[ip].v = (in->values.d/*taucs_values*/)[ip];
700 }
701 }
702
703 return(out);
704 }
705 #endif /* 0, we don't need this routine */
706
707 static
graph_resize(graph * a,int new_size)708 int graph_resize(graph *a,int new_size)
709 {
710 wedge* edges;
711
712 assert(new_size > a->max_size);
713
714 edges = (wedge*) taucs_malloc(new_size*sizeof(wedge));
715 if (!edges) {
716 return -1;
717 }
718
719 memcpy(edges,a->edges,a->max_size*sizeof(wedge));
720
721 taucs_free(a->edges);
722
723 a->edges=edges;
724
725 a->max_size = new_size;
726
727 return 0;
728 }
729
730 /************ LINKED LISTS *************/
731
732 typedef struct edg {
733 int entry_no;
734 struct edg *next;
735 } edge;
736
737 typedef struct linke {
738 edge **point;
739 edge *array;
740 } linked;
741
742 typedef struct thre {
743 int group_1;
744 int group_2;
745 int a;
746 int b;
747 double c;
748 byte already_connected;
749 byte completed_to_basis;
750 struct thre *next;
751 } three;
752
753 typedef struct si {
754 int group_1;
755 int group_2;
756 int a[2];
757 int b[2];
758 double c[2];
759 byte cross[2];
760 byte no_edges; /* number of edges connecting group_1 and group_2 (0,1 or 2) */
761 struct si *next;
762 } six;
763
taucs_check_diag_dominant_matrix(graph * A,int force_diagonal_dominance)764 int taucs_check_diag_dominant_matrix(graph *A, int force_diagonal_dominance)
765 {
766 int i;
767 double *sum;
768 int n;
769 int diagonally_dominant, all_nonpositive;
770
771 n = A->n;
772
773 sum = (double *)taucs_calloc(n,sizeof(double));
774 if (!sum) return -1;
775
776 Do(i,A->nent)
777 {
778 if ((A->edges)[i].i != (A->edges)[i].j)
779 {
780 sum[(A->edges)[i].i]-=fabs((A->edges)[i].v);
781 sum[(A->edges)[i].j]-=fabs((A->edges)[i].v);
782 }
783 else
784 {
785 sum[(A->edges)[i].i]+=(A->edges)[i].v;
786 if ((A->edges)[i].v < 0)
787 {
788 taucs_printf("ERROR! This matrix is not diagonally dominant. It has negative diagonals.\n");
789 /* taucs_free(sum); */
790 /* return -2; */
791 }
792 }
793 }
794
795 diagonally_dominant = 1; /* until proven otherwise */
796 all_nonpositive = 1;
797 Do(i,n)
798 {
799 if (sum[i] < -EPSILON) diagonally_dominant = 0;
800 if (sum[i] > EPSILON) all_nonpositive = 0;
801 }
802
803 if ((force_diagonal_dominance)&&(diagonally_dominant == 0)) {
804 int first_time = 1;
805 for(i=0;i<A->nent;i++)
806 {
807 if ((A->edges)[i].i == (A->edges)[i].j && sum[(A->edges)[i].i] <= EPSILON)
808 {
809 if (first_time) {
810 first_time=0;
811 taucs_printf("\t\tAMWB warning: perturbing to force diagonal dominance\n");
812 }
813 (A->edges)[i].v -= sum[ (A->edges)[i].i ];
814 if (all_nonpositive && (A->edges)[i].i == 0) {
815 taucs_printf("taucs warning: perturbing to ensure strict diagonal dominance\n");
816 (A->edges)[i].v += 0.1; /* arbitrary perturbation */
817 }
818 }
819 }
820 } else
821 if (diagonally_dominant == 0)
822 {
823 taucs_printf("ERROR! This matrix is not diagonally dominant. sum[%d] = %lf\n",i,sum[i]);
824 taucs_free(sum);
825 return -2;
826 }
827
828 taucs_free(sum);
829 return 0;
830 }
831
832 #if 0
833 static
834 double *analyze_graph(graph *A)
835 {
836 int i;
837 int t1,t2;
838 int n;
839 double *diag,t3;
840
841 n = A->n;
842 diag = (double *)taucs_calloc(n,sizeof(double));
843 if (!diag) return NULL;
844
845 Do(i,A->nent)
846 {
847 t1=(A->edges)[i].i;
848 t2=(A->edges)[i].j;
849 t3=(A->edges)[i].v;
850
851 if (t1 == t2)
852 diag[t1] += fabs(t3);
853 else
854 {
855 diag[t1] -= fabs(t3);
856 diag[t2] -= fabs(t3);
857 }
858 }
859 return(diag);
860 }
861 #endif /* 0, we don't need this routine */
862
863 /*********************************************************/
864 /* ccs diagnostics, row sums, and conversion to a graph */
865 /*********************************************************/
866
867 #define TAUCS_SYM_NOT_SYMLOWER 1
868 #define TAUCS_SYM_POS_OFFDIAGONALS 2
869 #define TAUCS_SYM_NEG_DIAGONALS 4
870 #define TAUCS_SYM_NOT_DIAGDOMINANT 8
871
872 static
ccs_matrix_to_graph_plus(taucs_ccs_matrix * in,int * diagnostics,double diag[],int force_diagonal_dominance)873 graph *ccs_matrix_to_graph_plus(taucs_ccs_matrix *in,
874 int* diagnostics,
875 double diag[],
876 int force_diagonal_dominance)
877 {
878 graph *out;
879 int nent,n;
880 int i,j,k,ip;
881 double v;
882 int negative_on_diagonal;
883 int positive_off_diagonal;
884 int not_diagonally_dominant;
885
886 *diagnostics = 0;
887
888 if (!(in->flags & TAUCS_SYMMETRIC) || !(in->flags & TAUCS_LOWER)) {
889 *diagnostics = TAUCS_SYM_NOT_SYMLOWER;
890 return NULL;
891 }
892
893 nent = in->colptr[in->n];
894 out = construct_graph(nent);
895 if (!out) return NULL;
896
897 n = in->n;
898
899 out->n = n;
900 out->nent = nent;
901
902 for (i=0; i<n; i++) diag[i] = 0.0;
903
904 negative_on_diagonal = 0;
905 positive_off_diagonal = 0;
906
907 for(j=0;j<n;j++) {
908 for(ip=in->colptr[j];ip<in->colptr[j+1];ip++) {
909 i = (in->rowind)[ip];
910 v = (in->values.d/*taucs_values*/)[ip];
911 (out->edges)[ip].i = i;
912 (out->edges)[ip].j = j;
913 (out->edges)[ip].v = v;
914
915 if (i == j) {
916 negative_on_diagonal |= (v < 0.0);
917
918 diag[i] += fabs(v);
919 } else {
920 positive_off_diagonal |= (v > 0.0);
921
922 diag[i] -= fabs(v);
923 diag[j] -= fabs(v);
924 }
925 }
926 }
927
928 if (force_diagonal_dominance) {
929 int strict_diagdominance = 0;
930 int first_time = 1;
931
932 for (i=0; i<n; i++)
933 strict_diagdominance |= (diag[i] > 0.0);
934
935 for(k=0;k<out->nent;k++) {
936 i = (out->edges)[k].i;
937 j = (out->edges)[k].j;
938 v = (out->edges)[k].v;
939
940 if (i == j && diag[i] < 0.0) {
941 if (first_time) {
942 first_time=0;
943 taucs_printf("taucs warning: perturbing to force diagonal dominance\n");
944 }
945 (out->edges)[k].v -= diag[i];
946 diag[i] = 0.0;
947 if (strict_diagdominance == 0 && i == 0) {
948 taucs_printf("taucs warning: perturbing to ensure strict diagonal dominance\n");
949 (out->edges)[k].v += 1e-8; /* arbitrary perturbation */
950 }
951 }
952 }
953
954 not_diagonally_dominant = 0;
955
956 } else {
957
958 not_diagonally_dominant = 0;
959 for (i=i; i<n; i++)
960 not_diagonally_dominant |= (diag[i] < -1e-12); /* arbitrary threashold */
961
962 }
963
964 *diagnostics = 0;
965
966 if (not_diagonally_dominant) *diagnostics |= TAUCS_SYM_NOT_DIAGDOMINANT;
967 if (negative_on_diagonal ) *diagnostics |= TAUCS_SYM_NEG_DIAGONALS;
968 if (positive_off_diagonal ) *diagnostics |= TAUCS_SYM_POS_OFFDIAGONALS;
969
970 return(out);
971 }
972
973
974 /*********************************************************/
975 /* */
976 /*********************************************************/
977
978
979 static
free_linked_list(linked * a)980 void free_linked_list(linked* a)
981 {
982 taucs_free(a->point);
983 taucs_free(a->array);
984 taucs_free(a);
985 }
986
987 static
free_linked_list_2(six * a)988 void free_linked_list_2(six *a)
989 {
990 if (a!=NULL)
991 {
992 free_linked_list_2(a->next);
993 taucs_free(a);
994 }
995
996 }
997
998 static
create_linked_list(graph * A,int n,int Anent,double * min,double * max)999 linked* create_linked_list(graph *A,int n,int Anent,double *min,double *max)
1000 {
1001 /* creates linked list which holds the off-diagonal entries of the sparse graph. */
1002 int i;
1003 edge *tmp;
1004 /*linked out; */
1005 linked* out;
1006 int free_place = 0;
1007
1008 *min = INF;
1009 *max = -INF;
1010
1011 out = (linked*) taucs_malloc(sizeof(linked));
1012 if (!out) {
1013 return NULL;
1014 }
1015 out->point = (edge **)taucs_calloc(n,sizeof(edge *));
1016 out->array = (edge *)taucs_calloc(2*Anent,sizeof(edge));
1017 if (!(out->point) || !(out->array)) {
1018 taucs_free(out->point);
1019 taucs_free(out->array);
1020 taucs_free(out);
1021 return NULL;
1022 }
1023
1024 Do(i,Anent)
1025 {
1026 if ((A->edges)[i].i != (A->edges)[i].j)
1027 {
1028 if (-(A->edges)[i].v < *min)
1029 *min = -(A->edges)[i].v;
1030 if (-(A->edges)[i].v > *max)
1031 *max = -(A->edges)[i].v;
1032
1033 tmp = &(out->array[free_place++]);
1034 tmp->entry_no = i;
1035 tmp->next = out->point[(A->edges)[i].i];
1036 out->point[(A->edges)[i].i] = tmp;
1037
1038 tmp = &(out->array[free_place++]);
1039 tmp->entry_no = i;
1040 tmp->next = out->point[(A->edges)[i].j];
1041 out->point[(A->edges)[i].j] = tmp;
1042 }
1043 }
1044
1045 return(out);
1046 }
1047
1048 static
create_linked_list_cluster(graph * A,int n,int Anent,double * min,double * max,int * partition,int * new_partition)1049 linked* create_linked_list_cluster(graph *A,int n,int Anent,double *min,double *max,int *partition,int *new_partition)
1050 {
1051 /* creates linked list which holds off-diagonal entries of the sparse graph.
1052 This linked list contains all the edges which connect vertices whose endpoints are
1053 in different sections in partition, but in the same section in new_partition.
1054 This will help us build trees within each section in new_partition. Each vertex
1055 in these trees will be a contracted section in partition */
1056 int i;
1057 edge *tmp;
1058
1059 linked* out = NULL;
1060 int free_place = 0;
1061
1062 *min = INF;
1063 *max = -INF;
1064
1065 out = (linked*) taucs_malloc(sizeof(linked));
1066 if (!out) {
1067 return NULL;
1068 }
1069 out->point = (edge **)taucs_calloc(n,sizeof(edge *));
1070 out->array = (edge *)taucs_calloc(2*Anent,sizeof(edge));
1071 if (!(out->point) || !(out->array)) {
1072 taucs_free(out->point);
1073 taucs_free(out->array);
1074 taucs_free(out);
1075 return NULL;
1076 }
1077
1078 Do(i,Anent)
1079 {
1080 if ((partition[(A->edges)[i].i] != partition[(A->edges)[i].j]) &&
1081 (new_partition[(A->edges)[i].i] == new_partition[(A->edges)[i].j]))
1082 {
1083 if (-(A->edges)[i].v < *min)
1084 *min = -(A->edges)[i].v;
1085 if (-(A->edges)[i].v > *max)
1086 *max = -(A->edges)[i].v;
1087
1088 tmp = &(out->array[free_place++]);
1089 tmp->entry_no = i;
1090 tmp->next = out->point[partition[(A->edges)[i].i]];
1091 out->point[partition[(A->edges)[i].i]] = tmp;
1092
1093 tmp = &(out->array[free_place++]);
1094 tmp->entry_no = i;
1095 tmp->next = out->point[partition[(A->edges)[i].j]];
1096 out->point[partition[(A->edges)[i].j]] = tmp;
1097 }
1098 }
1099
1100 return(out);
1101 }
1102
1103
1104 static
graph_to_ccs_matrix(graph * A)1105 taucs_ccs_matrix *graph_to_ccs_matrix(graph *A)
1106 {
1107 taucs_ccs_matrix *out;
1108 int n,nent,i,j1,j2;
1109 /*int count=0;*/
1110 int *tmp;
1111
1112 n = A->n;
1113 nent = A->nent;
1114
1115 tmp = (int *)taucs_malloc(n*sizeof(int));
1116 if (!tmp) return NULL;
1117
1118 out=construct_ccs_matrix(nent,n);
1119 if (!out) {
1120 taucs_free(tmp);
1121 return NULL;
1122 }
1123 out->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
1124
1125 Do(i,n)
1126 tmp[i] = 0;
1127 Do(i,nent)
1128 tmp[min((A->edges)[i].i,(A->edges)[i].j)]++;
1129 out->colptr[0] = 0;
1130 Do(i,n)
1131 out->colptr[i+1] = out->colptr[i] + tmp[i];
1132
1133 Do(i,n)
1134 tmp[i] = out->colptr[i];
1135
1136 Do(i,nent)
1137 {
1138 j1 = min((A->edges)[i].i , (A->edges)[i].j);
1139 j2 = max((A->edges)[i].i , (A->edges)[i].j);
1140 out->rowind[tmp[j1]]=j2;
1141 out->values.d/*taucs_values*/[tmp[j1]]=(A->edges)[i].v;
1142 tmp[j1]++;
1143 }
1144
1145 taucs_free(tmp);
1146 return(out);
1147 }
1148
1149 static
compute_sub_tree_sizes(int ver,int * first_child,int * next_child,int * sizes)1150 int compute_sub_tree_sizes(int ver,int *first_child,int *next_child,int *sizes)
1151 {
1152 int sum = 1,v;
1153
1154 if (first_child[ver] == -1)
1155 {
1156 sizes[ver] = 1;
1157 return(1);
1158 }
1159 else
1160 {
1161 v=first_child[ver];
1162 while(v != -1)
1163 {
1164 sum += compute_sub_tree_sizes(v,first_child,next_child,sizes);
1165 v = next_child[v];
1166 }
1167 }
1168 sizes[ver] = sum;
1169 return(sum);
1170 }
1171
1172 static
assign_group(int ver,int gr,int * first_child,int * next_child,int * groups)1173 void assign_group(int ver,int gr,int *first_child,int *next_child,int *groups)
1174 {
1175 int v;
1176
1177 groups[ver] = gr;
1178 if (first_child[ver] != -1)
1179 {
1180 v=first_child[ver];
1181 while(v != -1)
1182 {
1183 assign_group(v,gr,first_child,next_child,groups);
1184 v = next_child[v];
1185 }
1186 }
1187 }
1188
1189 static
create_children_arrays(int * pi,int n,int ** fc,int ** nc)1190 int create_children_arrays(int *pi,int n,int **fc,int **nc)
1191 {
1192 int *first_child,*next_child;
1193 int father,child,ch;
1194 int i;
1195
1196 first_child = (int *)taucs_malloc(n*sizeof(int));
1197 next_child = (int *)taucs_malloc(n*sizeof(int));
1198 if (!first_child || !next_child) {
1199 taucs_free(first_child);
1200 taucs_free(next_child);
1201 return -1;
1202 }
1203
1204 Do(i,n)
1205 first_child[i] = next_child[i] = -1;
1206
1207
1208 Do(i,n)
1209 {
1210 child = i;
1211 father = pi[i];
1212
1213 if (father != -1)
1214 {
1215 if (first_child[father] == -1)
1216 first_child[father] = child;
1217 else
1218 {
1219 ch = first_child[father];
1220 while(next_child[ch] != -1)
1221 ch = next_child[ch];
1222 next_child[ch] = child;
1223
1224 }
1225 }
1226 }
1227
1228
1229 *fc = first_child;
1230 *nc = next_child;
1231
1232 return 0; /* success */
1233 }
1234
1235 static
disconnect(int father,int child,int * first_child,int * next_child,int * pi)1236 void disconnect(int father,int child,
1237 int *first_child,int *next_child,int *pi)
1238 /* disconnect subtree whose root is 'child', from the tree */
1239 {
1240
1241 int oldfather;
1242 int v;
1243 /* int tmp;*/
1244
1245 oldfather = father;
1246
1247 assert(first_child[father] != -1);
1248
1249 if (first_child[father] == child)
1250 first_child[father] = next_child[child];
1251 else
1252 {
1253 v = first_child[father];
1254 while(next_child[v] != child)
1255 v = next_child[v];
1256
1257 next_child[v] = next_child[next_child[v]];
1258
1259 }
1260
1261 }
1262
1263 static
divide_to_groups(int r,int * first_child,int * next_child,int * pi,int * curr_group,int * groups,int * sub_tree_sizes,int root,double subgraphs,int n)1264 void divide_to_groups(int r,int *first_child,int *next_child,int *pi,int *curr_group,
1265 int *groups,int *sub_tree_sizes,int root,double subgraphs,
1266 int n)
1267 /* divides the vertices into different groups (divides the tree is subtrees) */
1268 {
1269 int v;
1270 double low;
1271
1272 low = max(1,((double)n/subgraphs));
1273
1274 if(first_child[r] != -1)
1275 {
1276 v=first_child[r];
1277 sub_tree_sizes[r] = 1;
1278 while(v != -1)
1279 {
1280 if (sub_tree_sizes[v] > low)
1281 divide_to_groups(v,first_child,next_child,pi,curr_group,groups,
1282 sub_tree_sizes,root,subgraphs,n);
1283
1284 if (sub_tree_sizes[v] >= low)
1285 {
1286 assign_group(v,*curr_group,first_child,next_child,groups);
1287 disconnect(r,v,first_child,next_child,pi);
1288 (*curr_group)++;
1289 }
1290 else
1291 sub_tree_sizes[r] += sub_tree_sizes[v];
1292 v = next_child[v];
1293 }
1294
1295 }
1296
1297 }
1298
1299
1300 static
DFS_visit(graph * precond,int r,byte * color,linked l,int * pi,int * visited)1301 void DFS_visit(graph *precond,int r,byte *color,linked l,int *pi,int *visited)
1302 {
1303 edge *p;
1304 int r1;
1305
1306 color[r] = 1;
1307 (*visited)++;
1308
1309 p = l.point[r];
1310 while (p != NULL)
1311 {
1312 /* this looks strange. Sivan */
1313 /*r1 = ivec1[ p->entry_no ]+ivec2[ p->entry_no ]-r;*/
1314 r1 = (precond->edges)[ p->entry_no ].i + (precond->edges)[ p->entry_no ].j - r;
1315 if (color[r1]==0)
1316 {
1317 DFS_visit(precond,r1,color,l,pi,visited);
1318 pi[r1] = r;
1319 }
1320 p = p->next;
1321 }
1322 }
1323
1324 static
make_perm(int perm[],int k)1325 void make_perm(int perm[],int k)
1326 {
1327 int i,tmp,tmp1;
1328
1329 for(i=0;i<k;i++)
1330 perm[i]=i;
1331
1332 for(i=0;i<k;i++)
1333 {
1334 tmp = rand()%(k-i);
1335 tmp1 = perm[i+tmp];
1336 perm[i+tmp]=perm[i];
1337 perm[i]=tmp1;
1338 }
1339 }
1340
1341 #define USE_HEAPSORT
1342
1343 static
1344 taucs_ccs_matrix*
amwb_preconditioner_create(graph * mtxA,double * diag,int rnd,double subgraphs)1345 amwb_preconditioner_create(graph *mtxA, double* diag,
1346 int rnd,
1347 double subgraphs)
1348 /*
1349 amwb_preconditioner_create(taucs_ccs_matrix *taucs_ccs_mtxA,
1350 int rnd,
1351 double subgraphs)
1352 */
1353 {
1354 taucs_ccs_matrix *out;
1355
1356 /* graph *mtxA,*mtxA_tmp;*/
1357 graph *mtxA_tmp;
1358 graph *precond;
1359
1360 int i;
1361 /*int j,tmp,entry;*/
1362 /*
1363 int *ivec1, *ivec2 ;
1364 int *ivec1_p, *ivec2_p ;
1365 double *dvec;
1366 double *dvec_p;
1367 */
1368
1369 int size;
1370 int *pi;
1371 /*edge **array;*/
1372 /*edge *p;*/
1373 #ifdef USE_HEAPSORT
1374 heap Ah;
1375 heap Bh;
1376 heap h;
1377 #endif
1378 int Bent;
1379 /*int row,col;*/
1380 double weight;
1381 int *first_child,*next_child;
1382 int *groups,*sub_tree_sizes;
1383 int curr_group;
1384 int n,Anent;
1385 int precond_maxnent,chunk;
1386 /*double *diag;*/
1387 byte *closed_cycle,closed_cycle_x,closed_cycle_y;
1388 byte *color;
1389 byte *already_added;
1390 int count=0;
1391 int edge_sign; /* byte edge_sign;*/
1392 int u,v,x,y,un_root,r;
1393 linked* l;
1394 int *perm;
1395 int visited,*roots,rcount;
1396 int basis_Bent=0,step2_Bent=0;
1397 three *complete_subgraph;
1398 six **pairs;
1399 /*char bool;*/
1400 /* FILE *graph_file0,*graph_file1,*graph_file2,*graph_file3,*group_file; */
1401
1402 double wtime;
1403 double wtime_sort = 0.0;
1404 double wtime_global_basis;
1405 double wtime_treepartition;
1406 double wtime_component_bases;
1407 double wtime_pair_bases;
1408 double wtime_total;
1409 double dummy;
1410
1411 wtime_total = taucs_wtime();
1412
1413 /* graph_file0 = fopen("graphfile0.txt","w"); */
1414 /* graph_file1 = fopen("graphfile1.txt","w"); */
1415 /* graph_file2 = fopen("graphfile2.txt","w"); */
1416 /* graph_file3 = fopen("graphfile3.txt","w"); */
1417 /* group_file = fopen("groupfile.txt" ,"w"); */
1418
1419
1420 /********************************************************/
1421 /* */
1422 /********************************************************/
1423
1424 /********************************************************/
1425 /* convert matrix to a graph */
1426 /********************************************************/
1427
1428 /*** ALLOCATED: NONE ***/
1429
1430 /*
1431 wtime = taucs_wtime();
1432 mtxA = ccs_matrix_to_graph(taucs_ccs_mtxA);
1433 if (!mtxA) {
1434 return NULL;
1435 }
1436 wtime = taucs_wtime() - wtime;
1437 taucs_printf("\t\tAMWB matrix-to-graph = %.3f seconds\n",wtime);
1438 */
1439
1440 /********************************************************/
1441 /* check that the matrix is diagonally dominant */
1442 /********************************************************/
1443
1444 /*** ALLOCATED: mtxA ***/
1445
1446 #if 0
1447 wtime = taucs_wtime();
1448 i = taucs_check_diag_dominant_matrix(mtxA,1 /* force diagonal dominance */);
1449 if (i == -1) {
1450 free_graph(mtxA);
1451 return NULL;
1452 }
1453 if (i == -2) {
1454 free_graph(mtxA);
1455 return taucs_ccs_mtxA; /* not diagonally dominant */
1456 }
1457 wtime = taucs_wtime() - wtime;
1458 taucs_printf("\t\tAMWB check-diag-dominance = %.3f seconds\n",wtime);
1459 #endif
1460
1461 n = mtxA->n;
1462
1463 /********************************************************/
1464 /* generate random permutation and permute vertices */
1465 /********************************************************/
1466
1467 wtime = taucs_wtime();
1468 perm = (int *)taucs_malloc(mtxA->nent*sizeof(int));
1469 if (!perm) {
1470 free_graph(mtxA);
1471 return NULL;
1472 }
1473
1474 /*** ALLOCATED: mtxA,perm ***/
1475
1476 make_perm(perm,mtxA->nent);
1477
1478 mtxA_tmp = construct_graph(mtxA->max_size);
1479 if (!mtxA_tmp) {
1480 free_graph(mtxA);
1481 taucs_free(perm);
1482 return NULL;
1483 }
1484
1485 /*** ALLOCATED: mtxA,perm,mtxA_tmp ***/
1486
1487 mtxA_tmp->nent = mtxA->nent;
1488 mtxA_tmp->n = mtxA->n;
1489 Do(i,mtxA->nent)
1490 {
1491 /*
1492 mtxA_tmp->ivec1[i] = mtxA->ivec1[perm[i]];
1493 mtxA_tmp->ivec2[i] = mtxA->ivec2[perm[i]];
1494 mtxA_tmp->dvec[i] = mtxA->dvec[perm[i]];
1495 */
1496
1497 mtxA_tmp->edges[i].i = mtxA->edges[perm[i]].i;
1498 mtxA_tmp->edges[i].j = mtxA->edges[perm[i]].j;
1499 mtxA_tmp->edges[i].v = mtxA->edges[perm[i]].v;
1500 }
1501
1502 taucs_free(perm);
1503 free_graph(mtxA);
1504
1505 wtime = taucs_wtime() - wtime;
1506 taucs_printf("\t\tAMWB random permute = %.3f seconds\n",wtime);
1507
1508 /********************************************************/
1509 /* compute and remember row weights */
1510 /********************************************************/
1511
1512 wtime = taucs_wtime();
1513
1514 /*** ALLOCATED: mtxA_tmp ***/
1515
1516 /*
1517 diag = analyze_graph(mtxA_tmp);
1518 if (!diag) {
1519 free_graph(mtxA_tmp);
1520 return NULL;
1521 }
1522 */
1523
1524 wtime = taucs_wtime() - wtime;
1525 taucs_printf("\t\tAMWB row weights = %.3f seconds\n",wtime);
1526
1527 /********************************************************/
1528 /* allocate vectors */
1529 /********************************************************/
1530
1531 /*** ALLOCATED: mtxA_tmp,diag ***/
1532
1533 Anent = mtxA_tmp->nent;
1534
1535 already_added = (byte *)taucs_calloc(Anent,sizeof(byte));
1536 pi = (int *) taucs_malloc(n*sizeof(int));
1537 if (!already_added || !pi) {
1538 free_graph(mtxA_tmp);
1539 taucs_free(diag);
1540 taucs_free(already_added);
1541 taucs_free(pi);
1542 return NULL;
1543 }
1544
1545 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi ***/
1546
1547 Do(i,n)
1548 pi[i] = -1;
1549
1550 /********************************************************/
1551 /* construct empty preconditioner */
1552 /********************************************************/
1553
1554 precond_maxnent = 3*n;
1555
1556 taucs_printf("allocating space for %d entries in precond\n",precond_maxnent);fflush(stdout);
1557
1558 precond = construct_graph(precond_maxnent);
1559 if (!precond) {
1560 free_graph(mtxA_tmp);
1561 taucs_free(diag);
1562 taucs_free(already_added);
1563 taucs_free(pi);
1564 return NULL;
1565 }
1566
1567 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond ***/
1568
1569 precond->n=mtxA_tmp->n;
1570
1571 /*
1572 ivec1_p = precond->ivec1 ;
1573 ivec2_p = precond->ivec2 ;
1574 dvec_p = precond->dvec ;
1575 */
1576
1577 Bent = 0;
1578
1579 /********************************************************/
1580 /* allocate vectors */
1581 /********************************************************/
1582
1583 wtime_global_basis = taucs_wtime();
1584
1585 closed_cycle = (byte *)taucs_calloc(n,sizeof(byte));
1586 if (!closed_cycle) {
1587 free_graph(mtxA_tmp);
1588 free_graph(precond);
1589 taucs_free(diag);
1590 taucs_free(already_added);
1591 taucs_free(pi);
1592 return NULL;
1593 }
1594
1595 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,closed_cycle ***/
1596
1597 /* Variation on Kruskal - Introduction to Algorithms page 505 */
1598
1599 /********************************************************/
1600 /* initialize union-find */
1601 /********************************************************/
1602
1603 if (unionfind_init(n) == -1) {
1604 free_graph(mtxA_tmp);
1605 free_graph(precond);
1606 taucs_free(diag);
1607 taucs_free(already_added);
1608 taucs_free(pi);
1609 taucs_free(closed_cycle);
1610 return NULL;
1611 }
1612
1613 /********************************************************/
1614 /* sort edges of matrix */
1615 /********************************************************/
1616
1617 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,closed_cycle,UF ***/
1618
1619 wtime = taucs_wtime();
1620 #ifdef USE_HEAPSORT
1621
1622
1623 /*
1624 size = heap_sort(Anent,&h,mtxA_tmp);
1625 if (size == -1) {
1626 */
1627 if (pqueue_create(&Ah,Anent) == -1) {
1628 free_graph(mtxA_tmp);
1629 free_graph(precond);
1630 taucs_free(diag);
1631 taucs_free(already_added);
1632 taucs_free(pi);
1633 taucs_free(closed_cycle);
1634 unionfind_free();
1635 return NULL;
1636 }
1637 size = pqueue_fill(&Ah,mtxA_tmp);
1638 #else
1639 assert(Anent == mtxA_tmp->nent);
1640 size = mtxA_tmp->nent;
1641 graph_sort(mtxA_tmp);
1642 #endif
1643 wtime = taucs_wtime() - wtime;
1644 taucs_printf("\t\tAMWB sort(%d) = %.3f seconds\n",Anent,wtime);
1645 wtime_sort += wtime;
1646
1647 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,closed_cycle,UF,heap ***/
1648
1649 /********************************************************/
1650 /* build a basis for the entire graph */
1651 /********************************************************/
1652
1653 Do(i,size)
1654 {
1655 if (count == n)
1656 break;
1657
1658 #ifdef USE_HEAPSORT
1659 u = mtxA_tmp->edges[Ah.edges[i]].i;
1660 v = mtxA_tmp->edges[Ah.edges[i]].j;
1661 weight = mtxA_tmp->edges[Ah.edges[i]].v;
1662 #else
1663 u = mtxA_tmp->edges[i].i;
1664 v = mtxA_tmp->edges[i].j;
1665 weight = mtxA_tmp->edges[i].v;
1666 if (u==v) continue;
1667 #endif
1668
1669 edge_sign = (weight>0);
1670
1671 x = find_set(u);
1672 y = find_set(v);
1673
1674 if (x!=y)
1675 {
1676 /* printf("different trees\n"); */
1677 if (!((closed_cycle[x])&&(closed_cycle[y])))
1678 {
1679 count++;
1680 /* printf("(%d,%d) - %lf\n",u,v,weight); */
1681 /*
1682 ivec1_p[Bent] = u;
1683 ivec2_p[Bent] = v;
1684 dvec_p[Bent] = weight;
1685 */
1686 precond->edges[Bent].i = u;
1687 precond->edges[Bent].j = v;
1688 precond->edges[Bent].v = weight;
1689
1690 Bent++;
1691
1692 diag[u] += fabs(weight);
1693 diag[v]+=fabs(weight);
1694
1695 #ifdef USE_HEAPSORT
1696 already_added[Ah.edges[i]] = 1;
1697 #else
1698 already_added[i] = 1;
1699 #endif
1700 un_root = Union(u,v,x,y,edge_sign);
1701 closed_cycle[un_root] = closed_cycle[x] | closed_cycle[y];
1702 }
1703 /* else
1704 {
1705 printf("cannot add (%d,%d) - %lf - both trees already have a cycle\n",u,v,weight);
1706 }
1707 */
1708 }
1709 else
1710 {
1711 /* printf("same tree\n"); */
1712 if ((edge_sign != (label[u]^label[v])) && (closed_cycle[x]==0))
1713 {
1714 count++;
1715 /* printf("(%d,%d) - %lf\n",u,v,weight); */
1716 /*
1717 ivec1_p[Bent] = u;
1718 ivec2_p[Bent] = v;
1719 dvec_p[Bent] = weight;
1720 */
1721 precond->edges[Bent].i = u;
1722 precond->edges[Bent].j = v;
1723 precond->edges[Bent].v = weight;
1724
1725 Bent++;
1726
1727 diag[u] += fabs(weight);
1728 diag[v]+=fabs(weight);
1729
1730 #ifdef USE_HEAPSORT
1731 already_added[Ah.edges[i]] = 1;
1732 #else
1733 already_added[i] = 1;
1734 #endif
1735 closed_cycle[x] = 1;
1736 }
1737 /* else
1738 {
1739 if (closed_cycle[x]==1)
1740 printf("cannot add (%d,%d) - %lf - tree already contains cycle\n",u,v,weight);
1741 else
1742 printf("cannot add (%d,%d) - %lf - it closes a positive cycle\n",u,v,weight);
1743 }
1744 */
1745 }
1746 }
1747
1748 /********************************************************/
1749 /* the preconditioner is now a max-weight-basis */
1750 /********************************************************/
1751
1752 wtime_global_basis = taucs_wtime() - wtime_global_basis;
1753 taucs_printf("\t\tAMWB global basis = %.3f seconds\n",wtime_global_basis);
1754
1755 taucs_free(closed_cycle);
1756 #ifdef USE_HEAPSORT
1757 /*free_heap(h);*/
1758 #endif
1759 unionfind_free();
1760
1761 precond->nent = Bent;
1762 basis_Bent = Bent;
1763
1764 /********************************************************/
1765 /* break into subgraphs */
1766 /********************************************************/
1767
1768 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond ***/
1769
1770 l = create_linked_list(precond,n,Bent,&dummy,&dummy);
1771 if (!l) {
1772 free_graph(mtxA_tmp);
1773 free_graph(precond);
1774 taucs_free(diag);
1775 taucs_free(already_added);
1776 taucs_free(pi);
1777 return NULL;
1778 }
1779
1780 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1781
1782 color = (byte *)taucs_calloc(n,sizeof(byte));
1783 roots = (int *)taucs_malloc(n*sizeof(int));
1784 if (!color || !roots) {
1785 free_graph(mtxA_tmp);
1786 free_graph(precond);
1787 taucs_free(diag);
1788 taucs_free(already_added);
1789 taucs_free(pi);
1790 free_linked_list(l);
1791 return NULL;
1792 }
1793
1794 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,color,roots ***/
1795
1796 wtime_treepartition = taucs_wtime();
1797
1798 visited = 0;
1799 rcount = 0;
1800 /*
1801 while(visited<n)
1802 {
1803 r = rand()%n;
1804 while (color[r])
1805 r = rand()%n;
1806 roots[rcount++] = r;
1807 pi[r] = -1;
1808 DFS_visit(precond,r,color,*l,pi,&visited);
1809 }
1810 */
1811
1812 for (r=0; r<n; r++) {
1813 if (color[r] != 0) continue;
1814 roots[rcount++] = r;
1815 pi[r] = -1;
1816 DFS_visit(precond,r,color,*l,pi,&visited);
1817 }
1818
1819 taucs_free(color);
1820
1821 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,roots ***/
1822
1823 /* pi now contains the parent array of the tree */
1824
1825 if (create_children_arrays(pi,n,&first_child,&next_child) == -1) {
1826 free_graph(mtxA_tmp);
1827 free_graph(precond);
1828 taucs_free(diag);
1829 taucs_free(already_added);
1830 taucs_free(pi);
1831 taucs_free(roots);
1832 free_linked_list(l);
1833 return NULL;
1834 }
1835
1836 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,roots,FC,NC ***/
1837
1838 /* the first_child/next_child arrays enable us to find the children
1839 of any vertex in the tree */
1840
1841 groups = (int *)taucs_malloc(n*sizeof(int));
1842 sub_tree_sizes = (int *)taucs_malloc(n*sizeof(int));
1843 if(!groups || !sub_tree_sizes) {
1844 free_graph(mtxA_tmp);
1845 free_graph(precond);
1846 taucs_free(diag);
1847 taucs_free(already_added);
1848 taucs_free(pi);
1849 taucs_free(roots);
1850 taucs_free(groups);
1851 taucs_free(sub_tree_sizes);
1852 taucs_free(first_child);
1853 taucs_free(next_child);
1854 free_linked_list(l);
1855 return NULL;
1856 }
1857
1858 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,roots,FC,NC ***/
1859 /*** ALLOCATED: groups,sub_tree_sizes ***/
1860
1861 Do(i,n)
1862 groups[i] = -1;
1863
1864 curr_group = 0;
1865
1866 Do(i,rcount)
1867 {
1868 r = roots[i];
1869 compute_sub_tree_sizes(r,first_child,next_child,sub_tree_sizes);
1870 /* now for every vertex v in the tree, sub_tree_sizes[v] is the size
1871 of the subtree whose root is v */
1872
1873 divide_to_groups(r,first_child,next_child,pi,&curr_group,groups,
1874 sub_tree_sizes,r,subgraphs,n);
1875 if ((sub_tree_sizes[r]<((double)n/subgraphs))&&(curr_group>0))
1876 curr_group--;
1877 assign_group(r,curr_group,first_child,next_child,groups);
1878 curr_group++;
1879 }
1880
1881 taucs_printf("actual number of subgraphs = %ld\n",curr_group);fflush(stdout);
1882 /* now the tree is divided into linked groups */
1883
1884 chunk = max(min((curr_group*(curr_group-1)/2)/10,5000),10000);
1885
1886 taucs_free(roots);
1887 taucs_free(first_child);
1888 taucs_free(next_child);
1889 taucs_free(sub_tree_sizes);
1890
1891 wtime_treepartition = taucs_wtime() - wtime_treepartition;
1892 taucs_printf("\t\tAMWB treepartition = %.3f seconds\n",wtime_treepartition);
1893
1894
1895
1896
1897 wtime_component_bases = taucs_wtime();
1898
1899 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1900 /*** ALLOCATED: groups ***/
1901
1902 precond->nent = Bent;
1903
1904 /********************************************************/
1905 /* complete each subgraph into a basis */
1906 /********************************************************/
1907
1908 complete_subgraph = (three *)taucs_calloc(curr_group,sizeof(three));
1909 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1910 /*** ALLOCATED: groups ***/
1911 if (!complete_subgraph) {
1912 free_graph(mtxA_tmp);
1913 free_graph(precond);
1914 taucs_free(diag);
1915 taucs_free(already_added);
1916 taucs_free(pi);
1917 taucs_free(groups);
1918 free_linked_list(l);
1919 return NULL;
1920 }
1921
1922 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1923 /*** ALLOCATED: groups,complete_subgraph ***/
1924
1925 /* For each subgraph, complete_subgraph will store the edge needed
1926 to complete the subgraph into a basis (if such an edge exists) */
1927
1928 /* Variation on Kruskal - Introduction to Algorithms page 505 */
1929 if (unionfind_init(n) == -1) {
1930 free_graph(mtxA_tmp);
1931 free_graph(precond);
1932 free_linked_list(l);
1933 taucs_free(diag);
1934 taucs_free(already_added);
1935 taucs_free(pi);
1936 taucs_free(groups);
1937 taucs_free(complete_subgraph);
1938 return NULL;
1939 }
1940
1941 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1942 /*** ALLOCATED: groups,complete_subgraph,UF ***/
1943
1944 closed_cycle = (byte *)taucs_calloc(n,sizeof(byte));
1945 if (!closed_cycle) {
1946 free_graph(mtxA_tmp);
1947 free_graph(precond);
1948 free_linked_list(l);
1949 unionfind_free();
1950 taucs_free(diag);
1951 taucs_free(already_added);
1952 taucs_free(pi);
1953 taucs_free(groups);
1954 taucs_free(complete_subgraph);
1955 return NULL;
1956 }
1957
1958 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
1959 /*** ALLOCATED: groups,complete_subgraph,UF,closed_cycle ***/
1960
1961 wtime = taucs_wtime();
1962 #ifdef USE_HEAPSORT
1963 /*
1964 size = heap_sort(Bent,&h,precond);
1965 if (size == -1) {
1966 free_graph(mtxA_tmp);
1967 free_graph(precond);
1968 free_linked_list(l);
1969 unionfind_free();
1970 taucs_free(diag);
1971 taucs_free(already_added);
1972 taucs_free(pi);
1973 taucs_free(groups);
1974 taucs_free(complete_subgraph);
1975 taucs_free(closed_cycle);
1976 return NULL;
1977 }
1978 */
1979 if (pqueue_create(&Bh,Anent) == -1) {
1980 free_graph(mtxA_tmp);
1981 free_graph(precond);
1982 taucs_free(diag);
1983 taucs_free(already_added);
1984 taucs_free(pi);
1985 taucs_free(closed_cycle);
1986 unionfind_free();
1987 return NULL;
1988 }
1989 size = pqueue_fill(&Bh,precond);
1990 #else
1991 assert(Bent == precond->nent);
1992 size = precond->nent;
1993 graph_sort(precond);
1994 #endif
1995 wtime = taucs_wtime() - wtime;
1996 taucs_printf("\t\tAMWB sort 2(%d) = %.3f seconds\n",Bent,wtime);
1997 wtime_sort += wtime;
1998
1999 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
2000 /*** ALLOCATED: groups,complete_subgraph,UF,closed_cycle,heap ***/
2001
2002 Do(i,size)
2003 {
2004 #ifdef USE_HEAPSORT
2005 u = precond->edges[Bh.edges[i]].i;
2006 v = precond->edges[Bh.edges[i]].j;
2007 #else
2008 u = precond->edges[i].i;
2009 v = precond->edges[i].j;
2010 if (u==v) continue;
2011 #endif
2012
2013
2014 if (groups[u] == groups[v])
2015 {
2016 #ifdef USE_HEAPSORT
2017 weight = precond->edges[Bh.edges[i]].v;
2018 #else
2019 weight = precond->edges[i].v;
2020 #endif
2021
2022 edge_sign = (weight>0);
2023
2024 x = find_set(u);
2025 y = find_set(v);
2026
2027 if (x!=y)
2028 {
2029 if (!((closed_cycle[x])&&(closed_cycle[y])))
2030 {
2031 un_root = Union(u,v,x,y,edge_sign);
2032 closed_cycle[un_root] = closed_cycle[x] | closed_cycle[y];
2033 }
2034 else
2035 assert(0); /* this is a subgraph of a basis */
2036 }
2037 else
2038 {
2039 if ((edge_sign != (label[u]^label[v])) && (closed_cycle[x]==0))
2040 closed_cycle[x] = 1;
2041 }
2042
2043 }
2044 }
2045
2046 wtime = taucs_wtime();
2047 #ifdef USE_HEAPSORT
2048 /*
2049 free_heap(h);
2050 size = heap_sort(Anent,&h,mtxA_tmp);
2051 if (size == -1) {
2052 free_graph(mtxA_tmp);
2053 free_graph(precond);
2054 free_linked_list(l);
2055 unionfind_free();
2056 taucs_free(diag);
2057 taucs_free(already_added);
2058 taucs_free(pi);
2059 taucs_free(groups);
2060 taucs_free(complete_subgraph);
2061 taucs_free(closed_cycle);
2062 return NULL;
2063 }
2064 */
2065 /*size = pqueue_fill(&Ah,mtxA_tmp);*/
2066 size = Ah.heap_size;
2067 #else
2068 assert(Anent == mtxA_tmp->nent);
2069 size = mtxA_tmp->nent;
2070 graph_sort(mtxA_tmp);
2071 #endif
2072 wtime = taucs_wtime() - wtime;
2073 taucs_printf("\t\tAMWB sort(%d) = %.3f seconds\n",Anent,wtime);
2074 wtime_sort += wtime;
2075
2076 Do(i,size)
2077 {
2078 #ifdef USE_HEAPSORT
2079 u = mtxA_tmp->edges[Ah.edges[i]].i;
2080 v = mtxA_tmp->edges[Ah.edges[i]].j;
2081 #else
2082 u = mtxA_tmp->edges[i].i;
2083 v = mtxA_tmp->edges[i].j;
2084 if (u==v) continue;
2085 #endif
2086
2087 if (groups[u] == groups[v])
2088 {
2089 #ifdef USE_HEAPSORT
2090 weight = mtxA_tmp->edges[Ah.edges[i]].v;
2091 #else
2092 weight = mtxA_tmp->edges[i].v;
2093 #endif
2094
2095 edge_sign = (weight>0);
2096
2097 x = find_set(u);
2098 y = find_set(v);
2099
2100 if (x!=y)
2101 {
2102 if (!((closed_cycle[x])&&(closed_cycle[y])))
2103 {
2104 /*
2105 ivec1_p[Bent] = u;
2106 ivec2_p[Bent] = v;
2107 dvec_p[Bent] = weight;
2108 */
2109
2110 precond->edges[Bent].i = u;
2111 precond->edges[Bent].j = v;
2112 precond->edges[Bent].v = weight;
2113
2114 Bent++;
2115
2116 diag[u] += fabs(weight);
2117 diag[v]+=fabs(weight);
2118
2119 assert(complete_subgraph[groups[u]].completed_to_basis==0);
2120 complete_subgraph[groups[u]].completed_to_basis = 1;
2121 complete_subgraph[groups[u]].a = u;
2122 complete_subgraph[groups[u]].b = v;
2123 complete_subgraph[groups[u]].c = weight;
2124
2125 un_root = Union(u,v,x,y,edge_sign);
2126 closed_cycle[un_root] = closed_cycle[x] | closed_cycle[y];
2127 }
2128 }
2129 else
2130 {
2131 if ((edge_sign != (label[u]^label[v])) && (closed_cycle[x]==0))
2132 {
2133 /*
2134 ivec1_p[Bent] = u;
2135 ivec2_p[Bent] = v;
2136 dvec_p[Bent] = weight;
2137 */
2138 precond->edges[Bent].i = u;
2139 precond->edges[Bent].j = v;
2140 precond->edges[Bent].v = weight;
2141
2142 Bent++;
2143
2144 diag[u] += fabs(weight);
2145 diag[v] += fabs(weight);
2146
2147 assert(complete_subgraph[groups[u]].completed_to_basis==0);
2148 complete_subgraph[groups[u]].completed_to_basis = 1;
2149 complete_subgraph[groups[u]].a = u;
2150 complete_subgraph[groups[u]].b = v;
2151 complete_subgraph[groups[u]].c = weight;
2152
2153 closed_cycle[x] = 1;
2154 }
2155 }
2156
2157 }
2158 }
2159 #ifdef USE_HEAPSORT
2160 /* free_heap(h);*/
2161 #endif
2162 taucs_free(closed_cycle);
2163 unionfind_free();
2164
2165 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
2166 /*** ALLOCATED: groups,complete_subgraph ***/
2167
2168 precond->nent = Bent;
2169
2170 wtime_component_bases = taucs_wtime() - wtime_component_bases;
2171 taucs_printf("\t\tAMWB component bases = %.3f seconds\n",wtime_component_bases);
2172
2173
2174
2175 wtime_pair_bases = taucs_wtime();
2176
2177 step2_Bent = Bent;
2178
2179 /* COMPLETE EACH PAIR OF SUBGRAPHS INTO A BASIS */
2180 if (curr_group>1) {
2181 pairs = (six **)taucs_calloc(curr_group,sizeof(six *));
2182 closed_cycle = (byte *)taucs_calloc(n,sizeof(byte));
2183 if(!pairs || !closed_cycle) {
2184 free_graph(mtxA_tmp);
2185 free_graph(precond);
2186 free_linked_list(l);
2187 taucs_free(diag);
2188 taucs_free(already_added);
2189 taucs_free(pi);
2190 taucs_free(groups);
2191 taucs_free(complete_subgraph);
2192 taucs_free(pairs);
2193 taucs_free(closed_cycle);
2194 }
2195
2196 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
2197 /*** ALLOCATED: groups,complete_subgraph,pairs,closed_cycle ***/
2198
2199 if (unionfind_init(n) == -1) {
2200 free_graph(mtxA_tmp);
2201 free_graph(precond);
2202 free_linked_list(l);
2203 taucs_free(diag);
2204 taucs_free(already_added);
2205 taucs_free(pi);
2206 taucs_free(groups);
2207 taucs_free(complete_subgraph);
2208 taucs_free(pairs);
2209 taucs_free(closed_cycle);
2210 return NULL;
2211 }
2212
2213 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,UF ***/
2214 /*** ALLOCATED: groups,complete_subgraph,pairs,closed_cycle ***/
2215
2216 wtime = taucs_wtime();
2217 #ifdef USE_HEAPSORT
2218 /*
2219 size = heap_sort(basis_Bent,&h,precond);
2220 if (size == -1) {
2221 free_graph(mtxA_tmp);
2222 free_graph(precond);
2223 free_linked_list(l);
2224 unionfind_free();
2225 taucs_free(diag);
2226 taucs_free(already_added);
2227 taucs_free(pi);
2228 taucs_free(groups);
2229 taucs_free(complete_subgraph);
2230 taucs_free(pairs);
2231 taucs_free(closed_cycle);
2232 return NULL;
2233 }
2234 */
2235 size = pqueue_fill(&Bh,precond);
2236
2237 #else
2238 assert(basis_Bent == precond->nent);
2239 size = precond->nent;
2240 graph_sort(precond);
2241 #endif
2242 wtime = taucs_wtime() - wtime;
2243 taucs_printf("\t\tAMWB sort(%d) = %.3f seconds\n",basis_Bent,wtime);
2244 wtime_sort += wtime;
2245
2246 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked,UF ***/
2247 /*** ALLOCATED: groups,complete_subgraph,pairs,closed_cycle,heap ***/
2248
2249 Do(i,size)
2250 {
2251 #ifdef USE_HEAPSORT
2252 u = precond->edges[Bh.edges[i]].i;
2253 v = precond->edges[Bh.edges[i]].j;
2254 #else
2255 u = precond->edges[i].i;
2256 v = precond->edges[i].j;
2257 if (u==v) continue;
2258 #endif
2259
2260 if (groups[u] == groups[v])
2261 {
2262 #ifdef USE_HEAPSORT
2263 weight = precond->edges[Bh.edges[i]].v;
2264 #else
2265 weight = precond->edges[i].v;
2266 #endif
2267
2268 edge_sign = (weight>0);
2269
2270 x = find_set(u);
2271 y = find_set(v);
2272
2273 if (x!=y)
2274 {
2275 if (!((closed_cycle[x])&&(closed_cycle[y])))
2276 {
2277 un_root = Union(u,v,x,y,edge_sign);
2278 closed_cycle[un_root] = closed_cycle[x] | closed_cycle[y];
2279 }
2280 }
2281 else
2282 {
2283 if ((edge_sign != (label[u]^label[v])) && (closed_cycle[x]==0))
2284 closed_cycle[x] = 1;
2285 }
2286
2287 }
2288 }
2289
2290 wtime = taucs_wtime();
2291 #ifdef USE_HEAPSORT
2292 /*
2293 free_heap(h);
2294 size = heap_sort(Anent,&h,mtxA_tmp);
2295 if (size == -1) {
2296 free_graph(mtxA_tmp);
2297 free_graph(precond);
2298 free_linked_list(l);
2299 unionfind_free();
2300 taucs_free(diag);
2301 taucs_free(already_added);
2302 taucs_free(pi);
2303 taucs_free(groups);
2304 taucs_free(complete_subgraph);
2305 taucs_free(pairs);
2306 taucs_free(closed_cycle);
2307 return NULL;
2308 }
2309 */
2310 /*size = pqueue_fill(&Ah,mtxA_tmp);*/
2311 size = Ah.heap_size;
2312 #else
2313 assert(Anent == mtxA_tmp->nent);
2314 size = mtxA_tmp->nent;
2315 graph_sort(mtxA_tmp);
2316 #endif
2317 wtime = taucs_wtime() - wtime;
2318 taucs_printf("\t\tAMWB sort(%d) = %.3f seconds\n",Anent,wtime);
2319 wtime_sort += wtime;
2320
2321 Do(i,size)
2322 {
2323 int g1,g2;
2324 six *p;
2325 six *tmp;
2326
2327
2328 #ifdef USE_HEAPSORT
2329 u = mtxA_tmp->edges[Ah.edges[i]].i;
2330 v = mtxA_tmp->edges[Ah.edges[i]].j;
2331 #else
2332 u = mtxA_tmp->edges[i].i;
2333 v = mtxA_tmp->edges[i].j;
2334 if (u==v) continue;
2335 #endif
2336
2337 if (Bent == precond_maxnent)
2338 {
2339 precond_maxnent += chunk;
2340 taucs_printf("adding space for %d entries in precond\n",chunk);fflush(stdout);
2341 if (graph_resize(precond,precond_maxnent) == -1) {
2342 int i_local;
2343 Do(i_local,curr_group)
2344 free_linked_list_2(pairs[i_local]);
2345 /* precond has not been freed */
2346 free_graph(mtxA_tmp);
2347 free_graph(precond);
2348 free_linked_list(l);
2349 unionfind_free();
2350 #ifdef USE_HEAPSORT
2351 free_heap(h);
2352 #endif
2353 taucs_free(diag);
2354 taucs_free(already_added);
2355 taucs_free(pi);
2356 taucs_free(groups);
2357 taucs_free(complete_subgraph);
2358 taucs_free(pairs);
2359 taucs_free(closed_cycle);
2360 return NULL;
2361 }
2362
2363 /*
2364 ivec1_p = precond->ivec1;
2365 ivec2_p = precond->ivec2;
2366 dvec_p = precond->dvec;
2367 */
2368 }
2369
2370 g1 = min(groups[u],groups[v]); g2 = max(groups[u],groups[v]);
2371
2372 if (g1 != g2)
2373 {
2374 p = pairs[(g1+g2)%curr_group];
2375 while (p!=NULL)
2376 {
2377 if ((p->group_1 == g1)&&(p->group_2 == g2))
2378 goto after2;
2379 p = p->next;
2380 }
2381 after2:
2382
2383 if (p == NULL)
2384 {
2385 tmp = (six *)taucs_calloc(1,sizeof(six));
2386 if (tmp == NULL) {
2387 int i_local;
2388 Do(i_local,curr_group)
2389 free_linked_list_2(pairs[i_local]);
2390 free_graph(mtxA_tmp);
2391 free_graph(precond);
2392 free_linked_list(l);
2393 unionfind_free();
2394 #ifdef USE_HEAPSORT
2395 free_heap(h);
2396 #endif
2397 taucs_free(diag);
2398 taucs_free(already_added);
2399 taucs_free(pi);
2400 taucs_free(groups);
2401 taucs_free(complete_subgraph);
2402 taucs_free(pairs);
2403 taucs_free(closed_cycle);
2404 return NULL;
2405 }
2406 tmp->group_1 = g1;
2407 tmp->group_2 = g2;
2408 tmp->no_edges = 0;
2409 tmp->next = pairs[(g1+g2)%curr_group];
2410 pairs[(g1+g2)%curr_group] = tmp;
2411 p = tmp;
2412 }
2413
2414 /* first check if the next edges to be added are inner to one of the subgraphs
2415 (not cross edges) */
2416 if (p->no_edges < 2)
2417 {
2418 #ifdef USE_HEAPSORT
2419 weight = mtxA_tmp->edges[Ah.edges[i]].v;
2420 #else
2421 weight = mtxA_tmp->edges[i].v;
2422 #endif
2423
2424 edge_sign = (weight>0);
2425
2426 x = find_set(u);
2427 y = find_set(v);
2428 closed_cycle_x = closed_cycle[x];
2429 closed_cycle_y = closed_cycle[y];
2430
2431 if (complete_subgraph[g1].completed_to_basis)
2432 if (fabs(weight)<=fabs(complete_subgraph[g1].c))
2433 if ((p->no_edges==0)||
2434 ((groups[p->a[0]]!=g1)||
2435 (groups[p->b[0]]!=g1)))
2436 {
2437 p->a[p->no_edges] = complete_subgraph[g1].a;
2438 p->b[p->no_edges] = complete_subgraph[g1].b;
2439 p->c[p->no_edges] = complete_subgraph[g1].c;
2440 p->cross[p->no_edges] = 0;
2441 (p->no_edges)++;
2442 }
2443
2444 if (p->no_edges < 2)
2445 if (complete_subgraph[g2].completed_to_basis)
2446 if (fabs(weight)<=fabs(complete_subgraph[g2].c))
2447 if ((p->no_edges==0)||
2448 ((groups[p->a[0]]!=g2)||
2449 (groups[p->b[0]]!=g2)))
2450 {
2451 p->a[p->no_edges] = complete_subgraph[g2].a;
2452 p->b[p->no_edges] = complete_subgraph[g2].b;
2453 p->c[p->no_edges] = complete_subgraph[g2].c;
2454 p->cross[p->no_edges] = 0;
2455 (p->no_edges)++;
2456 }
2457
2458 /* if p->no_edges == 2 get out of if */
2459
2460 if (p->no_edges == 1)
2461 {
2462 if (groups[p->a[0]] == groups[p->b[0]])
2463 {
2464 if (groups[p->a[0]] == groups[u])
2465 closed_cycle_x = 1;
2466 else
2467 closed_cycle_y = 1;
2468 }
2469
2470 }
2471
2472
2473 if (p->no_edges == 0) {
2474 /* sivan: added this brace below to avoid warning. */
2475 /* I have no idea why there are two identical ifs. */
2476 /* I hope I added the matching brace in the right place */
2477 if (p->no_edges == 0) {
2478 if (x != y) {
2479 if (!((closed_cycle[x])&&(closed_cycle[y]))) {
2480 #ifdef USE_HEAPSORT
2481 if (already_added[Ah.edges[i]]==0)
2482 #else
2483 if (already_added[i]==0)
2484 #endif
2485 {
2486 /*
2487 ivec1_p[Bent] = u;
2488 ivec2_p[Bent] = v;
2489 dvec_p[Bent] = weight;
2490 */
2491 precond->edges[Bent].i = u;
2492 precond->edges[Bent].j = v;
2493 precond->edges[Bent].v = weight;
2494
2495 Bent++;
2496
2497 diag[u] += fabs(weight);
2498 diag[v]+=fabs(weight);
2499 }
2500 p->a[p->no_edges] = u;
2501 p->b[p->no_edges] = v;
2502 p->c[p->no_edges] = weight;
2503 p->cross[p->no_edges] = 1;
2504 (p->no_edges)++;
2505 }
2506 } else {
2507 if ((edge_sign != (label[u]^label[v])) && (closed_cycle[x]==0)) {
2508 #ifdef USE_HEAPSORT
2509 if (already_added[Ah.edges[i]]==0)
2510 #else
2511 if (already_added[i]==0)
2512 #endif
2513 {
2514 /*
2515 ivec1_p[Bent] = u;
2516 ivec2_p[Bent] = v;
2517 dvec_p[Bent] = weight;
2518 */
2519 precond->edges[Bent].i = u;
2520 precond->edges[Bent].j = v;
2521 precond->edges[Bent].v = weight;
2522
2523 Bent++;
2524
2525 diag[u] += fabs(weight);
2526 diag[v] += fabs(weight);
2527 }
2528 p->a[p->no_edges] = u;
2529 p->b[p->no_edges] = v;
2530 p->c[p->no_edges] = weight;
2531 p->cross[p->no_edges] = 1;
2532 (p->no_edges)++;
2533 }
2534 }
2535 }
2536 }
2537
2538 if (p->no_edges == 1) {
2539 if (p->cross[0]==0) {
2540 if (groups[p->a[0]]==groups[u])
2541 closed_cycle_x = 1;
2542 else
2543 closed_cycle_y = 1;
2544 }
2545
2546 if ((!((closed_cycle_x)&&(closed_cycle_y)))&&(p->cross[0]==0))
2547 {
2548 #ifdef USE_HEAPSORT
2549 if (already_added[Ah.edges[i]]==0)
2550 #else
2551 if (already_added[i]==0)
2552 #endif
2553 {
2554 /*
2555 ivec1_p[Bent] = u;
2556 ivec2_p[Bent] = v;
2557 dvec_p[Bent] = weight;
2558 */
2559 precond->edges[Bent].i = u;
2560 precond->edges[Bent].j = v;
2561 precond->edges[Bent].v = weight;
2562
2563 Bent++;
2564
2565 diag[u] += fabs(weight);
2566 diag[v]+=fabs(weight);
2567 }
2568 p->a[p->no_edges] = u;
2569 p->b[p->no_edges] = v;
2570 p->c[p->no_edges] = weight;
2571 p->cross[p->no_edges] = 1;
2572 (p->no_edges)++;
2573 }
2574
2575 if ((!closed_cycle_x)&&(!closed_cycle_y)&&(p->cross[0]==1))
2576 {
2577 int x1,y1,edge_sign1;
2578 x1 = find_set(p->a[0]);
2579 y1 = find_set(p->b[0]);
2580 edge_sign1 = (p->c[0]>0);
2581
2582 if ((edge_sign^edge_sign1^label[u]^label[v]^label[p->a[0]]^label[p->b[0]]) ==1)
2583 {
2584 #ifdef USE_HEAPSORT
2585 if (already_added[Ah.edges[i]]==0)
2586 #else
2587 if (already_added[i]==0)
2588 #endif
2589 {
2590 /*
2591 ivec1_p[Bent] = u;
2592 ivec2_p[Bent] = v;
2593 dvec_p[Bent] = weight;
2594 */
2595
2596 precond->edges[Bent].i = u;
2597 precond->edges[Bent].j = v;
2598 precond->edges[Bent].v = weight;
2599
2600 Bent++;
2601
2602 diag[u] += fabs(weight);
2603 diag[v]+=fabs(weight);
2604 }
2605 p->a[p->no_edges] = u;
2606 p->b[p->no_edges] = v;
2607 p->c[p->no_edges] = weight;
2608 p->cross[p->no_edges] = 1;
2609 (p->no_edges)++;
2610 }
2611 }
2612 }
2613 }
2614 }
2615 }
2616
2617 Do(i,curr_group)
2618 free_linked_list_2(pairs[i]);
2619
2620 taucs_free(pairs);
2621 #ifdef USE_HEAPSORT
2622 /* free_heap(h);*/
2623 free_heap(Ah);
2624 free_heap(Bh);
2625 #endif
2626 taucs_free(closed_cycle);
2627 unionfind_free();
2628
2629 precond->nent = Bent;
2630
2631 /*** ALLOCATED: mtxA_tmp,diag,already_added,pi,precond,linked ***/
2632 /*** ALLOCATED: groups,complete_subgraph ***/
2633 }
2634
2635 taucs_free(complete_subgraph);
2636 taucs_free(already_added);
2637
2638 wtime_pair_bases = taucs_wtime() - wtime_pair_bases;
2639 taucs_printf("\t\tAMWB pair bases = %.3f seconds\n",wtime_pair_bases);
2640
2641
2642 /*** ALLOCATED: mtxA_tmp,diag,pi,precond,linked ***/
2643 /*** ALLOCATED: groups ***/
2644
2645 /* allocate more memory to the preconditioner if needed */
2646
2647 wtime = taucs_wtime();
2648 if (precond_maxnent < Bent + n)
2649 {
2650 taucs_printf("adding space for %d entries in precond for diagonal entries\n",Bent+n-precond_maxnent);
2651 precond_maxnent = Bent + n;
2652
2653
2654 if (graph_resize(precond,precond_maxnent) == -1) {
2655 /* precond has not been freed */
2656 free_graph(mtxA_tmp);
2657 free_graph(precond);
2658 free_linked_list(l);
2659 taucs_free(diag);
2660 taucs_free(pi);
2661 taucs_free(groups);
2662 return NULL;
2663 }
2664
2665 precond->nent = Bent;
2666 /*
2667 ivec1_p = precond->ivec1 ;
2668 ivec2_p = precond->ivec2 ;
2669 dvec_p = precond->dvec ;
2670 */
2671 }
2672
2673 Do(i,n)
2674 {
2675 /*
2676 ivec1_p[Bent] = i;
2677 ivec2_p[Bent] = i;
2678 dvec_p[Bent] = diag[i];
2679 */
2680
2681 precond->edges[Bent].i = i;
2682 precond->edges[Bent].j = i;
2683 precond->edges[Bent].v = diag[i];
2684
2685 Bent++;
2686 }
2687 precond->nent = Bent;
2688
2689 wtime = taucs_wtime() - wtime;
2690 taucs_printf("\t\tAMWB precond resize = %.3f seconds\n",wtime);
2691
2692
2693 taucs_printf("actual number of entries in preconditioner = %d\n",Bent);fflush(stdout);
2694
2695 taucs_free(diag);
2696 taucs_free(groups);
2697 taucs_free(pi);
2698 free_linked_list(l);
2699 free_graph(mtxA_tmp);
2700
2701 /*** ALLOCATED: precond ***/
2702
2703 wtime = taucs_wtime();
2704 out = graph_to_ccs_matrix(precond);
2705 if (!out) {
2706 free_graph(precond);
2707 return NULL;
2708 }
2709 wtime = taucs_wtime() - wtime;
2710 taucs_printf("\t\tAMWB graph-to-matrix = %.3f seconds\n",wtime);
2711
2712 free_graph(precond);
2713
2714
2715 wtime_total = taucs_wtime() - wtime_total;
2716 taucs_printf("\t\tAMWB time = %.3f seconds (%.3f sort)\n",
2717 wtime_total,
2718 wtime_sort);
2719
2720 return out;
2721 }
2722
2723 /*********************************************************/
2724 /* MST-specific routines */
2725 /*********************************************************/
2726
2727 typedef struct msthea {
2728 int heap_size;
2729 int* vertices;
2730 double* key;
2731 } mstheap;
2732
2733 #define INF 100000000.0
2734
2735 #define Parent(i) ((((i)+1)/2) - 1)
2736 #define Left(i) ((((i)+1)*2) - 1)
2737 #define Right(i) ((((i)+1)*2 + 1) - 1)
2738
2739 static
mstheap_exchange(mstheap A,int a,int b,int * point_to_heap)2740 void mstheap_exchange(mstheap A,int a,int b,int *point_to_heap)
2741 {
2742 int tmp1;
2743 double tmp2;
2744
2745 tmp1 = A.vertices[a];
2746 A.vertices[a] = A.vertices[b];
2747 A.vertices[b] = tmp1;
2748
2749 tmp2 = A.key[a];
2750 A.key[a] = A.key[b];
2751 A.key[b] = tmp2;
2752
2753 point_to_heap[A.vertices[a]] = a;
2754 point_to_heap[A.vertices[b]] = b;
2755 }
2756
2757 static
Mstheapify(mstheap A,int i,int * point_to_heap)2758 void Mstheapify(mstheap A,int i,int *point_to_heap)
2759 {
2760 int l,r,largest;
2761
2762 l = Left(i);
2763 r = Right(i);
2764
2765 if ((l < A.heap_size) && (A.key[l] > A.key[i]))
2766 largest = l;
2767 else
2768 largest = i;
2769
2770 if ((r < A.heap_size) && (A.key[r] > A.key[largest]))
2771 largest = r;
2772
2773 if (largest != i)
2774 {
2775 mstheap_exchange(A,i,largest,point_to_heap);
2776 Mstheapify(A,largest,point_to_heap);
2777 }
2778 }
2779
2780 static
build_mstheap(int r,int size,mstheap * h,int * point_to_heap,char alloc_flag)2781 int build_mstheap(int r,int size,mstheap *h,int *point_to_heap,char alloc_flag)
2782 {
2783 int i;
2784
2785 h->heap_size = size;
2786
2787 if (alloc_flag)
2788 {
2789 h->vertices = (int *)taucs_malloc(size * sizeof(int));
2790 h->key = (double *)taucs_malloc(size * sizeof(double));
2791
2792 if (h->key == NULL || h->vertices == NULL) return -1;
2793 }
2794
2795 Do(i,size)
2796 {
2797 h->vertices[i] = i;
2798 h->key[i] = (- INF);
2799 point_to_heap[i] = i;
2800 }
2801
2802 h->key[r] = 0;
2803
2804 mstheap_exchange((*h),r,0,point_to_heap);
2805
2806 return 0;
2807 }
2808
2809 static
Mstheap_Extract_Max(mstheap * A,int * point_to_heap)2810 int Mstheap_Extract_Max(mstheap *A,int *point_to_heap)
2811 {
2812 int out;
2813
2814 assert (A->heap_size >= 1);
2815
2816 out = A->vertices[0];
2817 A->vertices[0] = A->vertices[A->heap_size - 1];
2818 A->key[0] = A->key[A->heap_size - 1];
2819 point_to_heap[A->vertices[0]] = 0;
2820
2821 A->heap_size --;
2822
2823 Mstheapify((*A),0,point_to_heap);
2824
2825 return(out);
2826 }
2827
2828 static
mstheap_increase_key(mstheap h,int vv,double val,int * point_to_heap)2829 void mstheap_increase_key(mstheap h,int vv,double val,int *point_to_heap)
2830 {
2831 int i;
2832 /*int count=0;*/
2833 double key;
2834 int ver, v;
2835
2836 v = point_to_heap[vv];
2837
2838 h.key[v] = val;
2839
2840 key = h.key[v];
2841 ver = h.vertices[v];
2842
2843 i = v;
2844
2845 while((i>0) && (h.key[Parent(i)] < key))
2846 {
2847 h.key[i] = h.key[Parent(i)];
2848 h.vertices[i] = h.vertices[Parent(i)];
2849 point_to_heap[h.vertices[i]] = i;
2850 i = Parent(i);
2851 }
2852
2853 h.key[i] = key;
2854 h.vertices[i] = ver;
2855 point_to_heap[h.vertices[i]] = i;
2856
2857 }
2858
2859 static
free_mstheap(mstheap h)2860 void free_mstheap(mstheap h)
2861 {
2862 taucs_free(h.vertices);
2863 taucs_free(h.key);
2864 }
2865
2866 static
add_heavy_edges(graph * mtxA,graph * precond,int Bent,edge ** array,int * groups,int no_groups,int * pi,double * diag)2867 int add_heavy_edges(graph* mtxA,
2868 graph* precond,
2869 int Bent,
2870 edge **array,
2871 int *groups,
2872 int no_groups,
2873 int *pi,
2874 double *diag)
2875 {
2876 three **pairs;
2877 int i,j,a,b,k;
2878 double w;
2879 /*
2880 int *ivec1_p, *ivec2_p ;
2881 double *dvec_p;
2882 int *ivec1, *ivec2 ;
2883 double *dvec;
2884 */
2885 int orig_Bent;
2886 three *p,*tmp;
2887 int n,nent;
2888 int precond_maxnent,chunk;
2889
2890 three *pool;
2891 int next_in_pool;
2892
2893 n = mtxA->n;
2894 nent = mtxA->nent;
2895
2896 precond_maxnent = precond->max_size;
2897 chunk = max(min((no_groups*(no_groups-1)/2)/10,5000),10000);
2898 orig_Bent = Bent;
2899
2900 pairs = (three **) taucs_calloc(no_groups,sizeof(three *));
2901 pool = (three*) taucs_malloc((n+nent) * sizeof(three));
2902 if (!pairs || !pool) {
2903 taucs_free(pairs);
2904 taucs_free(pool);
2905 return -1;
2906 }
2907 next_in_pool = 0;
2908
2909 /*
2910 ivec1 = mtxA->ivec1;
2911 ivec2 = mtxA->ivec2;
2912 dvec = mtxA->dvec;
2913 */
2914
2915 Do(k,n) {
2916 i = k;
2917 j = pi[k];
2918
2919 if (j != (-1)) {
2920 a = min(groups[i],groups[j]);
2921 b = max(groups[i],groups[j]);
2922
2923 if (a != b) {
2924 p = pairs[(a+b)%no_groups];
2925 while (p!=NULL) {
2926 if ((p->group_1 == a)&&(p->group_2 == b))
2927 goto after;
2928 p = p->next;
2929 }
2930 after:
2931 if (p == NULL) {
2932 /*tmp = (three *)taucs_malloc(sizeof(three));*/
2933 /*if (tmp == NULL) {taucs_printf("ERROR! OUT OF MEMORY\n");exit(234);}*/
2934 assert(next_in_pool < n+nent);
2935 tmp = pool + next_in_pool;
2936 next_in_pool ++;
2937 tmp->group_1 = a;
2938 tmp->group_2 = b;
2939 tmp->already_connected = 1;
2940 tmp->next = pairs[(a+b)%no_groups];
2941 pairs[(a+b)%no_groups] = tmp;
2942 }
2943 }
2944 }
2945 }
2946
2947 Do(k,nent) {
2948 i = (mtxA->edges)[k].i;
2949 j = (mtxA->edges)[k].j;
2950
2951 a = min(groups[i],groups[j]);
2952 b = max(groups[i],groups[j]);
2953
2954 if (a != b) {
2955 w = - (mtxA->edges)[k].v;
2956 if (w) {
2957 p = pairs[(a+b)%no_groups];
2958 while (p!=NULL) {
2959 if ((p->group_1 == a)&&(p->group_2 == b)) {
2960 if (p->already_connected==0)
2961 if (w > p->c) {
2962 p->a = i;
2963 p->b = j;
2964 p->c = w;
2965 }
2966 goto after1;
2967 }
2968 p = p->next;
2969 }
2970 after1:
2971 if (p == NULL) {
2972 /*tmp = (three *)taucs_malloc(sizeof(three));*/
2973 /*if (tmp == NULL) {taucs_printf("ERROR! OUT OF MEMORY\n");exit(234);}*/
2974 assert(next_in_pool < n+nent);
2975 tmp = pool + next_in_pool;
2976 next_in_pool ++;
2977 tmp->group_1 = a;
2978 tmp->group_2 = b;
2979 tmp->a = i;
2980 tmp->b = j;
2981 tmp->c = w;
2982 tmp->already_connected = 0;
2983 tmp->next = pairs[(a+b)%no_groups];
2984 pairs[(a+b)%no_groups] = tmp;
2985 }
2986 }
2987 }
2988 }
2989
2990 /*
2991 ivec1_p = precond->ivec1;
2992 ivec2_p = precond->ivec2;
2993 dvec_p = precond->dvec;
2994 */
2995
2996 Do(i,no_groups) {
2997 p = pairs[i];
2998 while(p!=NULL) {
2999 if(p->already_connected == 0) {
3000 if (p->a > p->b)
3001 swap(p->a,p->b);
3002
3003 (precond->edges)[Bent].i = p->a;
3004 (precond->edges)[Bent].j = p->b;
3005 (precond->edges)[Bent].v = -(p->c);
3006
3007 Bent++;
3008
3009 if (Bent == precond_maxnent) {
3010 precond_maxnent += chunk;
3011 taucs_printf("adding space for %d entries in precond\n",chunk);
3012 graph_resize(precond,precond_maxnent);
3013 precond->nent = orig_Bent;
3014
3015 /*
3016 ivec1_p = precond->ivec1;
3017 ivec2_p = precond->ivec2;
3018 dvec_p = precond->dvec;
3019 */
3020
3021 }
3022 diag[p->a] -= (-p->c);
3023 diag[p->b] -= (-p->c);
3024 }
3025
3026 p = p->next;
3027 }
3028 }
3029
3030 /*Do(i,no_groups) free_linked_list_2(pairs[i]);*/
3031 taucs_free(pool);
3032 taucs_free(pairs);
3033
3034 return(Bent);
3035 }
3036
3037 static int Dijkstra(graph *mtxA,int r,int *pi,linked *l,int *d,int *maxdist,int *partition,double min,double y,int j);
3038
3039
Av_Part_W(graph * mtxA,int * partition,int * new_partition,int * parts,graph * out,int nparts)3040 static int Av_Part_W(graph *mtxA,int *partition,int *new_partition,int *parts,graph *out,int nparts)
3041 {
3042 /* Peleg, Noga et al.'s partition. From Peleg's book, page 217 */
3043
3044 linked *l = NULL, *l_c = NULL;
3045 int *pi=NULL,*d=NULL,i,k,*findrho=NULL,minrho,maxdist,classes,n,nent,root,j,curr_partition=0;
3046 int row, col;
3047 int *pi1 = 0; /* warning */
3048 double x, y, min, max, not;
3049 byte bool=1;
3050 edge *p,*dummy, *pe ,*max_pe;
3051 int count = 0;
3052
3053 n = mtxA->n;
3054 nent = mtxA->nent;
3055
3056 x = exp(sqrt(log(n)*log(log(n))))/3;
3057 y = x * 9*log(n) * (floor(3*log(n)/log(x))+1);
3058
3059 pi = (int *)taucs_malloc(n*sizeof(int));
3060 d = (int *)taucs_malloc(n*sizeof(int));
3061 l = create_linked_list(mtxA,n,mtxA->nent,&min,&max);
3062 if (!pi || !d || !l)
3063 goto exit_Av_Part_W;
3064
3065 classes = (int)(log(max/min)/log(y))+1;
3066
3067 for(i=0;i<n;i++)
3068 new_partition[i] = -1;
3069
3070 j = 1;
3071
3072 while (count < n)
3073 {
3074 root = rand() % n;
3075 if (new_partition[root] == -1)
3076 {
3077 for(i=0;i<n;i++)
3078 d[i] = -1;
3079 Dijkstra(mtxA,root,pi,l,d,&maxdist,partition,min,y,j);
3080
3081 findrho = (int *)taucs_calloc((maxdist+1)*classes,sizeof(int));
3082 if (!findrho)
3083 goto exit_Av_Part_W;
3084
3085 for(i=0;i<n;i++)
3086 {
3087 if (d[i] != -1)
3088 {
3089 p = (l->point)[i];
3090 while (p != NULL)
3091 {
3092 if ((d[mtxA->edges[p->entry_no].i]!=-1) && (d[mtxA->edges[p->entry_no].j]!=-1))
3093 if ((d[mtxA->edges[p->entry_no].i] - d[mtxA->edges[p->entry_no].j] >= -1) &&
3094 (d[mtxA->edges[p->entry_no].i] - d[mtxA->edges[p->entry_no].j] <= 1))
3095 findrho[(max(d[mtxA->edges[p->entry_no].i],d[mtxA->edges[p->entry_no].j]))*classes+
3096 /*(int)(log(abs(mtxA->edges[p->entry_no].v)/min)/log(y))]++; omer*/
3097 (int)(log(abs((int)(mtxA->edges[p->entry_no].v))/min)/log(y))]++;
3098 p = p->next;
3099 }
3100 }
3101 }
3102
3103 for(i=0;i<min(j,classes);i++)
3104 findrho[i] = 0; /* ignore edges connecting two vertices whose distance is 0 */
3105
3106 /* At this point, findrho[k,i], or findrho[k*classes+i], contains the number of edges in E_i,
3107 connecting two vertices whose distance is k, or a vertex of distance k with a vertex of distance k-1 */
3108
3109 for(k=1;k<maxdist;k++)
3110 for(i=0;i<min(j,classes);i++)
3111 findrho[k*classes+i] += findrho[(k-1)*classes+i];
3112
3113 /* At this point, findrho[k,i], or findrho[k*classes+i], contains the number of edges in E_i,
3114 connecting two vertices whose distance is j, or a vertex of distance j with a vertex of distance j-1
3115 for j=1,...,k */
3116
3117 for(minrho=1;minrho<maxdist;minrho++)
3118 {
3119 bool = 1;
3120 for(k=0;k<min(j,classes);k++)
3121 {
3122 if ((double)(findrho[(minrho+1)*classes+k]-findrho[minrho*classes+k]) > (findrho[minrho*classes+k])/x)
3123 bool = 0;
3124 }
3125 if (bool)
3126 goto afterr;
3127 }
3128
3129 afterr:
3130 if (bool)
3131 {
3132 for(i=0;i<n;i++)
3133 if ((d[i] <= minrho) && (d[i] != -1) )
3134 if (new_partition[i] == -1)
3135 {
3136 count ++;
3137 new_partition[i] = curr_partition;
3138 }
3139 }
3140 else
3141 {
3142 for(i=0;i<n;i++)
3143 if ((new_partition[i] == -1) && (d[i] != -1))
3144 {
3145 count ++;
3146 new_partition[i] = curr_partition;
3147 }
3148 }
3149
3150 for(i=0;i<n;i++)
3151 {
3152 if (new_partition[i] == curr_partition)
3153 l->point[i] = NULL;
3154 else
3155 {
3156 p = l->point[i];
3157 l->point[i] = NULL;
3158 while (p != NULL)
3159 {
3160 if ((new_partition[mtxA->edges[p->entry_no].i] != curr_partition) && (new_partition[mtxA->edges[p->entry_no].j] != curr_partition))
3161 {
3162 dummy = l->point[i];
3163 l->point[i] = p;
3164 p = p->next;
3165 (l->point[i])->next = dummy;
3166 }
3167 else
3168 p = p->next;
3169 }
3170 }
3171 }
3172
3173 curr_partition ++;
3174 j++;
3175 taucs_free(findrho);
3176 }
3177 }
3178
3179 *parts = curr_partition;
3180
3181 l_c = create_linked_list_cluster(mtxA,nparts,mtxA->nent,¬,¬,partition,new_partition);
3182 pi1 = (int *)taucs_malloc(nparts*sizeof(int));
3183
3184 if (!l_c || !pi1)
3185 goto exit_Av_Part_W;
3186
3187 for(i=0;i<n;i++)
3188 if (pi[i] != -1)
3189 {
3190 if (partition[i] != partition[pi[i]])
3191 pi1[partition[i]] = partition[pi[i]];
3192 }
3193 else
3194 pi1[partition[i]] = -1;
3195
3196 Do(i,nparts) {
3197 row = i;
3198 col = pi1[i];
3199 if (col != (-1)) {
3200 pe = l_c->point[row];
3201 max_pe = NULL;
3202 while (pe != NULL) {
3203 if ( (partition[(mtxA->edges)[pe->entry_no].j] == col)
3204 || (partition[(mtxA->edges)[pe->entry_no].i] == col))
3205 {
3206 if (!max_pe)
3207 max_pe = pe;
3208 else
3209 if (-mtxA->edges[pe->entry_no].v > -mtxA->edges[max_pe->entry_no].v)
3210 max_pe = pe;
3211 }
3212 pe = pe->next;
3213 }
3214
3215 assert(max_pe);
3216 out->edges[out->nent].i = (mtxA->edges)[max_pe->entry_no].i;
3217 out->edges[out->nent].j = (mtxA->edges)[max_pe->entry_no].j;
3218 out->edges[out->nent].v = (mtxA->edges)[max_pe->entry_no].v;
3219 out->nent++;
3220
3221 }
3222 }
3223
3224
3225 taucs_free(pi);taucs_free(d);free_linked_list(l);free_linked_list(l_c);taucs_free(pi1);
3226 return 1;
3227
3228 exit_Av_Part_W:
3229 taucs_free(pi);taucs_free(d);free_linked_list(l);taucs_free(l);taucs_free(findrho);free_linked_list(l_c);taucs_free(l_c);taucs_free(pi1);
3230 return 0;
3231
3232
3233 }
3234
3235 static taucs_ccs_matrix *amst_preconditioner_create(graph *mtxA, double* diag,int rnd,double subgraphs,int stretch_flag);
3236 static int Prim(graph *mtxA,int r,int *pi,int *d,linked *l);
3237 /*static int Prim_cluster(graph *mtxA,int r,int *pi,linked *l,int *partition,int *new_partition,int nparts,int *point_to_heap,char *in_Q,mstheap h);*/
3238
Dijkstra(graph * mtxA,int r,int * pi,linked * l,int * d,int * maxdist,int * partition,double min,double y,int j)3239 static int Dijkstra(graph *mtxA,int r,int *pi,linked *l,int *d,int *maxdist,int *partition,double min,double y,int j)
3240 {
3241 /* Dijkstra is used in order to compute the distance of vertices from the root.
3242 The distance of two vertices within the same partition is 0 */
3243
3244 int n,entry,u,v,i;
3245 mstheap h;
3246 int *point_to_heap = NULL;
3247 char *in_Q = NULL;
3248 int size_of_Q;
3249 edge *p;
3250 double weight;
3251
3252 *maxdist=0;
3253
3254 n = mtxA->n;
3255
3256 point_to_heap = (int*) taucs_malloc(n*sizeof(int));
3257 in_Q = (char*) taucs_malloc(n*sizeof(char));
3258
3259 Do(i,n) in_Q[i] = 1;
3260 size_of_Q = n;
3261
3262 /* Prim's Algorithm - Introduction to Algorithms page 509 */
3263
3264 if (build_mstheap(r,n,&h,point_to_heap,1) == -1)
3265 goto exit_Dijkstra;
3266
3267 pi[r] = -1;
3268
3269 while(size_of_Q > 0) {
3270 if (h.key[0] == -INF)
3271 goto after_Dijkstra;
3272 assert(h.key[0] != -INF);
3273 weight = -h.key[0];
3274 u = Mstheap_Extract_Max(&h,point_to_heap);
3275 d[u] = (int)weight;
3276 if (weight>*maxdist)
3277 *maxdist = (int)weight;
3278
3279 in_Q[u] = 0;
3280 size_of_Q --;
3281
3282 p = (l->point)[u];
3283 while (p != NULL) {
3284 entry = p->entry_no;
3285 /*if ((int)(log(abs(mtxA->edges[entry].v/min))/log(y)) < j) omer*/
3286 if ((int)(log(abs((int)(mtxA->edges[entry].v/min)))/log(y)) < j)
3287 {
3288 /* v belongs to Adj[u] */
3289
3290 if ((mtxA->edges)[entry].j != u)
3291 v = (mtxA->edges)[entry].j;
3292 else
3293 v = (mtxA->edges)[entry].i;
3294
3295 if (in_Q[v] && (-h.key[point_to_heap[v]] > ((partition[v]==partition[u])?0:1) + weight))
3296 {
3297 pi[v] = u;
3298 mstheap_increase_key(h,v,-(((partition[v]==partition[u])?0:1) + weight),point_to_heap);
3299 }
3300 }
3301 p = p->next;
3302 }
3303 }
3304
3305 after_Dijkstra:
3306 free_mstheap(h);
3307 taucs_free(point_to_heap);
3308 taucs_free(in_Q);
3309 return 1;
3310
3311 exit_Dijkstra:
3312 free_mstheap(h);
3313 taucs_free(point_to_heap);
3314 taucs_free(in_Q);
3315 return 0;
3316
3317 }
3318
stupid_part(int * partition,int n,int j,int * nparts)3319 void stupid_part(int *partition,int n,int j,int *nparts)
3320 {
3321 int i,k,q;
3322
3323 k = 1<<j;
3324
3325 q = ((n%k == 0)?(n/k):(n/k+1));
3326
3327 for(i=0;i<n;i++)
3328 for(j=0;j<n;j++)
3329 {
3330 partition[i*n+j] = q*(i/k)+(j/k);
3331 }
3332 *nparts=partition[n*n-1]+1;
3333
3334 /* for(i=0;i<n*n;i++) */
3335 /* printf("ASD %d %d - %d %d\n",i,partition[i],n,k); */
3336 /* exit(345); */
3337 }
3338
low_stretch(graph * mtxA)3339 graph *low_stretch(graph *mtxA)
3340 {
3341 int *partition = NULL,*new_partition = NULL,*choose_root = NULL;
3342 int i,n,nparts,new_nparts;/* r omer*/
3343 /*double dummy; omer*/
3344 int *pi = NULL;
3345 /*linked *l = NULL;*/
3346 /*int k=0,j=0;*/ /* row,col, omer*/
3347 /*edge *pe,*max_pe; omer*/
3348 graph *out;
3349 /*int *point_to_heap = NULL;*/
3350 /*char *in_Q = NULL;*/
3351 /*mstheap h; omer*/
3352
3353 n = mtxA->n;
3354
3355 out = construct_graph(n-1);
3356 partition = (int *)taucs_malloc(n*sizeof(int));
3357 new_partition = (int *)taucs_malloc(n*sizeof(int));
3358 choose_root = (int *)taucs_malloc(n*sizeof(int));
3359 pi = (int *)taucs_malloc(n*sizeof(int));
3360
3361 if (!out || !partition || !new_partition || !choose_root || !pi)
3362 goto exit_low_stretch;
3363
3364 out->n = n;
3365 out->nent = 0;
3366
3367 for (i=0;i<n;i++)
3368 partition[i] = i;
3369
3370 nparts = n;
3371
3372 while (nparts > 1)
3373 {
3374 if (Av_Part_W(mtxA,partition,new_partition,&new_nparts,out,nparts) == 0)
3375 goto exit_low_stretch;
3376
3377 #if 0
3378 for(i=0;i<n;i++)
3379 pi[i] = -1;
3380
3381 out->n = n;
3382 out->nent = n-1;
3383
3384 j++;
3385
3386 /* stupid_part(new_partition,(int)(sqrt(n)),j,&new_nparts); */
3387
3388 for(i=0;i<new_nparts;i++)
3389 choose_root[i] = 0;
3390
3391 l = create_linked_list_cluster(mtxA,nparts,mtxA->nent,&dummy,&dummy,partition,new_partition);
3392 /* This linked list contains all the edges which connect vertices whose endpoints are
3393 in different sections in partition, but in the same section in new_partition.
3394 This will help us build trees within each section in new_partition. */
3395
3396 if (l == NULL)
3397 goto exit_low_stretch;
3398
3399
3400 point_to_heap = (int *)taucs_malloc(n*sizeof(int));
3401 in_Q = (char *)taucs_malloc(n*sizeof(char));
3402 h.vertices = NULL;
3403 h.key = NULL;
3404 h.vertices = (int *)taucs_malloc(n * sizeof(int));
3405 h.key = (double *)taucs_malloc(n * sizeof(double));
3406
3407 if (!point_to_heap || !in_Q || !h.vertices || !h.key)
3408 {
3409 taucs_free(point_to_heap);taucs_free(in_Q);taucs_free(h.vertices);taucs_free(h.key);
3410 goto exit_low_stretch;
3411 }
3412
3413 for(i=0;i<n;i++)
3414 {
3415 if (choose_root[new_partition[i]] == 0)
3416 {
3417 /* new_partition[i] is a section for which a tree has not yet been found */
3418 choose_root[new_partition[i]] = 1;
3419 r = partition[i];
3420 if (Prim_cluster(mtxA,r,pi,l,partition,new_partition,nparts,point_to_heap,in_Q,h) == 0)
3421 {
3422 free_linked_list(l);
3423 goto exit_low_stretch;
3424 }
3425 }
3426 }
3427
3428 taucs_free(point_to_heap);
3429 taucs_free(in_Q);
3430 free_mstheap(h);
3431
3432 Do(i,nparts) {
3433 row = i;
3434 col = pi[i];
3435 if (col != (-1)) {
3436 pe = l->point[row];
3437 max_pe = NULL;
3438 while (pe != NULL) {
3439 if ( (partition[(mtxA->edges)[pe->entry_no].j] == col)
3440 || (partition[(mtxA->edges)[pe->entry_no].i] == col))
3441 {
3442 if (!max_pe)
3443 max_pe = pe;
3444 else
3445 if (-mtxA->edges[pe->entry_no].v > -mtxA->edges[max_pe->entry_no].v)
3446 max_pe = pe;
3447 }
3448 pe = pe->next;
3449 }
3450
3451 assert(max_pe);
3452 out->edges[k].i = (mtxA->edges)[max_pe->entry_no].i;
3453 out->edges[k].j = (mtxA->edges)[max_pe->entry_no].j;
3454 out->edges[k].v = (mtxA->edges)[max_pe->entry_no].v;
3455 k++;
3456
3457 }
3458 }
3459
3460 out->nent = k;
3461 free_linked_list(l);
3462 #endif
3463
3464 for(i=0;i<n;i++)
3465 partition[i] = new_partition[i];
3466
3467 nparts = new_nparts;
3468
3469 }
3470
3471 assert(out->nent==(n-1)); /* helps verify that out is a tree */
3472
3473 taucs_free(partition);taucs_free(new_partition);taucs_free(choose_root);taucs_free(pi);
3474 return out;
3475
3476 exit_low_stretch:
3477 free_graph(out);taucs_free(partition);taucs_free(new_partition);taucs_free(choose_root);taucs_free(pi);
3478 return 0;
3479 }
3480
3481 static
Prim(graph * mtxA,int r,int * pi,int * d,linked * l)3482 int Prim(graph *mtxA,int r,int *pi,int *d,linked *l)
3483 {
3484 /* Prim's Algorithm - Introduction to Algorithms page 509 */
3485
3486 int n,entry,u,v,i;
3487 mstheap h;
3488 int *point_to_heap;
3489 char *in_Q;
3490 int size_of_Q;
3491 edge *p;
3492
3493 n = mtxA->n;
3494
3495 point_to_heap = (int*) taucs_malloc(n*sizeof(int));
3496 in_Q = (char*) taucs_malloc(n*sizeof(char));
3497
3498 Do(i,n) in_Q[i] = 1;
3499 size_of_Q = n;
3500
3501 /* Prim's Algorithm - Introduction to Algorithms page 509 */
3502
3503 if ((build_mstheap(r,n,&h,point_to_heap,1) == -1)) {
3504 taucs_free(in_Q);
3505 taucs_free(pi);
3506 taucs_free(point_to_heap);
3507 /* free linked_list; */
3508 return 0;
3509 }
3510
3511 pi[r] = -1;
3512 d[r] = 0;
3513
3514 while(size_of_Q > 0) {
3515 u = Mstheap_Extract_Max(&h,point_to_heap);
3516 in_Q[u] = 0;
3517 size_of_Q --;
3518
3519 p = (l->point)[u];
3520 while (p != NULL) {
3521 entry = p->entry_no;
3522 /* v belongs to Adj[u] */
3523
3524 if ((mtxA->edges)[entry].j != u)
3525 v = (mtxA->edges)[entry].j;
3526 else
3527 v = (mtxA->edges)[entry].i;
3528
3529 if (in_Q[v] && ((-((mtxA->edges)[entry].v)) > h.key[point_to_heap[v]])) {
3530 pi[v] = u;
3531 d[v] = d[u] + 1;
3532 mstheap_increase_key(h,v,-((mtxA->edges)[entry].v),point_to_heap);
3533 }
3534 p = p->next;
3535 }
3536 }
3537
3538 free_mstheap(h);
3539 taucs_free(point_to_heap);
3540 taucs_free(in_Q);
3541 return 1;
3542 }
3543
3544 #if 0
3545 static
3546 int Prim_cluster(graph *mtxA,int r,int *pi,linked *l,int *partition,int *new_partition,int nparts,int *point_to_heap,char *in_Q,mstheap h)
3547 {
3548 /* Prim's Algorithm - Introduction to Algorithms page 509 */
3549
3550 int n,entry,u,v,i;
3551 /* mstheap h; */
3552 /* int *point_to_heap = NULL; */
3553 /* char *in_Q = NULL; */
3554 int size_of_Q;
3555 edge *p;
3556
3557 n = mtxA->n;
3558
3559 /* point_to_heap = (int*) taucs_malloc(nparts*sizeof(int)); */
3560 /* in_Q = (char*) taucs_malloc(nparts*sizeof(char)); */
3561 /* if (!point_to_heap || !in_Q) */
3562 /* goto exit_Prim_cluster; */
3563
3564 Do(i,nparts) in_Q[i] = 1;
3565 size_of_Q = nparts;
3566
3567 /* Prim's Algorithm - Introduction to Algorithms page 509 */
3568
3569 if (build_mstheap(r,nparts,&h,point_to_heap,0) == -1)
3570 goto exit_Prim_cluster;
3571
3572 pi[r] = -1;
3573
3574 while(size_of_Q > 0) {
3575 if (h.key[0] == -INF)
3576 goto after_Prim_cluster;
3577 u = Mstheap_Extract_Max(&h,point_to_heap);
3578
3579 in_Q[u] = 0;
3580 size_of_Q --;
3581
3582 p = (l->point)[u];
3583 while (p != NULL) {
3584 entry = p->entry_no;
3585 /* v belongs to Adj[u] */
3586
3587 if (partition[(mtxA->edges)[entry].j] != u)
3588 v = partition[(mtxA->edges)[entry].j];
3589 else
3590 v = partition[(mtxA->edges)[entry].i];
3591
3592 if (in_Q[v] && ((-((mtxA->edges)[entry].v)) > h.key[point_to_heap[v]])) {
3593 pi[v] = u;
3594 mstheap_increase_key(h,v,-((mtxA->edges)[entry].v),point_to_heap);
3595 }
3596 p = p->next;
3597 }
3598 }
3599
3600 after_Prim_cluster:
3601 /* free_mstheap(h); */
3602 /* taucs_free(point_to_heap); */
3603 /* taucs_free(in_Q); */
3604 return 1;
3605
3606 exit_Prim_cluster:
3607 /* free_mstheap(h); */
3608 /* taucs_free(point_to_heap); */
3609 /* taucs_free(in_Q); */
3610 return 0;
3611 }
3612 #endif /* 0, we don't need this routine */
3613
dist(int i,int j,double w,graph * mtxA,int * pi,int * d,linked * l,double * dilation,double * congestion)3614 static double dist(int i, int j, double w, graph *mtxA,int *pi,int *d,linked *l,double *dilation,double *congestion)
3615 {
3616 double out=0;
3617 int tmp;
3618 edge *e;
3619 double q = 0;
3620
3621 if (d[i] < d[j])
3622 {tmp = i;i=j;j=tmp;}
3623
3624 /* now we know that d[i] >= d[j] */
3625
3626 while (d[i] > d[j])
3627 {
3628 e = l->point[i];
3629 while ((mtxA->edges[e->entry_no].i != pi[i]) && (mtxA->edges[e->entry_no].j != pi[i]))
3630 e = e->next;
3631 out += mtxA->edges[e->entry_no].v;
3632 q += mtxA->edges[e->entry_no].v;
3633 congestion[i] += mtxA->edges[e->entry_no].v/w;
3634 i = pi[i];
3635 }
3636
3637 /* now we know that d[i] == d[j] */
3638
3639 while (i != j)
3640 {
3641 e = l->point[i];
3642 while ((mtxA->edges[e->entry_no].i != pi[i]) && (mtxA->edges[e->entry_no].j != pi[i]))
3643 e = e->next;
3644 out += mtxA->edges[e->entry_no].v;
3645 q += mtxA->edges[e->entry_no].v;
3646 congestion[i] += mtxA->edges[e->entry_no].v/w;
3647 i = pi[i];
3648
3649 e = l->point[j];
3650 while ((mtxA->edges[e->entry_no].i != pi[j]) && (mtxA->edges[e->entry_no].j != pi[j]))
3651 e = e->next;
3652 out += mtxA->edges[e->entry_no].v;
3653 q += mtxA->edges[e->entry_no].v;
3654 congestion[j] += mtxA->edges[e->entry_no].v/w;
3655 j = pi[j];
3656 }
3657
3658 *dilation = max(*dilation,q/w);
3659 return(out);
3660 }
3661
find_stretch(graph * mtxA,int * pi,int * d)3662 static double find_stretch(graph *mtxA,int *pi,int *d)
3663 {
3664 int i,E=0;
3665 double stretch=0,dummy;
3666 linked *l;
3667 double dilation;
3668 double *congestion,cong=0;
3669
3670 congestion = (double *)taucs_calloc(mtxA->n,sizeof(double));
3671
3672 l = create_linked_list(mtxA,mtxA->n,mtxA->nent,&dummy,&dummy);
3673 assert(l && congestion);
3674
3675 for (i=0;i<mtxA->nent;i++)
3676 {
3677 if (mtxA->edges[i].i != mtxA->edges[i].j)
3678 {
3679 E++;
3680 stretch += dist(mtxA->edges[i].i,mtxA->edges[i].j,mtxA->edges[i].v,mtxA,pi,d,l,&dilation,congestion)/mtxA->edges[i].v;
3681 }
3682 }
3683
3684 Do(i,mtxA->n)
3685 cong = max(cong,congestion[i]);
3686
3687 printf("Cong-Dil = %f\n",cong*dilation);
3688
3689 free_linked_list(l);
3690 taucs_free(congestion);
3691 return(stretch/E);
3692 }
3693
3694 static
3695 taucs_ccs_matrix*
amst_preconditioner_create(graph * mtxA,double * diag,int rnd,double subgraphs,int stretch_flag)3696 amst_preconditioner_create(graph *mtxA, double* diag,
3697 int rnd,
3698 double subgraphs,
3699 int stretch_flag)
3700 {
3701 taucs_ccs_matrix* out;
3702 /*
3703 sym_matrix *mtxA;
3704 sym_matrix *precond;
3705 */
3706 /*graph *mtxA;*/
3707 graph *precond;
3708 int i;
3709 /*int tmp;*/
3710 /*
3711 int *ivec1, *ivec2 ;
3712 int *ivec1_p, *ivec2_p ;
3713 double *dvec;
3714 double *dvec_p;
3715 */
3716 int *pi,*d;
3717 edge **array;
3718 edge *p;
3719 int r;
3720 int Bent,row,col;
3721 double weight = 0;
3722 int *first_child,*next_child;
3723 int *groups,*sub_tree_sizes;
3724 int curr_group;
3725 int n;
3726 int precond_maxnent,chunk;
3727 /*double *diag;*/
3728 /*linked l;*/
3729 linked* lp;
3730 double dummy;
3731 double wtime;
3732
3733 n = mtxA->n;
3734
3735 wtime = taucs_wtime();
3736
3737 pi = (int*) taucs_malloc(n*sizeof(int));
3738 d = (int*) taucs_malloc(n*sizeof(int));
3739
3740 /*l = create_linked_list_old(mtxA,n,mtxA->nent); */ /* THIS MAY RUN OUT OF MEMORY ! */
3741 lp = create_linked_list(mtxA,n,mtxA->nent,&dummy,&dummy); /* THIS MAY RUN OUT OF MEMORY ! */
3742
3743 if (!pi || !d) {
3744 taucs_free(pi);
3745 taucs_free(d);
3746 taucs_free(diag);
3747 free_graph(mtxA);
3748 }
3749
3750 #if 0
3751 taucs_free(diag);
3752 diag = analyze_graph(mtxA); /* should change! */
3753 #endif
3754
3755 /* array is an array of linked lists, which hold the
3756 off-diagonal entries of mtxA */
3757
3758 array = lp->point;
3759
3760 Do(i,n) pi[i] = -1;
3761
3762 /*
3763 ivec1 = mtxA->ivec1 ;
3764 ivec2 = mtxA->ivec2 ;
3765 dvec = mtxA->dvec ;
3766 */
3767
3768
3769 wtime = taucs_wtime() - wtime;
3770 taucs_printf("\t\tAMST prepare for mst = %.3f seconds\n",wtime);
3771
3772 wtime = taucs_wtime();
3773
3774 r = rnd % n;
3775
3776 if (stretch_flag)
3777 {
3778 graph *low_stretch_tree;
3779 linked *lp_low;
3780
3781 if ((low_stretch_tree = low_stretch(mtxA)) == NULL)
3782 {
3783 taucs_free(pi);
3784 taucs_free(diag);
3785 free_graph(mtxA);
3786 return 0;
3787 }
3788
3789 lp_low = create_linked_list(low_stretch_tree,n,low_stretch_tree->nent,&dummy,&dummy); /* THIS MAY RUN OUT OF MEMORY ! */
3790 if (lp_low == NULL)
3791 {
3792 taucs_free(pi);
3793 taucs_free(diag);
3794 free_graph(mtxA);
3795 free_graph(low_stretch_tree);
3796 return 0;
3797 }
3798 Do(i,n)
3799 pi[i] = -2;
3800 Prim(low_stretch_tree,r,pi,d,lp_low);
3801
3802 free_graph(low_stretch_tree);
3803 free_linked_list(lp_low);
3804 }
3805
3806 else
3807 {
3808 Prim(mtxA,r,pi,d,lp);
3809 }
3810
3811 /* pi now contains the parent array of the tree, d the distance from the root */
3812 printf("Stretch = %f\n",find_stretch(mtxA,pi,d));
3813
3814
3815 groups = (int *) taucs_malloc(n*sizeof(int));
3816 first_child = (int *) taucs_malloc(n*sizeof(int));
3817 next_child = (int *) taucs_malloc(n*sizeof(int));
3818 sub_tree_sizes = (int *) taucs_malloc(n*sizeof(int));
3819
3820 if (!groups || !first_child || !next_child || !sub_tree_sizes) {
3821 taucs_free(groups);
3822 taucs_free(first_child);
3823 taucs_free(next_child);
3824 taucs_free(sub_tree_sizes);
3825
3826 taucs_free(pi);
3827 taucs_free(diag);
3828 free_graph(mtxA);
3829 /* free linked_list; */
3830 }
3831
3832
3833 if (create_children_arrays(pi,n,&first_child,&next_child) == -1) {
3834 taucs_free(groups);
3835 taucs_free(first_child);
3836 taucs_free(next_child);
3837 taucs_free(sub_tree_sizes);
3838
3839 taucs_free(pi);
3840 taucs_free(diag);
3841 free_graph(mtxA);
3842 /* free linked_list; */
3843 }
3844
3845 wtime = taucs_wtime() - wtime;
3846 taucs_printf("\t\tAMST mst = %.3f seconds\n",wtime);
3847
3848 wtime = taucs_wtime();
3849
3850 Do(i,n) groups[i] = -1;
3851
3852 curr_group = 0;
3853
3854 compute_sub_tree_sizes(r,first_child,next_child,sub_tree_sizes);
3855
3856 /* now for every vertex v in the tree, sub_tree_sizes[v] is the size
3857 of the subtree whose root is v */
3858
3859 divide_to_groups(r,first_child,next_child,pi,&curr_group,groups,
3860 sub_tree_sizes,r,subgraphs,n);
3861 assign_group(r,curr_group,first_child,next_child,groups);
3862 curr_group++;
3863
3864 taucs_printf("actual number of subgraphs = %ld\n",curr_group);
3865
3866
3867 taucs_free(first_child);
3868 taucs_free(next_child);
3869 taucs_free(sub_tree_sizes);
3870
3871 wtime = taucs_wtime() - wtime;
3872 taucs_printf("\t\tAMST partition = %.3f seconds\n",wtime);
3873
3874 /* now the tree is devided into linked groups */
3875
3876 wtime = taucs_wtime();
3877
3878 chunk = max(min((curr_group*(curr_group-1)/2)/10,5000),100);
3879 precond_maxnent = (n-1) + n + chunk;
3880
3881 taucs_printf("allocating space for %d entries in precond\n",precond_maxnent);
3882
3883 precond = construct_graph(precond_maxnent);
3884 if (!precond) {
3885 taucs_free(pi);
3886 taucs_free(diag);
3887 free_graph(mtxA);
3888 /* free linked_list; */
3889 return NULL;
3890 }
3891 precond->n=mtxA->n;
3892
3893 /*
3894 ivec1_p = precond->ivec1 ;
3895 ivec2_p = precond->ivec2 ;
3896 dvec_p = precond->dvec ;
3897 */
3898
3899 Bent = 0;
3900
3901 wtime = taucs_wtime() - wtime;
3902 taucs_printf("\t\tAMST allocating mst precond = %.3f seconds\n",wtime);
3903
3904 wtime = taucs_wtime();
3905
3906 /* adds the tree edges to the preconditioner */
3907 Do(i,n) {
3908 row = i;
3909 col = pi[i];
3910 if (col != (-1)) {
3911 p = array[row];
3912
3913 while (p != NULL) {
3914 if ( ((mtxA->edges)[p->entry_no].j == col)
3915 || ((mtxA->edges)[p->entry_no].i == col)) {
3916 weight = (mtxA->edges)[p->entry_no].v;
3917 break;
3918 }
3919 p = p->next;
3920 }
3921
3922 (precond->edges)[Bent].i = row;
3923 (precond->edges)[Bent].j = col;
3924 (precond->edges)[Bent].v = weight;
3925
3926 Bent++;
3927
3928 diag[row] -= weight;
3929 diag[col] -= weight;
3930 }
3931 }
3932
3933 precond->nent = Bent;
3934
3935
3936 wtime = taucs_wtime() - wtime;
3937 taucs_printf("\t\tAMST adding tree edges = %.3f seconds\n",wtime);
3938
3939 /*
3940 add the heavy edges between every two subgraphs, if such an edge
3941 exists, and if the subraphs were not already connected through the
3942 tree
3943 */
3944
3945 wtime = taucs_wtime();
3946
3947 if (curr_group>1) /* more than 1 group */
3948 Bent = add_heavy_edges(mtxA,precond,Bent,array,groups,curr_group,pi,diag);
3949
3950 if (Bent == -1) { /* memory allocation failure in add_heavy_edges */
3951 taucs_free(pi);
3952 taucs_free(diag);
3953 taucs_free(groups);
3954 free_linked_list(lp);
3955 /* taucs_free(point_to_heap); */
3956 free_graph(mtxA);
3957 free_graph(precond);
3958 return NULL;
3959 }
3960
3961 wtime = taucs_wtime() - wtime;
3962 taucs_printf("\t\tAMST finding heavy edges = %.3f seconds\n",wtime);
3963
3964 /* allocate more memory to the preconditioner if needed */
3965
3966 wtime = taucs_wtime();
3967
3968 if (precond_maxnent < Bent + n) {
3969 taucs_printf("adding space for %d entries in precond for diagonal entries\n",Bent+n-precond_maxnent);
3970 precond_maxnent = Bent + n;
3971
3972 graph_resize(precond,precond_maxnent);
3973 precond->nent = Bent;
3974 /*
3975 ivec1_p = precond->ivec1 ;
3976 ivec2_p = precond->ivec2 ;
3977 dvec_p = precond->dvec ;
3978 */
3979 }
3980
3981 Do(i,n) {
3982
3983 (precond->edges)[Bent].i = i;
3984 (precond->edges)[Bent].j = i;
3985 (precond->edges)[Bent].v = diag[i];
3986
3987 Bent++;
3988 }
3989
3990 precond->nent = Bent;
3991
3992 wtime = taucs_wtime() - wtime;
3993 taucs_printf("\t\tAMST resize and add heavy edges = %.3f seconds\n",wtime);
3994
3995 wtime = taucs_wtime();
3996
3997 taucs_free(pi);
3998 taucs_free(diag);
3999 taucs_free(groups);
4000 free_linked_list(lp);
4001 taucs_printf("actual number of entries in preconditioner = %d\n",Bent);
4002 free_graph(mtxA);
4003
4004 /* out could be NULL, but we return out anyway */
4005
4006 out = graph_to_ccs_matrix(precond);
4007
4008 assert(out);
4009
4010 free_graph(precond);
4011
4012 wtime = taucs_wtime() - wtime;
4013 taucs_printf("\t\tAMST free memory and convert to ccs = %.3f seconds\n",wtime);
4014
4015 #define TAUCS_VAIDYA_RELAX_NONONO
4016 #ifdef TAUCS_VAIDYA_RELAX
4017 {
4018 taucs_ccs_matrix* A = ccs_mtxA;
4019 taucs_ccs_matrix* M = out;
4020 int j;
4021
4022 for (j=0; j<A->n; j++) {
4023 double dA, dM, modification;
4024 int i,ip;
4025 for (ip=(A->colptr)[j]; ip<(A->colptr)[j+1]; ip++) {
4026 i = (A->rowind)[ ip ];
4027 if (i==j) {
4028 dA = (A->values.d/*taucs_values*/)[ ip ];
4029 break;
4030 }
4031 }
4032
4033 for (ip=(M->colptr)[j]; ip<(M->colptr)[j+1]; ip++) {
4034 i = (M->rowind)[ ip ];
4035 if (i==j) {
4036 dM = (M->values.d/*taucs_values*/)[ ip ];
4037 break;
4038 }
4039 }
4040
4041 assert(dA >= dM);
4042
4043 modification = dA - dM;
4044 if (j < 30) printf(">>> %.4e %.4e\n",dA,dM);
4045
4046 (M->values.d/*taucs_values*/)[ ip ] += 0.01 * modification;
4047 }
4048 }
4049 #endif
4050
4051 return out;
4052 }
4053
4054
4055 taucs_ccs_matrix*
taucs_amwb_preconditioner_create(taucs_ccs_matrix * A,int rnd,double subgraphs,int stretch_flag)4056 taucs_amwb_preconditioner_create(taucs_ccs_matrix *A,
4057 int rnd,
4058 double subgraphs,
4059 int stretch_flag)
4060 {
4061 double wtime;
4062 double* diag;
4063 graph* G_A;
4064 int diagnostics;
4065 int n;
4066
4067 if (!(A->flags & TAUCS_DOUBLE)) {
4068 taucs_printf("taucs_amwb_preconditioner_create: matrix must be double-precision real\n");
4069 return NULL;
4070 }
4071
4072 n = A->n;
4073
4074 diag = (double*) taucs_malloc(n*sizeof(double));
4075 if (diag == NULL) return NULL;
4076
4077 wtime = taucs_wtime();
4078 G_A = ccs_matrix_to_graph_plus(A,&diagnostics,diag,1 /* force diag dominance */);
4079 if (!G_A) {
4080 taucs_free(diag);
4081 return NULL;
4082 }
4083 wtime = taucs_wtime() - wtime;
4084 taucs_printf("\t\tAMWB matrix-to-graph + analysis = %.3f seconds\n",wtime);
4085
4086 if (diagnostics & TAUCS_SYM_NOT_SYMLOWER) {
4087 taucs_printf("taucs_amwb_preconditioner_create: matrix must be symmetrix & lower\n");
4088 /* in this case, G_A == NULL, no need to free */
4089 taucs_free(diag);
4090 return A;
4091 }
4092 if (diagnostics & TAUCS_SYM_NOT_DIAGDOMINANT) {
4093 taucs_printf("taucs_amwb_preconditioner_create: matrix not diagonally dominant\n");
4094 taucs_free(diag);
4095 free_graph(G_A);
4096 return A;
4097 }
4098 if (diagnostics & TAUCS_SYM_NEG_DIAGONALS) {
4099 taucs_printf("taucs_amwb_preconditioner_create: negative diagonal elements\n");
4100 taucs_free(diag);
4101 free_graph(G_A);
4102 return A;
4103 }
4104
4105 if (diagnostics & TAUCS_SYM_POS_OFFDIAGONALS)
4106 return amwb_preconditioner_create(G_A, diag, rnd, subgraphs);
4107 else
4108 return amst_preconditioner_create(G_A, diag, rnd, subgraphs,stretch_flag);
4109 }
4110
4111 #endif /* TAUCS_CORE_DOUBLE */
4112
4113 /*********************************************************/
4114 /* */
4115 /*********************************************************/
4116
4117