1/* $Id$Revision: */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include <stdio.h>
15#include <string.h>
16#include <math.h>
17#include <assert.h>
18#include "logic.h"
19#include "memory.h"
20#include "arith.h"
21#include "SparseMatrix.h"
22#include "BinaryHeap.h"
23#if PQ
24#include "LinkedList.h"
25#include "PriorityQueue.h"
26#endif
27
28static size_t size_of_matrix_type(int type){
29 int size = 0;
30 switch (type){
31 case MATRIX_TYPE_REAL:
32 size = sizeof(real);
33 break;
34 case MATRIX_TYPE_COMPLEX:
35 size = 2*sizeof(real);
36 break;
37 case MATRIX_TYPE_INTEGER:
38 size = sizeof(int);
39 break;
40 case MATRIX_TYPE_PATTERN:
41 size = 0;
42 break;
43 case MATRIX_TYPE_UNKNOWN:
44 size = 0;
45 break;
46 default:
47 size = 0;
48 break;
49 }
50
51 return size;
52}
53
54SparseMatrix SparseMatrix_sort(SparseMatrix A){
55 SparseMatrix B;
56 B = SparseMatrix_transpose(A);
57 SparseMatrix_delete(A);
58 A = SparseMatrix_transpose(B);
59 SparseMatrix_delete(B);
60 return A;
61}
62SparseMatrix SparseMatrix_make_undirected(SparseMatrix A){
63 /* make it strictly low diag only, and set flag to undirected */
64 SparseMatrix B;
65 B = SparseMatrix_symmetrize(A, FALSE);
66 SparseMatrix_set_undirected(B);
67 return SparseMatrix_remove_upper(B);
68}
69SparseMatrix SparseMatrix_transpose(SparseMatrix A){
70 if (!A) return NULL;
71
72 int *ia = A->ia, *ja = A->ja, *ib, *jb, nz = A->nz, m = A->m, n = A->n, type = A->type, format = A->format;
73 SparseMatrix B;
74 int i, j;
75
76 assert(A->format == FORMAT_CSR);/* only implemented for CSR right now */
77
78 B = SparseMatrix_new(n, m, nz, type, format);
79 B->nz = nz;
80 ib = B->ia;
81 jb = B->ja;
82
83 for (i = 0; i <= n; i++) ib[i] = 0;
84 for (i = 0; i < m; i++){
85 for (j = ia[i]; j < ia[i+1]; j++){
86 ib[ja[j]+1]++;
87 }
88 }
89
90 for (i = 0; i < n; i++) ib[i+1] += ib[i];
91
92 switch (A->type){
93 case MATRIX_TYPE_REAL:{
94 real *a = (real*) A->a;
95 real *b = (real*) B->a;
96 for (i = 0; i < m; i++){
97 for (j = ia[i]; j < ia[i+1]; j++){
98 jb[ib[ja[j]]] = i;
99 b[ib[ja[j]]++] = a[j];
100 }
101 }
102 break;
103 }
104 case MATRIX_TYPE_COMPLEX:{
105 real *a = (real*) A->a;
106 real *b = (real*) B->a;
107 for (i = 0; i < m; i++){
108 for (j = ia[i]; j < ia[i+1]; j++){
109 jb[ib[ja[j]]] = i;
110 b[2*ib[ja[j]]] = a[2*j];
111 b[2*ib[ja[j]]+1] = a[2*j+1];
112 ib[ja[j]]++;
113 }
114 }
115 break;
116 }
117 case MATRIX_TYPE_INTEGER:{
118 int *ai = (int*) A->a;
119 int *bi = (int*) B->a;
120 for (i = 0; i < m; i++){
121 for (j = ia[i]; j < ia[i+1]; j++){
122 jb[ib[ja[j]]] = i;
123 bi[ib[ja[j]]++] = ai[j];
124 }
125 }
126 break;
127 }
128 case MATRIX_TYPE_PATTERN:
129 for (i = 0; i < m; i++){
130 for (j = ia[i]; j < ia[i+1]; j++){
131 jb[ib[ja[j]]++] = i;
132 }
133 }
134 break;
135 case MATRIX_TYPE_UNKNOWN:
136 SparseMatrix_delete(B);
137 return NULL;
138 default:
139 SparseMatrix_delete(B);
140 return NULL;
141 }
142
143
144 for (i = n-1; i >= 0; i--) ib[i+1] = ib[i];
145 ib[0] = 0;
146
147
148 return B;
149}
150
151SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only){
152 SparseMatrix B;
153 if (SparseMatrix_is_symmetric(A, pattern_symmetric_only)) return SparseMatrix_copy(A);
154 B = SparseMatrix_transpose(A);
155 if (!B) return NULL;
156 A = SparseMatrix_add(A, B);
157 SparseMatrix_delete(B);
158 SparseMatrix_set_symmetric(A);
159 SparseMatrix_set_pattern_symmetric(A);
160 return A;
161}
162
163SparseMatrix SparseMatrix_symmetrize_nodiag(SparseMatrix A, int pattern_symmetric_only){
164 SparseMatrix B;
165 if (SparseMatrix_is_symmetric(A, pattern_symmetric_only)) {
166 B = SparseMatrix_copy(A);
167 return SparseMatrix_remove_diagonal(B);
168 }
169 B = SparseMatrix_transpose(A);
170 if (!B) return NULL;
171 A = SparseMatrix_add(A, B);
172 SparseMatrix_delete(B);
173 SparseMatrix_set_symmetric(A);
174 SparseMatrix_set_pattern_symmetric(A);
175 return SparseMatrix_remove_diagonal(A);
176}
177
178int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only){
179 if (!A) return FALSE;
180
181 /* assume no repeated entries! */
182 SparseMatrix B;
183 int *ia, *ja, *ib, *jb, type, m;
184 int *mask;
185 int res = FALSE;
186 int i, j;
187 assert(A->format == FORMAT_CSR);/* only implemented for CSR right now */
188
189 if (SparseMatrix_known_symmetric(A)) return TRUE;
190 if (test_pattern_symmetry_only && SparseMatrix_known_strucural_symmetric(A)) return TRUE;
191
192 if (A->m != A->n) return FALSE;
193
194 B = SparseMatrix_transpose(A);
195 if (!B) return FALSE;
196
197 ia = A->ia;
198 ja = A->ja;
199 ib = B->ia;
200 jb = B->ja;
201 m = A->m;
202
203 mask = MALLOC(sizeof(int)*((size_t) m));
204 for (i = 0; i < m; i++) mask[i] = -1;
205
206 type = A->type;
207 if (test_pattern_symmetry_only) type = MATRIX_TYPE_PATTERN;
208
209 switch (type){
210 case MATRIX_TYPE_REAL:{
211 real *a = (real*) A->a;
212 real *b = (real*) B->a;
213 for (i = 0; i <= m; i++) if (ia[i] != ib[i]) goto RETURN;
214 for (i = 0; i < m; i++){
215 for (j = ia[i]; j < ia[i+1]; j++){
216 mask[ja[j]] = j;
217 }
218 for (j = ib[i]; j < ib[i+1]; j++){
219 if (mask[jb[j]] < ia[i]) goto RETURN;
220 }
221 for (j = ib[i]; j < ib[i+1]; j++){
222 if (ABS(b[j] - a[mask[jb[j]]]) > SYMMETRY_EPSILON) goto RETURN;
223 }
224 }
225 res = TRUE;
226 break;
227 }
228 case MATRIX_TYPE_COMPLEX:{
229 real *a = (real*) A->a;
230 real *b = (real*) B->a;
231 for (i = 0; i <= m; i++) if (ia[i] != ib[i]) goto RETURN;
232 for (i = 0; i < m; i++){
233 for (j = ia[i]; j < ia[i+1]; j++){
234 mask[ja[j]] = j;
235 }
236 for (j = ib[i]; j < ib[i+1]; j++){
237 if (mask[jb[j]] < ia[i]) goto RETURN;
238 }
239 for (j = ib[i]; j < ib[i+1]; j++){
240 if (ABS(b[2*j] - a[2*mask[jb[j]]]) > SYMMETRY_EPSILON) goto RETURN;
241 if (ABS(b[2*j+1] - a[2*mask[jb[j]]+1]) > SYMMETRY_EPSILON) goto RETURN;
242 }
243 }
244 res = TRUE;
245 break;
246 }
247 case MATRIX_TYPE_INTEGER:{
248 int *ai = (int*) A->a;
249 int *bi = (int*) B->a;
250 for (i = 0; i < m; i++){
251 for (j = ia[i]; j < ia[i+1]; j++){
252 mask[ja[j]] = j;
253 }
254 for (j = ib[i]; j < ib[i+1]; j++){
255 if (mask[jb[j]] < ia[i]) goto RETURN;
256 }
257 for (j = ib[i]; j < ib[i+1]; j++){
258 if (bi[j] != ai[mask[jb[j]]]) goto RETURN;
259 }
260 }
261 res = TRUE;
262 break;
263 }
264 case MATRIX_TYPE_PATTERN:
265 for (i = 0; i < m; i++){
266 for (j = ia[i]; j < ia[i+1]; j++){
267 mask[ja[j]] = j;
268 }
269 for (j = ib[i]; j < ib[i+1]; j++){
270 if (mask[jb[j]] < ia[i]) goto RETURN;
271 }
272 }
273 res = TRUE;
274 break;
275 case MATRIX_TYPE_UNKNOWN:
276 goto RETURN;
277 break;
278 default:
279 goto RETURN;
280 break;
281 }
282
283 if (test_pattern_symmetry_only){
284 SparseMatrix_set_pattern_symmetric(A);
285 } else {
286 SparseMatrix_set_symmetric(A);
287 SparseMatrix_set_pattern_symmetric(A);
288 }
289 RETURN:
290 FREE(mask);
291
292 SparseMatrix_delete(B);
293 return res;
294}
295
296static SparseMatrix SparseMatrix_init(int m, int n, int type, size_t sz, int format){
297 SparseMatrix A;
298
299
300 A = MALLOC(sizeof(struct SparseMatrix_struct));
301 A->m = m;
302 A->n = n;
303 A->nz = 0;
304 A->nzmax = 0;
305 A->type = type;
306 A->size = sz;
307 switch (format){
308 case FORMAT_COORD:
309 A->ia = NULL;
310 break;
311 case FORMAT_CSC:
312 case FORMAT_CSR:
313 default:
314 A->ia = MALLOC(sizeof(int)*((size_t)(m+1)));
315 }
316 A->ja = NULL;
317 A->a = NULL;
318 A->format = format;
319 A->property = 0;
320 clear_flag(A->property, MATRIX_PATTERN_SYMMETRIC);
321 clear_flag(A->property, MATRIX_SYMMETRIC);
322 clear_flag(A->property, MATRIX_SKEW);
323 clear_flag(A->property, MATRIX_HERMITIAN);
324 return A;
325}
326
327static SparseMatrix SparseMatrix_alloc(SparseMatrix A, int nz){
328 int format = A->format;
329 size_t nz_t = (size_t) nz; /* size_t is 64 bit on 64 bit machine. Using nz*A->size can overflow. */
330
331 A->a = NULL;
332 switch (format){
333 case FORMAT_COORD:
334 A->ia = MALLOC(sizeof(int)*nz_t);
335 A->ja = MALLOC(sizeof(int)*nz_t);
336 A->a = MALLOC(A->size*nz_t);
337 break;
338 case FORMAT_CSR:
339 case FORMAT_CSC:
340 default:
341 A->ja = MALLOC(sizeof(int)*nz_t);
342 if (A->size > 0 && nz_t > 0) {
343 A->a = MALLOC(A->size*nz_t);
344 }
345 break;
346 }
347 A->nzmax = nz;
348 return A;
349}
350
351static SparseMatrix SparseMatrix_realloc(SparseMatrix A, int nz){
352 int format = A->format;
353 size_t nz_t = (size_t) nz; /* size_t is 64 bit on 64 bit machine. Using nz*A->size can overflow. */
354
355 switch (format){
356 case FORMAT_COORD:
357 A->ia = REALLOC(A->ia, sizeof(int)*nz_t);
358 A->ja = REALLOC(A->ja, sizeof(int)*nz_t);
359 if (A->size > 0) {
360 if (A->a){
361 A->a = REALLOC(A->a, A->size*nz_t);
362 } else {
363 A->a = MALLOC(A->size*nz_t);
364 }
365 }
366 break;
367 case FORMAT_CSR:
368 case FORMAT_CSC:
369 default:
370 A->ja = REALLOC(A->ja, sizeof(int)*nz_t);
371 if (A->size > 0) {
372 if (A->a){
373 A->a = REALLOC(A->a, A->size*nz_t);
374 } else {
375 A->a = MALLOC(A->size*nz_t);
376 }
377 }
378 break;
379 }
380 A->nzmax = nz;
381 return A;
382}
383
384SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format){
385 /* return a sparse matrix skeleton with row dimension m and storage nz. If nz == 0,
386 only row pointers are allocated */
387 SparseMatrix A;
388 size_t sz;
389
390 sz = size_of_matrix_type(type);
391 A = SparseMatrix_init(m, n, type, sz, format);
392
393 if (nz > 0) A = SparseMatrix_alloc(A, nz);
394 return A;
395
396}
397SparseMatrix SparseMatrix_general_new(int m, int n, int nz, int type, size_t sz, int format){
398 /* return a sparse matrix skeleton with row dimension m and storage nz. If nz == 0,
399 only row pointers are allocated. this is more general and allow elements to be
400 any data structure, not just real/int/complex etc
401 */
402 SparseMatrix A;
403
404 A = SparseMatrix_init(m, n, type, sz, format);
405
406 if (nz > 0) A = SparseMatrix_alloc(A, nz);
407 return A;
408
409}
410
411void SparseMatrix_delete(SparseMatrix A){
412 /* return a sparse matrix skeleton with row dimension m and storage nz. If nz == 0,
413 only row pointers are allocated */
414 if (!A) return;
415 if (A->ia) FREE(A->ia);
416 if (A->ja) FREE(A->ja);
417 if (A->a) FREE(A->a);
418 FREE(A);
419}
420void SparseMatrix_print_csr(char *c, SparseMatrix A){
421 int *ia, *ja;
422 real *a;
423 int *ai;
424 int i, j, m = A->m;
425
426 assert (A->format == FORMAT_CSR);
427 printf("%s\n SparseArray[{",c);
428 ia = A->ia;
429 ja = A->ja;
430 a = A->a;
431 switch (A->type){
432 case MATRIX_TYPE_REAL:
433 a = (real*) A->a;
434 for (i = 0; i < m; i++){
435 for (j = ia[i]; j < ia[i+1]; j++){
436 printf("{%d, %d}->%f",i+1, ja[j]+1, a[j]);
437 if (j != ia[m]-1) printf(",");
438 }
439 }
440 break;
441 case MATRIX_TYPE_COMPLEX:
442 a = (real*) A->a;
443 for (i = 0; i < m; i++){
444 for (j = ia[i]; j < ia[i+1]; j++){
445 printf("{%d, %d}->%f + %f I",i+1, ja[j]+1, a[2*j], a[2*j+1]);
446 if (j != ia[m]-1) printf(",");
447 }
448 }
449 printf("\n");
450 break;
451 case MATRIX_TYPE_INTEGER:
452 ai = (int*) A->a;
453 for (i = 0; i < m; i++){
454 for (j = ia[i]; j < ia[i+1]; j++){
455 printf("{%d, %d}->%d",i+1, ja[j]+1, ai[j]);
456 if (j != ia[m]-1) printf(",");
457 }
458 }
459 printf("\n");
460 break;
461 case MATRIX_TYPE_PATTERN:
462 for (i = 0; i < m; i++){
463 for (j = ia[i]; j < ia[i+1]; j++){
464 printf("{%d, %d}->_",i+1, ja[j]+1);
465 if (j != ia[m]-1) printf(",");
466 }
467 }
468 printf("\n");
469 break;
470 case MATRIX_TYPE_UNKNOWN:
471 return;
472 default:
473 return;
474 }
475 printf("},{%d, %d}]\n", m, A->n);
476
477}
478
479
480
481void SparseMatrix_print_coord(char *c, SparseMatrix A){
482 int *ia, *ja;
483 real *a;
484 int *ai;
485 int i, m = A->m;
486
487 assert (A->format == FORMAT_COORD);
488 printf("%s\n SparseArray[{",c);
489 ia = A->ia;
490 ja = A->ja;
491 a = A->a;
492 switch (A->type){
493 case MATRIX_TYPE_REAL:
494 a = (real*) A->a;
495 for (i = 0; i < A->nz; i++){
496 printf("{%d, %d}->%f",ia[i]+1, ja[i]+1, a[i]);
497 if (i != A->nz - 1) printf(",");
498 }
499 printf("\n");
500 break;
501 case MATRIX_TYPE_COMPLEX:
502 a = (real*) A->a;
503 for (i = 0; i < A->nz; i++){
504 printf("{%d, %d}->%f + %f I",ia[i]+1, ja[i]+1, a[2*i], a[2*i+1]);
505 if (i != A->nz - 1) printf(",");
506 }
507 printf("\n");
508 break;
509 case MATRIX_TYPE_INTEGER:
510 ai = (int*) A->a;
511 for (i = 0; i < A->nz; i++){
512 printf("{%d, %d}->%d",ia[i]+1, ja[i]+1, ai[i]);
513 if (i != A->nz) printf(",");
514 }
515 printf("\n");
516 break;
517 case MATRIX_TYPE_PATTERN:
518 for (i = 0; i < A->nz; i++){
519 printf("{%d, %d}->_",ia[i]+1, ja[i]+1);
520 if (i != A->nz - 1) printf(",");
521 }
522 printf("\n");
523 break;
524 case MATRIX_TYPE_UNKNOWN:
525 return;
526 default:
527 return;
528 }
529 printf("},{%d, %d}]\n", m, A->n);
530
531}
532
533
534
535
536void SparseMatrix_print(char *c, SparseMatrix A){
537 switch (A->format){
538 case FORMAT_CSR:
539 SparseMatrix_print_csr(c, A);
540 break;
541 case FORMAT_CSC:
542 assert(0);/* not implemented yet...
543 SparseMatrix_print_csc(c, A);*/
544 break;
545 case FORMAT_COORD:
546 SparseMatrix_print_coord(c, A);
547 break;
548 default:
549 assert(0);
550 }
551}
552
553
554
555
556
557static void SparseMatrix_export_csr(FILE *f, SparseMatrix A){
558 int *ia, *ja;
559 real *a;
560 int *ai;
561 int i, j, m = A->m;
562
563 switch (A->type){
564 case MATRIX_TYPE_REAL:
565 fprintf(f,"%%%%MatrixMarket matrix coordinate real general\n");
566 break;
567 case MATRIX_TYPE_COMPLEX:
568 fprintf(f,"%%%%MatrixMarket matrix coordinate complex general\n");
569 break;
570 case MATRIX_TYPE_INTEGER:
571 fprintf(f,"%%%%MatrixMarket matrix coordinate integer general\n");
572 break;
573 case MATRIX_TYPE_PATTERN:
574 fprintf(f,"%%%%MatrixMarket matrix coordinate pattern general\n");
575 break;
576 case MATRIX_TYPE_UNKNOWN:
577 return;
578 default:
579 return;
580 }
581
582 fprintf(f,"%d %d %d\n",A->m,A->n,A->nz);
583 ia = A->ia;
584 ja = A->ja;
585 a = A->a;
586 switch (A->type){
587 case MATRIX_TYPE_REAL:
588 a = (real*) A->a;
589 for (i = 0; i < m; i++){
590 for (j = ia[i]; j < ia[i+1]; j++){
591 fprintf(f, "%d %d %16.8g\n",i+1, ja[j]+1, a[j]);
592 }
593 }
594 break;
595 case MATRIX_TYPE_COMPLEX:
596 a = (real*) A->a;
597 for (i = 0; i < m; i++){
598 for (j = ia[i]; j < ia[i+1]; j++){
599 fprintf(f, "%d %d %16.8g %16.8g\n",i+1, ja[j]+1, a[2*j], a[2*j+1]);
600 }
601 }
602 break;
603 case MATRIX_TYPE_INTEGER:
604 ai = (int*) A->a;
605 for (i = 0; i < m; i++){
606 for (j = ia[i]; j < ia[i+1]; j++){
607 fprintf(f, "%d %d %d\n",i+1, ja[j]+1, ai[j]);
608 }
609 }
610 break;
611 case MATRIX_TYPE_PATTERN:
612 for (i = 0; i < m; i++){
613 for (j = ia[i]; j < ia[i+1]; j++){
614 fprintf(f, "%d %d\n",i+1, ja[j]+1);
615 }
616 }
617 break;
618 case MATRIX_TYPE_UNKNOWN:
619 return;
620 default:
621 return;
622 }
623
624}
625
626void SparseMatrix_export_binary_fp(FILE *f, SparseMatrix A){
627
628 fwrite(&(A->m), sizeof(int), 1, f);
629 fwrite(&(A->n), sizeof(int), 1, f);
630 fwrite(&(A->nz), sizeof(int), 1, f);
631 fwrite(&(A->nzmax), sizeof(int), 1, f);
632 fwrite(&(A->type), sizeof(int), 1, f);
633 fwrite(&(A->format), sizeof(int), 1, f);
634 fwrite(&(A->property), sizeof(int), 1, f);
635 fwrite(&(A->size), sizeof(size_t), 1, f);
636 if (A->format == FORMAT_COORD){
637 fwrite(A->ia, sizeof(int), A->nz, f);
638 } else {
639 fwrite(A->ia, sizeof(int), A->m + 1, f);
640 }
641 fwrite(A->ja, sizeof(int), A->nz, f);
642 if (A->size > 0) fwrite(A->a, A->size, A->nz, f);
643
644}
645
646void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag){
647 FILE *f;
648
649 *flag = 0;
650 f = fopen(name, "wb");
651 if (!f) {
652 *flag = 1;
653 return;
654 }
655 SparseMatrix_export_binary_fp(f, A);
656 fclose(f);
657
658}
659
660
661
662SparseMatrix SparseMatrix_import_binary_fp(FILE *f){
663 SparseMatrix A = NULL;
664 int m, n, nz, nzmax, type, format, property, iread;
665 size_t sz;
666
667 iread = fread(&m, sizeof(int), 1, f);
668 if (iread != 1) return NULL;
669 iread = fread(&n, sizeof(int), 1, f);
670 if (iread != 1) return NULL;
671 iread = fread(&nz, sizeof(int), 1, f);
672 if (iread != 1) return NULL;
673 iread = fread(&nzmax, sizeof(int), 1, f);
674 if (iread != 1) return NULL;
675 iread = fread(&type, sizeof(int), 1, f);
676 if (iread != 1) return NULL;
677 iread = fread(&format, sizeof(int), 1, f);
678 if (iread != 1) return NULL;
679 iread = fread(&property, sizeof(int), 1, f);
680 if (iread != 1) return NULL;
681 iread = fread(&sz, sizeof(size_t), 1, f);
682 if (iread != 1) return NULL;
683
684 A = SparseMatrix_general_new(m, n, nz, type, sz, format);
685 A->nz = nz;
686 A->property = property;
687
688 if (format == FORMAT_COORD){
689 iread = fread(A->ia, sizeof(int), A->nz, f);
690 if (iread != A->nz) return NULL;
691 } else {
692 iread = fread(A->ia, sizeof(int), A->m + 1, f);
693 if (iread != A->m + 1) return NULL;
694 }
695 iread = fread(A->ja, sizeof(int), A->nz, f);
696 if (iread != A->nz) return NULL;
697
698 if (A->size > 0) {
699 iread = fread(A->a, A->size, A->nz, f);
700 if (iread != A->nz) return NULL;
701 }
702 fclose(f);
703 return A;
704}
705
706
707SparseMatrix SparseMatrix_import_binary(char *name){
708 SparseMatrix A = NULL;
709 FILE *f;
710 f = fopen(name, "rb");
711
712 A = SparseMatrix_import_binary_fp(f);
713 return A;
714}
715
716static void SparseMatrix_export_coord(FILE *f, SparseMatrix A){
717 int *ia, *ja;
718 real *a;
719 int *ai;
720 int i;
721
722 switch (A->type){
723 case MATRIX_TYPE_REAL:
724 fprintf(f,"%%%%MatrixMarket matrix coordinate real general\n");
725 break;
726 case MATRIX_TYPE_COMPLEX:
727 fprintf(f,"%%%%MatrixMarket matrix coordinate complex general\n");
728 break;
729 case MATRIX_TYPE_INTEGER:
730 fprintf(f,"%%%%MatrixMarket matrix coordinate integer general\n");
731 break;
732 case MATRIX_TYPE_PATTERN:
733 fprintf(f,"%%%%MatrixMarket matrix coordinate pattern general\n");
734 break;
735 case MATRIX_TYPE_UNKNOWN:
736 return;
737 default:
738 return;
739 }
740
741 fprintf(f,"%d %d %d\n",A->m,A->n,A->nz);
742 ia = A->ia;
743 ja = A->ja;
744 a = A->a;
745 switch (A->type){
746 case MATRIX_TYPE_REAL:
747 a = (real*) A->a;
748 for (i = 0; i < A->nz; i++){
749 fprintf(f, "%d %d %16.8g\n",ia[i]+1, ja[i]+1, a[i]);
750 }
751 break;
752 case MATRIX_TYPE_COMPLEX:
753 a = (real*) A->a;
754 for (i = 0; i < A->nz; i++){
755 fprintf(f, "%d %d %16.8g %16.8g\n",ia[i]+1, ja[i]+1, a[2*i], a[2*i+1]);
756 }
757 break;
758 case MATRIX_TYPE_INTEGER:
759 ai = (int*) A->a;
760 for (i = 0; i < A->nz; i++){
761 fprintf(f, "%d %d %d\n",ia[i]+1, ja[i]+1, ai[i]);
762 }
763 break;
764 case MATRIX_TYPE_PATTERN:
765 for (i = 0; i < A->nz; i++){
766 fprintf(f, "%d %d\n",ia[i]+1, ja[i]+1);
767 }
768 break;
769 case MATRIX_TYPE_UNKNOWN:
770 return;
771 default:
772 return;
773 }
774}
775
776
777
778void SparseMatrix_export(FILE *f, SparseMatrix A){
779
780 switch (A->format){
781 case FORMAT_CSR:
782 SparseMatrix_export_csr(f, A);
783 break;
784 case FORMAT_CSC:
785 assert(0);/* not implemented yet...
786 SparseMatrix_print_csc(c, A);*/
787 break;
788 case FORMAT_COORD:
789 SparseMatrix_export_coord(f, A);
790 break;
791 default:
792 assert(0);
793 }
794}
795
796
797SparseMatrix SparseMatrix_from_coordinate_format(SparseMatrix A){
798 /* convert a sparse matrix in coordinate form to one in compressed row form.*/
799 int *irn, *jcn;
800
801 void *a = A->a;
802
803 assert(A->format == FORMAT_COORD);
804 if (A->format != FORMAT_COORD) {
805 return NULL;
806 }
807 irn = A->ia;
808 jcn = A->ja;
809 return SparseMatrix_from_coordinate_arrays(A->nz, A->m, A->n, irn, jcn, a, A->type, A->size);
810
811}
812SparseMatrix SparseMatrix_from_coordinate_format_not_compacted(SparseMatrix A, int what_to_sum){
813 /* convert a sparse matrix in coordinate form to one in compressed row form.*/
814 int *irn, *jcn;
815
816 void *a = A->a;
817
818 assert(A->format == FORMAT_COORD);
819 if (A->format != FORMAT_COORD) {
820 return NULL;
821 }
822 irn = A->ia;
823 jcn = A->ja;
824 return SparseMatrix_from_coordinate_arrays_not_compacted(A->nz, A->m, A->n, irn, jcn, a, A->type, A->size, what_to_sum);
825
826}
827
828static SparseMatrix SparseMatrix_from_coordinate_arrays_internal(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz, int sum_repeated){
829 /* convert a sparse matrix in coordinate form to one in compressed row form.
830 nz: number of entries
831 irn: row indices 0-based
832 jcn: column indices 0-based
833 val values if not NULL
834 type: matrix type
835 */
836
837 SparseMatrix A = NULL;
838 int *ia, *ja;
839 real *a, *val;
840 int *ai, *vali;
841 int i;
842
843 assert(m > 0 && n > 0 && nz >= 0);
844
845 if (m <=0 || n <= 0 || nz < 0) return NULL;
846 A = SparseMatrix_general_new(m, n, nz, type, sz, FORMAT_CSR);
847 assert(A);
848 if (!A) return NULL;
849 ia = A->ia;
850 ja = A->ja;
851
852 for (i = 0; i <= m; i++){
853 ia[i] = 0;
854 }
855
856 switch (type){
857 case MATRIX_TYPE_REAL:
858 val = (real*) val0;
859 a = (real*) A->a;
860 for (i = 0; i < nz; i++){
861 if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
862 assert(0);
863 return NULL;
864 }
865 ia[irn[i]+1]++;
866 }
867 for (i = 0; i < m; i++) ia[i+1] += ia[i];
868 for (i = 0; i < nz; i++){
869 a[ia[irn[i]]] = val[i];
870 ja[ia[irn[i]]++] = jcn[i];
871 }
872 for (i = m; i > 0; i--) ia[i] = ia[i - 1];
873 ia[0] = 0;
874 break;
875 case MATRIX_TYPE_COMPLEX:
876 val = (real*) val0;
877 a = (real*) A->a;
878 for (i = 0; i < nz; i++){
879 if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
880 assert(0);
881 return NULL;
882 }
883 ia[irn[i]+1]++;
884 }
885 for (i = 0; i < m; i++) ia[i+1] += ia[i];
886 for (i = 0; i < nz; i++){
887 a[2*ia[irn[i]]] = *(val++);
888 a[2*ia[irn[i]]+1] = *(val++);
889 ja[ia[irn[i]]++] = jcn[i];
890 }
891 for (i = m; i > 0; i--) ia[i] = ia[i - 1];
892 ia[0] = 0;
893 break;
894 case MATRIX_TYPE_INTEGER:
895 vali = (int*) val0;
896 ai = (int*) A->a;
897 for (i = 0; i < nz; i++){
898 if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
899 assert(0);
900 return NULL;
901 }
902 ia[irn[i]+1]++;
903 }
904 for (i = 0; i < m; i++) ia[i+1] += ia[i];
905 for (i = 0; i < nz; i++){
906 ai[ia[irn[i]]] = vali[i];
907 ja[ia[irn[i]]++] = jcn[i];
908 }
909 for (i = m; i > 0; i--) ia[i] = ia[i - 1];
910 ia[0] = 0;
911 break;
912 case MATRIX_TYPE_PATTERN:
913 for (i = 0; i < nz; i++){
914 if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
915 assert(0);
916 return NULL;
917 }
918 ia[irn[i]+1]++;
919 }
920 for (i = 0; i < m; i++) ia[i+1] += ia[i];
921 for (i = 0; i < nz; i++){
922 ja[ia[irn[i]]++] = jcn[i];
923 }
924 for (i = m; i > 0; i--) ia[i] = ia[i - 1];
925 ia[0] = 0;
926 break;
927 case MATRIX_TYPE_UNKNOWN:
928 for (i = 0; i < nz; i++){
929 if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
930 assert(0);
931 return NULL;
932 }
933 ia[irn[i]+1]++;
934 }
935 for (i = 0; i < m; i++) ia[i+1] += ia[i];
936 MEMCPY(A->a, val0, A->size*((size_t)nz));
937 for (i = 0; i < nz; i++){
938 ja[ia[irn[i]]++] = jcn[i];
939 }
940 for (i = m; i > 0; i--) ia[i] = ia[i - 1];
941 ia[0] = 0;
942 break;
943 default:
944 assert(0);
945 return NULL;
946 }
947 A->nz = nz;
948
949
950
951 if(sum_repeated) A = SparseMatrix_sum_repeat_entries(A, sum_repeated);
952
953 return A;
954}
955
956
957SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz){
958 return SparseMatrix_from_coordinate_arrays_internal(nz, m, n, irn, jcn, val0, type, sz, SUM_REPEATED_ALL);
959}
960
961
962SparseMatrix SparseMatrix_from_coordinate_arrays_not_compacted(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz, int what_to_sum){
963 return SparseMatrix_from_coordinate_arrays_internal(nz, m, n, irn, jcn, val0, type, sz, what_to_sum);
964}
965
966SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B){
967 int m, n;
968 SparseMatrix C = NULL;
969 int *mask = NULL;
970 int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic, *jc;
971 int i, j, nz, nzmax;
972
973 assert(A && B);
974 assert(A->format == B->format && A->format == FORMAT_CSR);/* other format not yet supported */
975 assert(A->type == B->type);
976 m = A->m;
977 n = A->n;
978 if (m != B->m || n != B->n) return NULL;
979
980 nzmax = A->nz + B->nz;/* just assume that no entries overlaps for speed */
981
982 C = SparseMatrix_new(m, n, nzmax, A->type, FORMAT_CSR);
983 if (!C) goto RETURN;
984 ic = C->ia;
985 jc = C->ja;
986
987 mask = MALLOC(sizeof(int)*((size_t) n));
988
989 for (i = 0; i < n; i++) mask[i] = -1;
990
991 nz = 0;
992 ic[0] = 0;
993 switch (A->type){
994 case MATRIX_TYPE_REAL:{
995 real *a = (real*) A->a;
996 real *b = (real*) B->a;
997 real *c = (real*) C->a;
998 for (i = 0; i < m; i++){
999 for (j = ia[i]; j < ia[i+1]; j++){
1000 mask[ja[j]] = nz;
1001 jc[nz] = ja[j];
1002 c[nz] = a[j];
1003 nz++;
1004 }
1005 for (j = ib[i]; j < ib[i+1]; j++){
1006 if (mask[jb[j]] < ic[i]){
1007 jc[nz] = jb[j];
1008 c[nz++] = b[j];
1009 } else {
1010 c[mask[jb[j]]] += b[j];
1011 }
1012 }
1013 ic[i+1] = nz;
1014 }
1015 break;
1016 }
1017 case MATRIX_TYPE_COMPLEX:{
1018 real *a = (real*) A->a;
1019 real *b = (real*) B->a;
1020 real *c = (real*) C->a;
1021 for (i = 0; i < m; i++){
1022 for (j = ia[i]; j < ia[i+1]; j++){
1023 mask[ja[j]] = nz;
1024 jc[nz] = ja[j];
1025 c[2*nz] = a[2*j];
1026 c[2*nz+1] = a[2*j+1];
1027 nz++;
1028 }
1029 for (j = ib[i]; j < ib[i+1]; j++){
1030 if (mask[jb[j]] < ic[i]){
1031 jc[nz] = jb[j];
1032 c[2*nz] = b[2*j];
1033 c[2*nz+1] = b[2*j+1];
1034 nz++;
1035 } else {
1036 c[2*mask[jb[j]]] += b[2*j];
1037 c[2*mask[jb[j]]+1] += b[2*j+1];
1038 }
1039 }
1040 ic[i+1] = nz;
1041 }
1042 break;
1043 }
1044 case MATRIX_TYPE_INTEGER:{
1045 int *a = (int*) A->a;
1046 int *b = (int*) B->a;
1047 int *c = (int*) C->a;
1048 for (i = 0; i < m; i++){
1049 for (j = ia[i]; j < ia[i+1]; j++){
1050 mask[ja[j]] = nz;
1051 jc[nz] = ja[j];
1052 c[nz] = a[j];
1053 nz++;
1054 }
1055 for (j = ib[i]; j < ib[i+1]; j++){
1056 if (mask[jb[j]] < ic[i]){
1057 jc[nz] = jb[j];
1058 c[nz] = b[j];
1059 nz++;
1060 } else {
1061 c[mask[jb[j]]] += b[j];
1062 }
1063 }
1064 ic[i+1] = nz;
1065 }
1066 break;
1067 }
1068 case MATRIX_TYPE_PATTERN:{
1069 for (i = 0; i < m; i++){
1070 for (j = ia[i]; j < ia[i+1]; j++){
1071 mask[ja[j]] = nz;
1072 jc[nz] = ja[j];
1073 nz++;
1074 }
1075 for (j = ib[i]; j < ib[i+1]; j++){
1076 if (mask[jb[j]] < ic[i]){
1077 jc[nz] = jb[j];
1078 nz++;
1079 }
1080 }
1081 ic[i+1] = nz;
1082 }
1083 break;
1084 }
1085 case MATRIX_TYPE_UNKNOWN:
1086 break;
1087 default:
1088 break;
1089 }
1090 C->nz = nz;
1091
1092 RETURN:
1093 if (mask) FREE(mask);
1094
1095 return C;
1096}
1097
1098
1099
1100static void dense_transpose(real *v, int m, int n){
1101 /* transpose an m X n matrix in place. Well, we do no really do it without xtra memory. This is possibe, but too complicated for ow */
1102 int i, j;
1103 real *u;
1104 u = MALLOC(sizeof(real)*((size_t) m)*((size_t) n));
1105 MEMCPY(u,v, sizeof(real)*((size_t) m)*((size_t) n));
1106
1107 for (i = 0; i < m; i++){
1108 for (j = 0; j < n; j++){
1109 v[j*m+i] = u[i*n+j];
1110 }
1111 }
1112 FREE(u);
1113}
1114
1115
1116static void SparseMatrix_multiply_dense1(SparseMatrix A, real *v, real **res, int dim, int transposed, int res_transposed){
1117 /* A v or A^T v where v a dense matrix of second dimension dim. Real only for now. */
1118 int i, j, k, *ia, *ja, n, m;
1119 real *a, *u;
1120
1121 assert(A->format == FORMAT_CSR);
1122 assert(A->type == MATRIX_TYPE_REAL);
1123
1124 a = (real*) A->a;
1125 ia = A->ia;
1126 ja = A->ja;
1127 m = A->m;
1128 n = A->n;
1129 u = *res;
1130
1131 if (!transposed){
1132 if (!u) u = MALLOC(sizeof(real)*((size_t) m)*((size_t) dim));
1133 for (i = 0; i < m; i++){
1134 for (k = 0; k < dim; k++) u[i*dim+k] = 0.;
1135 for (j = ia[i]; j < ia[i+1]; j++){
1136 for (k = 0; k < dim; k++) u[i*dim+k] += a[j]*v[ja[j]*dim+k];
1137 }
1138 }
1139 if (res_transposed) dense_transpose(u, m, dim);
1140 } else {
1141 if (!u) u = MALLOC(sizeof(real)*((size_t) n)*((size_t) dim));
1142 for (i = 0; i < n*dim; i++) u[i] = 0.;
1143 for (i = 0; i < m; i++){
1144 for (j = ia[i]; j < ia[i+1]; j++){
1145 for (k = 0; k < dim; k++) u[ja[j]*dim + k] += a[j]*v[i*dim + k];
1146 }
1147 }
1148 if (res_transposed) dense_transpose(u, n, dim);
1149 }
1150
1151 *res = u;
1152
1153
1154}
1155
1156static void SparseMatrix_multiply_dense2(SparseMatrix A, real *v, real **res, int dim, int transposed, int res_transposed){
1157 /* A v^T or A^T v^T where v a dense matrix of second dimension n or m. Real only for now.
1158 transposed = FALSE: A*V^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension m x dim
1159 transposed = TRUE: A^T*V^T, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension n x dim
1160 */
1161 real *u, *rr;
1162 int i, m, n;
1163 assert(A->format == FORMAT_CSR);
1164 assert(A->type == MATRIX_TYPE_REAL);
1165 u = *res;
1166 m = A->m;
1167 n = A->n;
1168
1169 if (!transposed){
1170 if (!u) u = MALLOC(sizeof(real)*((size_t)m)*((size_t) dim));
1171 for (i = 0; i < dim; i++){
1172 rr = &(u[m*i]);
1173 SparseMatrix_multiply_vector(A, &(v[n*i]), &rr, transposed);
1174 }
1175 if (!res_transposed) dense_transpose(u, dim, m);
1176 } else {
1177 if (!u) u = MALLOC(sizeof(real)*((size_t)n)*((size_t)dim));
1178 for (i = 0; i < dim; i++){
1179 rr = &(u[n*i]);
1180 SparseMatrix_multiply_vector(A, &(v[m*i]), &rr, transposed);
1181 }
1182 if (!res_transposed) dense_transpose(u, dim, n);
1183 }
1184
1185 *res = u;
1186
1187
1188}
1189
1190
1191
1192void SparseMatrix_multiply_dense(SparseMatrix A, int ATransposed, real *v, int vTransposed, real **res, int res_transposed, int dim){
1193 /* depend on value of {ATranspose, vTransposed}, assume res_transposed == FALSE
1194 {FALSE, FALSE}: A * V, with A dimension m x n, with V of dimension n x dim. v[i*dim+j] gives V[i,j]. Result of dimension m x dim
1195 {TRUE, FALSE}: A^T * V, with A dimension m x n, with V of dimension m x dim. v[i*dim+j] gives V[i,j]. Result of dimension n x dim
1196 {FALSE, TRUE}: A*V^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension m x dim
1197 {TRUE, TRUE}: A^T*V^T, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension n x dim
1198
1199 furthermore, if res_transpose d== TRUE, then the result is transposed. Hence if res_transposed == TRUE
1200
1201 {FALSE, FALSE}: V^T A^T, with A dimension m x n, with V of dimension n x dim. v[i*dim+j] gives V[i,j]. Result of dimension dim x dim
1202 {TRUE, FALSE}: V^T A, with A dimension m x n, with V of dimension m x dim. v[i*dim+j] gives V[i,j]. Result of dimension dim x n
1203 {FALSE, TRUE}: V*A^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension dim x m
1204 {TRUE, TRUE}: V A, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension dim x n
1205 */
1206
1207 if (!vTransposed) {
1208 SparseMatrix_multiply_dense1(A, v, res, dim, ATransposed, res_transposed);
1209 } else {
1210 SparseMatrix_multiply_dense2(A, v, res, dim, ATransposed, res_transposed);
1211 }
1212
1213}
1214
1215
1216
1217void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int transposed){
1218 /* A v or A^T v. Real only for now. */
1219 int i, j, *ia, *ja, n, m;
1220 real *a, *u = NULL;
1221 int *ai;
1222 assert(A->format == FORMAT_CSR);
1223 assert(A->type == MATRIX_TYPE_REAL || A->type == MATRIX_TYPE_INTEGER);
1224
1225 ia = A->ia;
1226 ja = A->ja;
1227 m = A->m;
1228 n = A->n;
1229 u = *res;
1230
1231 switch (A->type){
1232 case MATRIX_TYPE_REAL:
1233 a = (real*) A->a;
1234 if (v){
1235 if (!transposed){
1236 if (!u) u = MALLOC(sizeof(real)*((size_t)m));
1237 for (i = 0; i < m; i++){
1238 u[i] = 0.;
1239 for (j = ia[i]; j < ia[i+1]; j++){
1240 u[i] += a[j]*v[ja[j]];
1241 }
1242 }
1243 } else {
1244 if (!u) u = MALLOC(sizeof(real)*((size_t)n));
1245 for (i = 0; i < n; i++) u[i] = 0.;
1246 for (i = 0; i < m; i++){
1247 for (j = ia[i]; j < ia[i+1]; j++){
1248 u[ja[j]] += a[j]*v[i];
1249 }
1250 }
1251 }
1252 } else {
1253 /* v is assumed to be all 1's */
1254 if (!transposed){
1255 if (!u) u = MALLOC(sizeof(real)*((size_t)m));
1256 for (i = 0; i < m; i++){
1257 u[i] = 0.;
1258 for (j = ia[i]; j < ia[i+1]; j++){
1259 u[i] += a[j];
1260 }
1261 }
1262 } else {
1263 if (!u) u = MALLOC(sizeof(real)*((size_t)n));
1264 for (i = 0; i < n; i++) u[i] = 0.;
1265 for (i = 0; i < m; i++){
1266 for (j = ia[i]; j < ia[i+1]; j++){
1267 u[ja[j]] += a[j];
1268 }
1269 }
1270 }
1271 }
1272 break;
1273 case MATRIX_TYPE_INTEGER:
1274 ai = (int*) A->a;
1275 if (v){
1276 if (!transposed){
1277 if (!u) u = MALLOC(sizeof(real)*((size_t)m));
1278 for (i = 0; i < m; i++){
1279 u[i] = 0.;
1280 for (j = ia[i]; j < ia[i+1]; j++){
1281 u[i] += ai[j]*v[ja[j]];
1282 }
1283 }
1284 } else {
1285 if (!u) u = MALLOC(sizeof(real)*((size_t)n));
1286 for (i = 0; i < n; i++) u[i] = 0.;
1287 for (i = 0; i < m; i++){
1288 for (j = ia[i]; j < ia[i+1]; j++){
1289 u[ja[j]] += ai[j]*v[i];
1290 }
1291 }
1292 }
1293 } else {
1294 /* v is assumed to be all 1's */
1295 if (!transposed){
1296 if (!u) u = MALLOC(sizeof(real)*((size_t)m));
1297 for (i = 0; i < m; i++){
1298 u[i] = 0.;
1299 for (j = ia[i]; j < ia[i+1]; j++){
1300 u[i] += ai[j];
1301 }
1302 }
1303 } else {
1304 if (!u) u = MALLOC(sizeof(real)*((size_t)n));
1305 for (i = 0; i < n; i++) u[i] = 0.;
1306 for (i = 0; i < m; i++){
1307 for (j = ia[i]; j < ia[i+1]; j++){
1308 u[ja[j]] += ai[j];
1309 }
1310 }
1311 }
1312 }
1313 break;
1314 default:
1315 assert(0);
1316 u = NULL;
1317 }
1318 *res = u;
1319
1320}
1321
1322
1323
1324SparseMatrix SparseMatrix_scaled_by_vector(SparseMatrix A, real *v, int apply_to_row){
1325 /* A SCALED BY VECOTR V IN ROW/COLUMN. Real only for now. */
1326 int i, j, *ia, *ja, m;
1327 real *a;
1328 assert(A->format == FORMAT_CSR);
1329 assert(A->type == MATRIX_TYPE_REAL);
1330
1331 a = (real*) A->a;
1332 ia = A->ia;
1333 ja = A->ja;
1334 m = A->m;
1335
1336
1337 if (!apply_to_row){
1338 for (i = 0; i < m; i++){
1339 for (j = ia[i]; j < ia[i+1]; j++){
1340 a[j] *= v[ja[j]];
1341 }
1342 }
1343 } else {
1344 for (i = 0; i < m; i++){
1345 if (v[i] != 0){
1346 for (j = ia[i]; j < ia[i+1]; j++){
1347 a[j] *= v[i];
1348 }
1349 }
1350 }
1351 }
1352 return A;
1353
1354}
1355SparseMatrix SparseMatrix_multiply_by_scaler(SparseMatrix A, real s){
1356 /* A scaled by a number */
1357 int i, j, *ia, m;
1358 real *a, *b = NULL;
1359 int *ai;
1360 assert(A->format == FORMAT_CSR);
1361
1362 switch (A->type){
1363 case MATRIX_TYPE_INTEGER:
1364 b = MALLOC(sizeof(real)*A->nz);
1365 ai = (int*) A->a;
1366 for (i = 0; i < A->nz; i++) b[i] = ai[i];
1367 FREE(A->a);
1368 A->a = b;
1369 A->type = MATRIX_TYPE_REAL;
1370 case MATRIX_TYPE_REAL:
1371 a = (real*) A->a;
1372 ia = A->ia;
1373 m = A->m;
1374 for (i = 0; i < m; i++){
1375 for (j = ia[i]; j < ia[i+1]; j++){
1376 a[j] *= s;
1377 }
1378 }
1379 break;
1380 case MATRIX_TYPE_COMPLEX:
1381 a = (real*) A->a;
1382 ia = A->ia;
1383 m = A->m;
1384 for (i = 0; i < m; i++){
1385 for (j = ia[i]; j < ia[i+1]; j++){
1386 a[2*j] *= s;
1387 a[2*j+1] *= s;
1388 }
1389 }
1390
1391 break;
1392 default:
1393 fprintf(stderr,"warning: scaling of matrix this type is not supported\n");
1394 }
1395
1396 return A;
1397
1398}
1399
1400
1401SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){
1402 int m;
1403 SparseMatrix C = NULL;
1404 int *mask = NULL;
1405 int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic, *jc;
1406 int i, j, k, jj, type, nz;
1407
1408 assert(A->format == B->format && A->format == FORMAT_CSR);/* other format not yet supported */
1409
1410 m = A->m;
1411 if (A->n != B->m) return NULL;
1412 if (A->type != B->type){
1413#ifdef DEBUG
1414 printf("in SparseMatrix_multiply, the matrix types do not match, right now only multiplication of matrices of the same type is supported\n");
1415#endif
1416 return NULL;
1417 }
1418 type = A->type;
1419
1420 mask = MALLOC(sizeof(int)*((size_t)(B->n)));
1421 if (!mask) return NULL;
1422
1423 for (i = 0; i < B->n; i++) mask[i] = -1;
1424
1425 nz = 0;
1426 for (i = 0; i < m; i++){
1427 for (j = ia[i]; j < ia[i+1]; j++){
1428 jj = ja[j];
1429 for (k = ib[jj]; k < ib[jj+1]; k++){
1430 if (mask[jb[k]] != -i - 2){
1431 if ((nz+1) <= nz) {
1432#ifdef DEBUG_PRINT
1433 fprintf(stderr,"overflow in SparseMatrix_multiply !!!\n");
1434#endif
1435 return NULL;
1436 }
1437 nz++;
1438 mask[jb[k]] = -i - 2;
1439 }
1440 }
1441 }
1442 }
1443
1444 C = SparseMatrix_new(m, B->n, nz, type, FORMAT_CSR);
1445 if (!C) goto RETURN;
1446 ic = C->ia;
1447 jc = C->ja;
1448
1449 nz = 0;
1450
1451 switch (type){
1452 case MATRIX_TYPE_REAL:
1453 {
1454 real *a = (real*) A->a;
1455 real *b = (real*) B->a;
1456 real *c = (real*) C->a;
1457 ic[0] = 0;
1458 for (i = 0; i < m; i++){
1459 for (j = ia[i]; j < ia[i+1]; j++){
1460 jj = ja[j];
1461 for (k = ib[jj]; k < ib[jj+1]; k++){
1462 if (mask[jb[k]] < ic[i]){
1463 mask[jb[k]] = nz;
1464 jc[nz] = jb[k];
1465 c[nz] = a[j]*b[k];
1466 nz++;
1467 } else {
1468 assert(jc[mask[jb[k]]] == jb[k]);
1469 c[mask[jb[k]]] += a[j]*b[k];
1470 }
1471 }
1472 }
1473 ic[i+1] = nz;
1474 }
1475 }
1476 break;
1477 case MATRIX_TYPE_COMPLEX:
1478 {
1479 real *a = (real*) A->a;
1480 real *b = (real*) B->a;
1481 real *c = (real*) C->a;
1482 a = (real*) A->a;
1483 b = (real*) B->a;
1484 c = (real*) C->a;
1485 ic[0] = 0;
1486 for (i = 0; i < m; i++){
1487 for (j = ia[i]; j < ia[i+1]; j++){
1488 jj = ja[j];
1489 for (k = ib[jj]; k < ib[jj+1]; k++){
1490 if (mask[jb[k]] < ic[i]){
1491 mask[jb[k]] = nz;
1492 jc[nz] = jb[k];
1493 c[2*nz] = a[2*j]*b[2*k] - a[2*j+1]*b[2*k+1];/*real part */
1494 c[2*nz+1] = a[2*j]*b[2*k+1] + a[2*j+1]*b[2*k];/*img part */
1495 nz++;
1496 } else {
1497 assert(jc[mask[jb[k]]] == jb[k]);
1498 c[2*mask[jb[k]]] += a[2*j]*b[2*k] - a[2*j+1]*b[2*k+1];/*real part */
1499 c[2*mask[jb[k]]+1] += a[2*j]*b[2*k+1] + a[2*j+1]*b[2*k];/*img part */
1500 }
1501 }
1502 }
1503 ic[i+1] = nz;
1504 }
1505 }
1506 break;
1507 case MATRIX_TYPE_INTEGER:
1508 {
1509 int *a = (int*) A->a;
1510 int *b = (int*) B->a;
1511 int *c = (int*) C->a;
1512 ic[0] = 0;
1513 for (i = 0; i < m; i++){
1514 for (j = ia[i]; j < ia[i+1]; j++){
1515 jj = ja[j];
1516 for (k = ib[jj]; k < ib[jj+1]; k++){
1517 if (mask[jb[k]] < ic[i]){
1518 mask[jb[k]] = nz;
1519 jc[nz] = jb[k];
1520 c[nz] = a[j]*b[k];
1521 nz++;
1522 } else {
1523 assert(jc[mask[jb[k]]] == jb[k]);
1524 c[mask[jb[k]]] += a[j]*b[k];
1525 }
1526 }
1527 }
1528 ic[i+1] = nz;
1529 }
1530 }
1531 break;
1532 case MATRIX_TYPE_PATTERN:
1533 ic[0] = 0;
1534 for (i = 0; i < m; i++){
1535 for (j = ia[i]; j < ia[i+1]; j++){
1536 jj = ja[j];
1537 for (k = ib[jj]; k < ib[jj+1]; k++){
1538 if (mask[jb[k]] < ic[i]){
1539 mask[jb[k]] = nz;
1540 jc[nz] = jb[k];
1541 nz++;
1542 } else {
1543 assert(jc[mask[jb[k]]] == jb[k]);
1544 }
1545 }
1546 }
1547 ic[i+1] = nz;
1548 }
1549 break;
1550 case MATRIX_TYPE_UNKNOWN:
1551 default:
1552 SparseMatrix_delete(C);
1553 C = NULL; goto RETURN;
1554 break;
1555 }
1556
1557 C->nz = nz;
1558
1559 RETURN:
1560 FREE(mask);
1561 return C;
1562
1563}
1564
1565
1566
1567SparseMatrix SparseMatrix_multiply3(SparseMatrix A, SparseMatrix B, SparseMatrix C){
1568 int m;
1569 SparseMatrix D = NULL;
1570 int *mask = NULL;
1571 int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic = C->ia, *jc = C->ja, *id, *jd;
1572 int i, j, k, l, ll, jj, type, nz;
1573
1574 assert(A->format == B->format && A->format == FORMAT_CSR);/* other format not yet supported */
1575
1576 m = A->m;
1577 if (A->n != B->m) return NULL;
1578 if (B->n != C->m) return NULL;
1579
1580 if (A->type != B->type || B->type != C->type){
1581#ifdef DEBUG
1582 printf("in SparseMatrix_multiply, the matrix types do not match, right now only multiplication of matrices of the same type is supported\n");
1583#endif
1584 return NULL;
1585 }
1586 type = A->type;
1587
1588 mask = MALLOC(sizeof(int)*((size_t)(C->n)));
1589 if (!mask) return NULL;
1590
1591 for (i = 0; i < C->n; i++) mask[i] = -1;
1592
1593 nz = 0;
1594 for (i = 0; i < m; i++){
1595 for (j = ia[i]; j < ia[i+1]; j++){
1596 jj = ja[j];
1597 for (l = ib[jj]; l < ib[jj+1]; l++){
1598 ll = jb[l];
1599 for (k = ic[ll]; k < ic[ll+1]; k++){
1600 if (mask[jc[k]] != -i - 2){
1601 if ((nz+1) <= nz) {
1602#ifdef DEBUG_PRINT
1603 fprintf(stderr,"overflow in SparseMatrix_multiply !!!\n");
1604#endif
1605 return NULL;
1606 }
1607 nz++;
1608 mask[jc[k]] = -i - 2;
1609 }
1610 }
1611 }
1612 }
1613 }
1614
1615 D = SparseMatrix_new(m, C->n, nz, type, FORMAT_CSR);
1616 if (!D) goto RETURN;
1617 id = D->ia;
1618 jd = D->ja;
1619
1620 nz = 0;
1621
1622 switch (type){
1623 case MATRIX_TYPE_REAL:
1624 {
1625 real *a = (real*) A->a;
1626 real *b = (real*) B->a;
1627 real *c = (real*) C->a;
1628 real *d = (real*) D->a;
1629 id[0] = 0;
1630 for (i = 0; i < m; i++){
1631 for (j = ia[i]; j < ia[i+1]; j++){
1632 jj = ja[j];
1633 for (l = ib[jj]; l < ib[jj+1]; l++){
1634 ll = jb[l];
1635 for (k = ic[ll]; k < ic[ll+1]; k++){
1636 if (mask[jc[k]] < id[i]){
1637 mask[jc[k]] = nz;
1638 jd[nz] = jc[k];
1639 d[nz] = a[j]*b[l]*c[k];
1640 nz++;
1641 } else {
1642 assert(jd[mask[jc[k]]] == jc[k]);
1643 d[mask[jc[k]]] += a[j]*b[l]*c[k];
1644 }
1645 }
1646 }
1647 }
1648 id[i+1] = nz;
1649 }
1650 }
1651 break;
1652 case MATRIX_TYPE_COMPLEX:
1653 {
1654 real *a = (real*) A->a;
1655 real *b = (real*) B->a;
1656 real *c = (real*) C->a;
1657 real *d = (real*) D->a;
1658 id[0] = 0;
1659 for (i = 0; i < m; i++){
1660 for (j = ia[i]; j < ia[i+1]; j++){
1661 jj = ja[j];
1662 for (l = ib[jj]; l < ib[jj+1]; l++){
1663 ll = jb[l];
1664 for (k = ic[ll]; k < ic[ll+1]; k++){
1665 if (mask[jc[k]] < id[i]){
1666 mask[jc[k]] = nz;
1667 jd[nz] = jc[k];
1668 d[2*nz] = (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k]
1669 - (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k+1];/*real part */
1670 d[2*nz+1] = (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k]
1671 + (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k+1];/*img part */
1672 nz++;
1673 } else {
1674 assert(jd[mask[jc[k]]] == jc[k]);
1675 d[2*mask[jc[k]]] += (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k]
1676 - (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k+1];/*real part */
1677 d[2*mask[jc[k]]+1] += (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k]
1678 + (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k+1];/*img part */
1679 }
1680 }
1681 }
1682 }
1683 id[i+1] = nz;
1684 }
1685 }
1686 break;
1687 case MATRIX_TYPE_INTEGER:
1688 {
1689 int *a = (int*) A->a;
1690 int *b = (int*) B->a;
1691 int *c = (int*) C->a;
1692 int *d = (int*) D->a;
1693 id[0] = 0;
1694 for (i = 0; i < m; i++){
1695 for (j = ia[i]; j < ia[i+1]; j++){
1696 jj = ja[j];
1697 for (l = ib[jj]; l < ib[jj+1]; l++){
1698 ll = jb[l];
1699 for (k = ic[ll]; k < ic[ll+1]; k++){
1700 if (mask[jc[k]] < id[i]){
1701 mask[jc[k]] = nz;
1702 jd[nz] = jc[k];
1703 d[nz] += a[j]*b[l]*c[k];
1704 nz++;
1705 } else {
1706 assert(jd[mask[jc[k]]] == jc[k]);
1707 d[mask[jc[k]]] += a[j]*b[l]*c[k];
1708 }
1709 }
1710 }
1711 }
1712 id[i+1] = nz;
1713 }
1714 }
1715 break;
1716 case MATRIX_TYPE_PATTERN:
1717 id[0] = 0;
1718 for (i = 0; i < m; i++){
1719 for (j = ia[i]; j < ia[i+1]; j++){
1720 jj = ja[j];
1721 for (l = ib[jj]; l < ib[jj+1]; l++){
1722 ll = jb[l];
1723 for (k = ic[ll]; k < ic[ll+1]; k++){
1724 if (mask[jc[k]] < id[i]){
1725 mask[jc[k]] = nz;
1726 jd[nz] = jc[k];
1727 nz++;
1728 } else {
1729 assert(jd[mask[jc[k]]] == jc[k]);
1730 }
1731 }
1732 }
1733 }
1734 id[i+1] = nz;
1735 }
1736 break;
1737 case MATRIX_TYPE_UNKNOWN:
1738 default:
1739 SparseMatrix_delete(D);
1740 D = NULL; goto RETURN;
1741 break;
1742 }
1743
1744 D->nz = nz;
1745
1746 RETURN:
1747 FREE(mask);
1748 return D;
1749
1750}
1751
1752/* For complex matrix:
1753 if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the x's if {i,j,Round(y)} are the same
1754 if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the y's if {i,j,Round(x)} are the same
1755 so a matrix like {{1,1,2+3i},{1,2,3+4i},{1,1,5+3i},{1,2,4+5i},{1,2,4+4i}} becomes
1756 {{1,1,2+5+3i},{1,2,3+4+4i},{1,2,4+5i}}.
1757
1758 Really this kind of thing is best handled using a 3D sparse matrix like
1759 {{{1,1,3},2},{{1,2,4},3},{{1,1,3},5},{{1,2,4},4}},
1760 which is then aggreted into
1761 {{{1,1,3},2+5},{{1,2,4},3+4},{{1,1,3},5}}
1762 but unfortunately I do not have such object implemented yet.
1763
1764
1765 For other matrix, what_to_sum = SUM_REPEATED_REAL_PART is the same as what_to_sum = SUM_REPEATED_IMAGINARY_PART
1766 or what_to_sum = SUM_REPEATED_ALL. In this implementation we assume that
1767 the {j,y} pairs are dense, so we usea 2D array for hashing
1768*/
1769SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){
1770 /* sum repeated entries in the same row, i.e., {1,1}->1, {1,1}->2 becomes {1,1}->3 */
1771 int *ia = A->ia, *ja = A->ja, type = A->type, n = A->n;
1772 int *mask = NULL, nz = 0, i, j, sta;
1773
1774 if (what_to_sum == SUM_REPEATED_NONE) return A;
1775
1776 mask = MALLOC(sizeof(int)*((size_t)n));
1777 for (i = 0; i < n; i++) mask[i] = -1;
1778
1779 switch (type){
1780 case MATRIX_TYPE_REAL:
1781 {
1782 real *a = (real*) A->a;
1783 nz = 0;
1784 sta = ia[0];
1785 for (i = 0; i < A->m; i++){
1786 for (j = sta; j < ia[i+1]; j++){
1787 if (mask[ja[j]] < ia[i]){
1788 ja[nz] = ja[j];
1789 a[nz] = a[j];
1790 mask[ja[j]] = nz++;
1791 } else {
1792 assert(ja[mask[ja[j]]] == ja[j]);
1793 a[mask[ja[j]]] += a[j];
1794 }
1795 }
1796 sta = ia[i+1];
1797 ia[i+1] = nz;
1798 }
1799 }
1800 break;
1801 case MATRIX_TYPE_COMPLEX:
1802 {
1803 real *a = (real*) A->a;
1804 if (what_to_sum == SUM_REPEATED_ALL){
1805 nz = 0;
1806 sta = ia[0];
1807 for (i = 0; i < A->m; i++){
1808 for (j = sta; j < ia[i+1]; j++){
1809 if (mask[ja[j]] < ia[i]){
1810 ja[nz] = ja[j];
1811 a[2*nz] = a[2*j];
1812 a[2*nz+1] = a[2*j+1];
1813 mask[ja[j]] = nz++;
1814 } else {
1815 assert(ja[mask[ja[j]]] == ja[j]);
1816 a[2*mask[ja[j]]] += a[2*j];
1817 a[2*mask[ja[j]]+1] += a[2*j+1];
1818 }
1819 }
1820 sta = ia[i+1];
1821 ia[i+1] = nz;
1822 }
1823 } else if (what_to_sum == SUM_IMGINARY_KEEP_LAST_REAL){
1824 /* merge {i,j,R1,I1} and {i,j,R2,I2} into {i,j,R1+R2,I2}*/
1825 nz = 0;
1826 sta = ia[0];
1827 for (i = 0; i < A->m; i++){
1828 for (j = sta; j < ia[i+1]; j++){
1829 if (mask[ja[j]] < ia[i]){
1830 ja[nz] = ja[j];
1831 a[2*nz] = a[2*j];
1832 a[2*nz+1] = a[2*j+1];
1833 mask[ja[j]] = nz++;
1834 } else {
1835 assert(ja[mask[ja[j]]] == ja[j]);
1836 a[2*mask[ja[j]]] += a[2*j];
1837 a[2*mask[ja[j]]+1] = a[2*j+1];
1838 }
1839 }
1840 sta = ia[i+1];
1841 ia[i+1] = nz;
1842 }
1843 } else if (what_to_sum == SUM_REPEATED_REAL_PART){
1844 int ymin, ymax, id;
1845 ymax = ymin = a[1];
1846 nz = 0;
1847 for (i = 0; i < A->m; i++){
1848 for (j = ia[i]; j < ia[i+1]; j++){
1849 ymax = MAX(ymax, (int) a[2*nz+1]);
1850 ymin = MIN(ymin, (int) a[2*nz+1]);
1851 nz++;
1852 }
1853 }
1854 FREE(mask);
1855 mask = MALLOC(sizeof(int)*((size_t)n)*((size_t)(ymax-ymin+1)));
1856 for (i = 0; i < n*(ymax-ymin+1); i++) mask[i] = -1;
1857
1858 nz = 0;
1859 sta = ia[0];
1860 for (i = 0; i < A->m; i++){
1861 for (j = sta; j < ia[i+1]; j++){
1862 id = ja[j] + ((int)a[2*j+1] - ymin)*n;
1863 if (mask[id] < ia[i]){
1864 ja[nz] = ja[j];
1865 a[2*nz] = a[2*j];
1866 a[2*nz+1] = a[2*j+1];
1867 mask[id] = nz++;
1868 } else {
1869 assert(id < n*(ymax-ymin+1));
1870 assert(ja[mask[id]] == ja[j]);
1871 a[2*mask[id]] += a[2*j];
1872 a[2*mask[id]+1] = a[2*j+1];
1873 }
1874 }
1875 sta = ia[i+1];
1876 ia[i+1] = nz;
1877 }
1878
1879 } else if (what_to_sum == SUM_REPEATED_IMAGINARY_PART){
1880 int xmin, xmax, id;
1881 xmax = xmin = a[1];
1882 nz = 0;
1883 for (i = 0; i < A->m; i++){
1884 for (j = ia[i]; j < ia[i+1]; j++){
1885 xmax = MAX(xmax, (int) a[2*nz]);
1886 xmin = MAX(xmin, (int) a[2*nz]);
1887 nz++;
1888 }
1889 }
1890 FREE(mask);
1891 mask = MALLOC(sizeof(int)*((size_t)n)*((size_t)(xmax-xmin+1)));
1892 for (i = 0; i < n*(xmax-xmin+1); i++) mask[i] = -1;
1893
1894 nz = 0;
1895 sta = ia[0];
1896 for (i = 0; i < A->m; i++){
1897 for (j = sta; j < ia[i+1]; j++){
1898 id = ja[j] + ((int)a[2*j] - xmin)*n;
1899 if (mask[id] < ia[i]){
1900 ja[nz] = ja[j];
1901 a[2*nz] = a[2*j];
1902 a[2*nz+1] = a[2*j+1];
1903 mask[id] = nz++;
1904 } else {
1905 assert(ja[mask[id]] == ja[j]);
1906 a[2*mask[id]] = a[2*j];
1907 a[2*mask[id]+1] += a[2*j+1];
1908 }
1909 }
1910 sta = ia[i+1];
1911 ia[i+1] = nz;
1912 }
1913
1914 }
1915 }
1916 break;
1917 case MATRIX_TYPE_INTEGER:
1918 {
1919 int *a = (int*) A->a;
1920 nz = 0;
1921 sta = ia[0];
1922 for (i = 0; i < A->m; i++){
1923 for (j = sta; j < ia[i+1]; j++){
1924 if (mask[ja[j]] < ia[i]){
1925 ja[nz] = ja[j];
1926 a[nz] = a[j];
1927 mask[ja[j]] = nz++;
1928 } else {
1929 assert(ja[mask[ja[j]]] == ja[j]);
1930 a[mask[ja[j]]] += a[j];
1931 }
1932 }
1933 sta = ia[i+1];
1934 ia[i+1] = nz;
1935 }
1936 }
1937 break;
1938 case MATRIX_TYPE_PATTERN:
1939 {
1940 nz = 0;
1941 sta = ia[0];
1942 for (i = 0; i < A->m; i++){
1943 for (j = sta; j < ia[i+1]; j++){
1944 if (mask[ja[j]] < ia[i]){
1945 ja[nz] = ja[j];
1946 mask[ja[j]] = nz++;
1947 } else {
1948 assert(ja[mask[ja[j]]] == ja[j]);
1949 }
1950 }
1951 sta = ia[i+1];
1952 ia[i+1] = nz;
1953 }
1954 }
1955 break;
1956 case MATRIX_TYPE_UNKNOWN:
1957 return NULL;
1958 break;
1959 default:
1960 return NULL;
1961 break;
1962 }
1963 A->nz = nz;
1964 FREE(mask);
1965 return A;
1966}
1967
1968SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, int nentries, int *irn, int *jcn, void *val){
1969 int nz, nzmax, i;
1970
1971 assert(A->format == FORMAT_COORD);
1972 if (nentries <= 0) return A;
1973 nz = A->nz;
1974 nzmax = A->nzmax;
1975
1976 if (nz + nentries >= A->nzmax){
1977 nzmax = nz + nentries;
1978 nzmax = MAX(10, (int) 0.2*nzmax) + nzmax;
1979 A = SparseMatrix_realloc(A, nzmax);
1980 }
1981 MEMCPY((char*) A->ia + ((size_t)nz)*sizeof(int)/sizeof(char), irn, sizeof(int)*((size_t)nentries));
1982 MEMCPY((char*) A->ja + ((size_t)nz)*sizeof(int)/sizeof(char), jcn, sizeof(int)*((size_t)nentries));
1983 if (A->size) MEMCPY((char*) A->a + ((size_t)nz)*A->size/sizeof(char), val, A->size*((size_t)nentries));
1984 for (i = 0; i < nentries; i++) {
1985 if (irn[i] >= A->m) A->m = irn[i]+1;
1986 if (jcn[i] >= A->n) A->n = jcn[i]+1;
1987 }
1988 A->nz += nentries;
1989 return A;
1990}
1991
1992
1993SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A){
1994 int i, j, *ia, *ja, nz, sta;
1995
1996 if (!A) return A;
1997
1998 nz = 0;
1999 ia = A->ia;
2000 ja = A->ja;
2001 sta = ia[0];
2002 switch (A->type){
2003 case MATRIX_TYPE_REAL:{
2004 real *a = (real*) A->a;
2005 for (i = 0; i < A->m; i++){
2006 for (j = sta; j < ia[i+1]; j++){
2007 if (ja[j] != i){
2008 ja[nz] = ja[j];
2009 a[nz++] = a[j];
2010 }
2011 }
2012 sta = ia[i+1];
2013 ia[i+1] = nz;
2014 }
2015 A->nz = nz;
2016 break;
2017 }
2018 case MATRIX_TYPE_COMPLEX:{
2019 real *a = (real*) A->a;
2020 for (i = 0; i < A->m; i++){
2021 for (j = sta; j < ia[i+1]; j++){
2022 if (ja[j] != i){
2023 ja[nz] = ja[j];
2024 a[2*nz] = a[2*j];
2025 a[2*nz+1] = a[2*j+1];
2026 nz++;
2027 }
2028 }
2029 sta = ia[i+1];
2030 ia[i+1] = nz;
2031 }
2032 A->nz = nz;
2033 break;
2034 }
2035 case MATRIX_TYPE_INTEGER:{
2036 int *a = (int*) A->a;
2037 for (i = 0; i < A->m; i++){
2038 for (j = sta; j < ia[i+1]; j++){
2039 if (ja[j] != i){
2040 ja[nz] = ja[j];
2041 a[nz++] = a[j];
2042 }
2043 }
2044 sta = ia[i+1];
2045 ia[i+1] = nz;
2046 }
2047 A->nz = nz;
2048 break;
2049 }
2050 case MATRIX_TYPE_PATTERN:{
2051 for (i = 0; i < A->m; i++){
2052 for (j = sta; j < ia[i+1]; j++){
2053 if (ja[j] != i){
2054 ja[nz++] = ja[j];
2055 }
2056 }
2057 sta = ia[i+1];
2058 ia[i+1] = nz;
2059 }
2060 A->nz = nz;
2061 break;
2062 }
2063 case MATRIX_TYPE_UNKNOWN:
2064 return NULL;
2065 default:
2066 return NULL;
2067 }
2068
2069 return A;
2070}
2071
2072
2073SparseMatrix SparseMatrix_remove_upper(SparseMatrix A){/* remove diag and upper diag */
2074 int i, j, *ia, *ja, nz, sta;
2075
2076 if (!A) return A;
2077
2078 nz = 0;
2079 ia = A->ia;
2080 ja = A->ja;
2081 sta = ia[0];
2082 switch (A->type){
2083 case MATRIX_TYPE_REAL:{
2084 real *a = (real*) A->a;
2085 for (i = 0; i < A->m; i++){
2086 for (j = sta; j < ia[i+1]; j++){
2087 if (ja[j] < i){
2088 ja[nz] = ja[j];
2089 a[nz++] = a[j];
2090 }
2091 }
2092 sta = ia[i+1];
2093 ia[i+1] = nz;
2094 }
2095 A->nz = nz;
2096 break;
2097 }
2098 case MATRIX_TYPE_COMPLEX:{
2099 real *a = (real*) A->a;
2100 for (i = 0; i < A->m; i++){
2101 for (j = sta; j < ia[i+1]; j++){
2102 if (ja[j] < i){
2103 ja[nz] = ja[j];
2104 a[2*nz] = a[2*j];
2105 a[2*nz+1] = a[2*j+1];
2106 nz++;
2107 }
2108 }
2109 sta = ia[i+1];
2110 ia[i+1] = nz;
2111 }
2112 A->nz = nz;
2113 break;
2114 }
2115 case MATRIX_TYPE_INTEGER:{
2116 int *a = (int*) A->a;
2117 for (i = 0; i < A->m; i++){
2118 for (j = sta; j < ia[i+1]; j++){
2119 if (ja[j] < i){
2120 ja[nz] = ja[j];
2121 a[nz++] = a[j];
2122 }
2123 }
2124 sta = ia[i+1];
2125 ia[i+1] = nz;
2126 }
2127 A->nz = nz;
2128 break;
2129 }
2130 case MATRIX_TYPE_PATTERN:{
2131 for (i = 0; i < A->m; i++){
2132 for (j = sta; j < ia[i+1]; j++){
2133 if (ja[j] < i){
2134 ja[nz++] = ja[j];
2135 }
2136 }
2137 sta = ia[i+1];
2138 ia[i+1] = nz;
2139 }
2140 A->nz = nz;
2141 break;
2142 }
2143 case MATRIX_TYPE_UNKNOWN:
2144 return NULL;
2145 default:
2146 return NULL;
2147 }
2148
2149 clear_flag(A->property, MATRIX_PATTERN_SYMMETRIC);
2150 clear_flag(A->property, MATRIX_SYMMETRIC);
2151 clear_flag(A->property, MATRIX_SKEW);
2152 clear_flag(A->property, MATRIX_HERMITIAN);
2153 return A;
2154}
2155
2156
2157
2158
2159SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A){
2160 int i, j, *ia, *ja;
2161 real deg;
2162
2163 if (!A) return A;
2164
2165 ia = A->ia;
2166 ja = A->ja;
2167 switch (A->type){
2168 case MATRIX_TYPE_REAL:{
2169 real *a = (real*) A->a;
2170 for (i = 0; i < A->m; i++){
2171 deg = ia[i+1] - ia[i];
2172 for (j = ia[i]; j < ia[i+1]; j++){
2173 a[j] = a[j]/deg;
2174 }
2175 }
2176 break;
2177 }
2178 case MATRIX_TYPE_COMPLEX:{
2179 real *a = (real*) A->a;
2180 for (i = 0; i < A->m; i++){
2181 deg = ia[i+1] - ia[i];
2182 for (j = ia[i]; j < ia[i+1]; j++){
2183 if (ja[j] != i){
2184 a[2*j] = a[2*j]/deg;
2185 a[2*j+1] = a[2*j+1]/deg;
2186 }
2187 }
2188 }
2189 break;
2190 }
2191 case MATRIX_TYPE_INTEGER:{
2192 assert(0);/* this operation would not make sense for int matrix */
2193 break;
2194 }
2195 case MATRIX_TYPE_PATTERN:{
2196 break;
2197 }
2198 case MATRIX_TYPE_UNKNOWN:
2199 return NULL;
2200 default:
2201 return NULL;
2202 }
2203
2204 return A;
2205}
2206
2207
2208SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A){
2209 /* symmetric, all entries to 1, diaginal removed */
2210 int i, *ia, *ja, nz, m, n;
2211 real *a;
2212 SparseMatrix B;
2213
2214 if (!A) return A;
2215
2216 nz = A->nz;
2217 ia = A->ia;
2218 ja = A->ja;
2219 n = A->n;
2220 m = A->m;
2221
2222 if (n != m) return NULL;
2223
2224 B = SparseMatrix_new(m, n, nz, MATRIX_TYPE_PATTERN, FORMAT_CSR);
2225
2226 MEMCPY(B->ia, ia, sizeof(int)*((size_t)(m+1)));
2227 MEMCPY(B->ja, ja, sizeof(int)*((size_t)nz));
2228 B->nz = A->nz;
2229
2230 A = SparseMatrix_symmetrize(B, TRUE);
2231 SparseMatrix_delete(B);
2232 A = SparseMatrix_remove_diagonal(A);
2233 A->a = MALLOC(sizeof(real)*((size_t)(A->nz)));
2234 a = (real*) A->a;
2235 for (i = 0; i < A->nz; i++) a[i] = 1.;
2236 A->type = MATRIX_TYPE_REAL;
2237 A->size = sizeof(real);
2238 return A;
2239}
2240
2241
2242
2243SparseMatrix SparseMatrix_normalize_to_rowsum1(SparseMatrix A){
2244 int i, j;
2245 real sum, *a;
2246
2247 if (!A) return A;
2248 if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) {
2249#ifdef DEBUG
2250 printf("only CSR and real matrix supported.\n");
2251#endif
2252 return A;
2253 }
2254
2255 a = (real*) A->a;
2256 for (i = 0; i < A->m; i++){
2257 sum = 0;
2258 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2259 sum += a[j];
2260 }
2261 if (sum != 0){
2262 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2263 a[j] /= sum;
2264 }
2265 }
2266 }
2267 return A;
2268}
2269
2270
2271
2272SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix A){
2273 int i, j;
2274 real max, *a;
2275
2276 if (!A) return A;
2277 if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) {
2278#ifdef DEBUG
2279 printf("only CSR and real matrix supported.\n");
2280#endif
2281 return A;
2282 }
2283
2284 a = (real*) A->a;
2285 for (i = 0; i < A->m; i++){
2286 max = 0;
2287 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2288 max = MAX(ABS(a[j]), max);
2289 }
2290 if (max != 0){
2291 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2292 a[j] /= max;
2293 }
2294 }
2295 }
2296 return A;
2297}
2298
2299
2300SparseMatrix SparseMatrix_to_complex(SparseMatrix A){
2301 int i, *ia, *ja;
2302
2303 if (!A) return A;
2304 if (A->format != FORMAT_CSR) {
2305#ifdef DEBUG
2306 printf("only CSR format supported.\n");
2307#endif
2308 return A;
2309 }
2310
2311 ia = A->ia;
2312 ja = A->ja;
2313 switch (A->type){
2314 case MATRIX_TYPE_REAL:{
2315 real *a = (real*) A->a;
2316 int nz = A->nz;
2317 A->a = a = REALLOC(a, 2*sizeof(real)*nz);
2318 for (i = nz - 1; i >= 0; i--){
2319 a[2*i] = a[i];
2320 a[2*i - 1] = 0;
2321 }
2322 A->type = MATRIX_TYPE_COMPLEX;
2323 A->size = 2*sizeof(real);
2324 break;
2325 }
2326 case MATRIX_TYPE_COMPLEX:{
2327 break;
2328 }
2329 case MATRIX_TYPE_INTEGER:{
2330 int *a = (int*) A->a;
2331 int nz = A->nz;
2332 real *aa = A->a = MALLOC(2*sizeof(real)*nz);
2333 for (i = nz - 1; i >= 0; i--){
2334 aa[2*i] = (real) a[i];
2335 aa[2*i - 1] = 0;
2336 }
2337 A->type = MATRIX_TYPE_COMPLEX;
2338 A->size = 2*sizeof(real);
2339 FREE(a);
2340 break;
2341 }
2342 case MATRIX_TYPE_PATTERN:{
2343 break;
2344 }
2345 case MATRIX_TYPE_UNKNOWN:
2346 return NULL;
2347 default:
2348 return NULL;
2349 }
2350
2351 return A;
2352}
2353
2354
2355SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun)(double x)){
2356 int i, j;
2357 real *a;
2358
2359
2360 if (!A) return A;
2361 if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) {
2362#ifdef DEBUG
2363 printf("only CSR and real matrix supported.\n");
2364#endif
2365 return A;
2366 }
2367
2368
2369 a = (real*) A->a;
2370 for (i = 0; i < A->m; i++){
2371 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2372 a[j] = fun(a[j]);
2373 }
2374 }
2375 return A;
2376}
2377
2378SparseMatrix SparseMatrix_apply_fun_general(SparseMatrix A, void (*fun)(int i, int j, int n, double *x)){
2379 int i, j;
2380 real *a;
2381 int len = 1;
2382
2383 if (!A) return A;
2384 if (A->format != FORMAT_CSR || (A->type != MATRIX_TYPE_REAL&&A->type != MATRIX_TYPE_COMPLEX)) {
2385#ifdef DEBUG
2386 printf("SparseMatrix_apply_fun: only CSR and real/complex matrix supported.\n");
2387#endif
2388 return A;
2389 }
2390 if (A->type == MATRIX_TYPE_COMPLEX) len = 2;
2391
2392 a = (real*) A->a;
2393 for (i = 0; i < A->m; i++){
2394 for (j = A->ia[i]; j < A->ia[i+1]; j++){
2395 fun(i, A->ja[j], len, &a[len*j]);
2396 }
2397 }
2398 return A;
2399}
2400
2401
2402SparseMatrix SparseMatrix_crop(SparseMatrix A, real epsilon){
2403 int i, j, *ia, *ja, nz, sta;
2404
2405 if (!A) return A;
2406
2407 nz = 0;
2408 ia = A->ia;
2409 ja = A->ja;
2410 sta = ia[0];
2411 switch (A->type){
2412 case MATRIX_TYPE_REAL:{
2413 real *a = (real*) A->a;
2414 for (i = 0; i < A->m; i++){
2415 for (j = sta; j < ia[i+1]; j++){
2416 if (ABS(a[j]) > epsilon){
2417 ja[nz] = ja[j];
2418 a[nz++] = a[j];
2419 }
2420 }
2421 sta = ia[i+1];
2422 ia[i+1] = nz;
2423 }
2424 A->nz = nz;
2425 break;
2426 }
2427 case MATRIX_TYPE_COMPLEX:{
2428 real *a = (real*) A->a;
2429 for (i = 0; i < A->m; i++){
2430 for (j = sta; j < ia[i+1]; j++){
2431 if (sqrt(a[2*j]*a[2*j]+a[2*j+1]*a[2*j+1]) > epsilon){
2432 ja[nz] = ja[j];
2433 a[2*nz] = a[2*j];
2434 a[2*nz+1] = a[2*j+1];
2435 nz++;
2436 }
2437 }
2438 sta = ia[i+1];
2439 ia[i+1] = nz;
2440 }
2441 A->nz = nz;
2442 break;
2443 }
2444 case MATRIX_TYPE_INTEGER:{
2445 int *a = (int*) A->a;
2446 for (i = 0; i < A->m; i++){
2447 for (j = sta; j < ia[i+1]; j++){
2448 if (ABS(a[j]) > epsilon){
2449 ja[nz] = ja[j];
2450 a[nz++] = a[j];
2451 }
2452 }
2453 sta = ia[i+1];
2454 ia[i+1] = nz;
2455 }
2456 A->nz = nz;
2457 break;
2458 }
2459 case MATRIX_TYPE_PATTERN:{
2460 break;
2461 }
2462 case MATRIX_TYPE_UNKNOWN:
2463 return NULL;
2464 default:
2465 return NULL;
2466 }
2467
2468 return A;
2469}
2470
2471SparseMatrix SparseMatrix_copy(SparseMatrix A){
2472 SparseMatrix B;
2473 if (!A) return A;
2474 B = SparseMatrix_general_new(A->m, A->n, A->nz, A->type, A->size, A->format);
2475 MEMCPY(B->ia, A->ia, sizeof(int)*((size_t)(A->m+1)));
2476 MEMCPY(B->ja, A->ja, sizeof(int)*((size_t)(A->ia[A->m])));
2477 if (A->a) MEMCPY(B->a, A->a, A->size*((size_t)A->nz));
2478 B->property = A->property;
2479 B->nz = A->nz;
2480 return B;
2481}
2482
2483int SparseMatrix_has_diagonal(SparseMatrix A){
2484
2485 int i, j, m = A->m, *ia = A->ia, *ja = A->ja;
2486
2487 for (i = 0; i < m; i++){
2488 for (j = ia[i]; j < ia[i+1]; j++){
2489 if (i == ja[j]) return TRUE;
2490 }
2491 }
2492 return FALSE;
2493}
2494
2495void SparseMatrix_level_sets_internal(int khops, SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask){
2496 /* mask is assumed to be initialized to negative if provided.
2497 . On exit, mask = levels for visited nodes (1 for root, 2 for its neighbors, etc),
2498 . unless reinitialize_mask = TRUE, in which case mask = -1.
2499 khops: max number of hops allowed. If khops < 0, no limit is applied.
2500 A: the graph, undirected
2501 root: starting node
2502 nlevel: max distance to root from any node (in the connected comp)
2503 levelset_ptr, levelset: the level sets
2504 */
2505 int i, j, sta = 0, sto = 1, nz, ii;
2506 int m = A->m, *ia = A->ia, *ja = A->ja;
2507
2508 if (!(*levelset_ptr)) *levelset_ptr = MALLOC(sizeof(int)*((size_t)(m+2)));
2509 if (!(*levelset)) *levelset = MALLOC(sizeof(int)*((size_t)m));
2510 if (!(*mask)) {
2511 *mask = malloc(sizeof(int)*((size_t)m));
2512 for (i = 0; i < m; i++) (*mask)[i] = UNMASKED;
2513 }
2514
2515 *nlevel = 0;
2516 assert(root >= 0 && root < m);
2517 (*levelset_ptr)[0] = 0;
2518 (*levelset_ptr)[1] = 1;
2519 (*levelset)[0] = root;
2520 (*mask)[root] = 1;
2521 *nlevel = 1;
2522 nz = 1;
2523 sta = 0; sto = 1;
2524 while (sto > sta && (khops < 0 || *nlevel <= khops)){
2525 for (i = sta; i < sto; i++){
2526 ii = (*levelset)[i];
2527 for (j = ia[ii]; j < ia[ii+1]; j++){
2528 if (ii == ja[j]) continue;
2529 if ((*mask)[ja[j]] < 0){
2530 (*levelset)[nz++] = ja[j];
2531 (*mask)[ja[j]] = *nlevel + 1;
2532 }
2533 }
2534 }
2535 (*levelset_ptr)[++(*nlevel)] = nz;
2536 sta = sto;
2537 sto = nz;
2538 }
2539 if (khops < 0 || *nlevel <= khops){
2540 (*nlevel)--;
2541 }
2542 if (reinitialize_mask) for (i = 0; i < (*levelset_ptr)[*nlevel]; i++) (*mask)[(*levelset)[i]] = UNMASKED;
2543}
2544
2545void SparseMatrix_level_sets(SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask){
2546
2547 int khops = -1;
2548
2549 return SparseMatrix_level_sets_internal(khops, A, root, nlevel, levelset_ptr, levelset, mask, reinitialize_mask);
2550
2551}
2552void SparseMatrix_level_sets_khops(int khops, SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask){
2553
2554 return SparseMatrix_level_sets_internal(khops, A, root, nlevel, levelset_ptr, levelset, mask, reinitialize_mask);
2555
2556}
2557
2558void SparseMatrix_weakly_connected_components(SparseMatrix A0, int *ncomp, int **comps, int **comps_ptr){
2559 SparseMatrix A = A0;
2560 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL, nlevel;
2561 int m = A->m, i, nn;
2562
2563 if (!SparseMatrix_is_symmetric(A, TRUE)){
2564 A = SparseMatrix_symmetrize(A, TRUE);
2565 }
2566 if (!(*comps_ptr)) *comps_ptr = MALLOC(sizeof(int)*((size_t)(m+1)));
2567
2568 *ncomp = 0;
2569 (*comps_ptr)[0] = 0;
2570 for (i = 0; i < m; i++){
2571 if (i == 0 || mask[i] < 0) {
2572 SparseMatrix_level_sets(A, i, &nlevel, &levelset_ptr, &levelset, &mask, FALSE);
2573 if (i == 0) *comps = levelset;
2574 nn = levelset_ptr[nlevel];
2575 levelset += nn;
2576 (*comps_ptr)[(*ncomp)+1] = (*comps_ptr)[(*ncomp)] + nn;
2577 (*ncomp)++;
2578 }
2579
2580 }
2581 if (A != A0) SparseMatrix_delete(A);
2582 if (levelset_ptr) FREE(levelset_ptr);
2583
2584 FREE(mask);
2585}
2586
2587
2588
2589struct nodedata_struct {
2590 real dist;/* distance to root */
2591 int id;/*node id */
2592};
2593typedef struct nodedata_struct* nodedata;
2594
2595static int cmp(void*i, void*j){
2596 nodedata d1, d2;
2597
2598 d1 = (nodedata) i;
2599 d2 = (nodedata) j;
2600 if (d1->dist > d2->dist){
2601 return 1;
2602 } else if (d1->dist == d2->dist){
2603 return 0;
2604 } else {
2605 return -1;
2606 }
2607}
2608
2609static int Dijkstra_internal(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max, int *mask){
2610 /* Find the shortest path distance of all nodes to root. If khops >= 0, the shortest ath is of distance <= khops,
2611
2612 A: the nxn connectivity matrix. Entries are assumed to be nonnegative. Absolute value will be taken if
2613 . entry value is negative.
2614 dist: length n. On on exit contain the distance from root to every other node. dist[root] = 0. dist[i] = distance from root to node i.
2615 . if the graph is disconnetced, unreachable node have a distance -1.
2616 . note: ||root - list[i]|| =!= dist[i] !!!, instead, ||root - list[i]|| == dist[list[i]]
2617 nlist: number of nodes visited
2618 list: length n. the list of node in order of their extraction from the heap.
2619 . The distance from root to last in the list should be the maximum
2620 dist_max: the maximum distance, should be realized at node list[nlist-1].
2621 mask: if NULL, not used. Othewise, only nodes i with mask[i] > 0 will be considered
2622 return: 0 if every node is reachable. -1 if not */
2623
2624 int m = A->m, i, j, jj, *ia = A->ia, *ja = A->ja, heap_id;
2625 BinaryHeap h;
2626 real *a = NULL, *aa;
2627 int *ai;
2628 nodedata ndata, ndata_min;
2629 enum {UNVISITED = -2, FINISHED = -1};
2630 int *heap_ids; /* node ID to heap ID array. Initialised to UNVISITED.
2631 Set to FINISHED after extracted as min from heap */
2632 int found = 0;
2633
2634 assert(SparseMatrix_is_symmetric(A, TRUE));
2635
2636 assert(m == A->n);
2637
2638 switch (A->type){
2639 case MATRIX_TYPE_COMPLEX:
2640 aa = (real*) A->a;
2641 a = MALLOC(sizeof(real)*((size_t)(A->nz)));
2642 for (i = 0; i < A->nz; i++) a[i] = aa[i*2];
2643 break;
2644 case MATRIX_TYPE_REAL:
2645 a = (real*) A->a;
2646 break;
2647 case MATRIX_TYPE_INTEGER:
2648 ai = (int*) A->a;
2649 a = MALLOC(sizeof(real)*((size_t)(A->nz)));
2650 for (i = 0; i < A->nz; i++) a[i] = (real) ai[i];
2651 break;
2652 case MATRIX_TYPE_PATTERN:
2653 a = MALLOC(sizeof(real)*((size_t)A->nz));
2654 for (i = 0; i < A->nz; i++) a[i] = 1.;
2655 break;
2656 default:
2657 assert(0);/* no such matrix type */
2658 }
2659
2660 heap_ids = MALLOC(sizeof(int)*((size_t)m));
2661 for (i = 0; i < m; i++) {
2662 dist[i] = -1;
2663 heap_ids[i] = UNVISITED;
2664 }
2665
2666 h = BinaryHeap_new(cmp);
2667 assert(h);
2668
2669 /* add root as the first item in the heap */
2670 ndata = MALLOC(sizeof(struct nodedata_struct));
2671 ndata->dist = 0;
2672 ndata->id = root;
2673 heap_ids[root] = BinaryHeap_insert(h, ndata);
2674
2675 assert(heap_ids[root] >= 0);/* by design heap ID from BinaryHeap_insert >=0*/
2676
2677 while ((ndata_min = BinaryHeap_extract_min(h))){
2678 i = ndata_min->id;
2679 dist[i] = ndata_min->dist;
2680 list[found++] = i;
2681 heap_ids[i] = FINISHED;
2682 //fprintf(stderr," =================\n min extracted is id=%d, dist=%f\n",i, ndata_min->dist);
2683 for (j = ia[i]; j < ia[i+1]; j++){
2684 jj = ja[j];
2685 heap_id = heap_ids[jj];
2686
2687 if (jj == i || heap_id == FINISHED || (mask && mask[jj] < 0)) continue;
2688
2689 if (heap_id == UNVISITED){
2690 ndata = MALLOC(sizeof(struct nodedata_struct));
2691 ndata->dist = ABS(a[j]) + ndata_min->dist;
2692 ndata->id = jj;
2693 heap_ids[jj] = BinaryHeap_insert(h, ndata);
2694 //fprintf(stderr," set neighbor id=%d, dist=%f, hid = %d, a[%d]=%f, dist=%f\n",jj, ndata->dist, heap_ids[jj], jj, a[j], ndata->dist);
2695
2696 } else {
2697 ndata = BinaryHeap_get_item(h, heap_id);
2698 ndata->dist = MIN(ndata->dist, ABS(a[j]) + ndata_min->dist);
2699 assert(ndata->id == jj);
2700 BinaryHeap_reset(h, heap_id, ndata);
2701 //fprintf(stderr," reset neighbor id=%d, dist=%f, hid = %d, a[%d]=%f, dist=%f\n",jj, ndata->dist,heap_id, jj, a[j], ndata->dist);
2702 }
2703 }
2704 FREE(ndata_min);
2705 }
2706 *nlist = found;
2707 *dist_max = dist[i];
2708
2709
2710 BinaryHeap_delete(h, FREE);
2711 FREE(heap_ids);
2712 if (a && a != A->a) FREE(a);
2713 if (found == m || mask){
2714 return 0;
2715 } else {
2716 return -1;
2717 }
2718}
2719
2720static int Dijkstra(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max){
2721 return Dijkstra_internal(A, root, dist, nlist, list, dist_max, NULL);
2722}
2723
2724static int Dijkstra_masked(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max, int *mask){
2725 /* this makes the algorithm only consider nodes that are masked.
2726 nodes are masked as 1, 2, ..., mask_max, which is (the number of hops from root)+1.
2727 Only paths consists of nodes that are masked are allowed. */
2728
2729 return Dijkstra_internal(A, root, dist, nlist, list, dist_max, mask);
2730}
2731
2732real SparseMatrix_pseudo_diameter_weighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ){
2733 /* weighted graph. But still assume to be undirected. unsymmetric matrix ill be symmetrized */
2734 SparseMatrix A = A0;
2735 int m = A->m, i, *list = NULL, nlist;
2736 int flag;
2737 real *dist = NULL, dist_max = -1, dist0 = -1;
2738 int roots[5], iroots, end11, end22;
2739
2740 if (!SparseMatrix_is_symmetric(A, TRUE)){
2741 A = SparseMatrix_symmetrize(A, TRUE);
2742 }
2743 assert(m == A->n);
2744
2745 dist = MALLOC(sizeof(real)*((size_t)m));
2746 list = MALLOC(sizeof(int)*((size_t)m));
2747 nlist = 1;
2748 list[nlist-1] = root;
2749
2750 assert(SparseMatrix_is_symmetric(A, TRUE));
2751
2752 do {
2753 dist0 = dist_max;
2754 root = list[nlist - 1];
2755 flag = Dijkstra(A, root, dist, &nlist, list, &dist_max);
2756 //fprintf(stderr,"after Dijkstra, {%d,%d}-%f\n",root, list[nlist-1], dist_max);
2757 assert(dist[list[nlist-1]] == dist_max);
2758 assert(root == list[0]);
2759 assert(nlist > 0);
2760 } while (dist_max > dist0);
2761
2762 *connectedQ = (!flag);
2763 assert((dist_max - dist0)/MAX(1, MAX(ABS(dist0), ABS(dist_max))) < 1.e-10);
2764
2765 *end1 = root;
2766 *end2 = list[nlist-1];
2767 //fprintf(stderr,"after search for diameter, diam = %f, ends = {%d,%d}\n", dist_max, *end1, *end2);
2768
2769 if (aggressive){
2770 iroots = 0;
2771 for (i = MAX(0, nlist - 6); i < nlist - 1; i++){
2772 roots[iroots++] = list[i];
2773 }
2774 for (i = 0; i < iroots; i++){
2775 root = roots[i];
2776 dist0 = dist_max;
2777 fprintf(stderr,"search for diameter again from root=%d\n", root);
2778 dist_max = SparseMatrix_pseudo_diameter_weighted(A, root, FALSE, &end11, &end22, connectedQ);
2779 if (dist_max > dist0){
2780 *end1 = end11; *end2 = end22;
2781 } else {
2782 dist_max = dist0;
2783 }
2784 }
2785 fprintf(stderr,"after aggressive search for diameter, diam = %f, ends = {%d,%d}\n", dist_max, *end1, *end2);
2786 }
2787
2788 FREE(dist);
2789 FREE(list);
2790
2791 if (A != A0) SparseMatrix_delete(A);
2792 return dist_max;
2793
2794}
2795
2796real SparseMatrix_pseudo_diameter_unweighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ){
2797 /* assume unit edge length! unsymmetric matrix ill be symmetrized */
2798 SparseMatrix A = A0;
2799 int m = A->m, i;
2800 int nlevel;
2801 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
2802 int nlevel0 = 0;
2803 int roots[5], iroots, enda, endb;
2804
2805 if (!SparseMatrix_is_symmetric(A, TRUE)){
2806 A = SparseMatrix_symmetrize(A, TRUE);
2807 }
2808
2809 assert(SparseMatrix_is_symmetric(A, TRUE));
2810
2811 SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
2812 // fprintf(stderr,"after level set, {%d,%d}=%d\n",levelset[0], levelset[levelset_ptr[nlevel]-1], nlevel);
2813
2814 *connectedQ = (levelset_ptr[nlevel] == m);
2815 while (nlevel0 < nlevel){
2816 nlevel0 = nlevel;
2817 root = levelset[levelset_ptr[nlevel] - 1];
2818 SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
2819 //fprintf(stderr,"after level set, {%d,%d}=%d\n",levelset[0], levelset[levelset_ptr[nlevel]-1], nlevel);
2820 }
2821 *end1 = levelset[0];
2822 *end2 = levelset[levelset_ptr[nlevel]-1];
2823
2824 if (aggressive){
2825 nlevel0 = nlevel;
2826 iroots = 0;
2827 for (i = levelset_ptr[nlevel-1]; i < MIN(levelset_ptr[nlevel], levelset_ptr[nlevel-1]+5); i++){
2828 iroots++;
2829 roots[i - levelset_ptr[nlevel-1]] = levelset[i];
2830 }
2831 for (i = 0; i < iroots; i++){
2832 root = roots[i];
2833 nlevel = (int) SparseMatrix_pseudo_diameter_unweighted(A, root, FALSE, &enda, &endb, connectedQ);
2834 if (nlevel > nlevel0) {
2835 nlevel0 = nlevel;
2836 *end1 = enda;
2837 *end2 = endb;
2838 }
2839 }
2840 }
2841
2842 FREE(levelset_ptr);
2843 FREE(levelset);
2844 FREE(mask);
2845 if (A != A0) SparseMatrix_delete(A);
2846 return (real) nlevel0 - 1;
2847}
2848
2849real SparseMatrix_pseudo_diameter_only(SparseMatrix A){
2850 int end1, end2, connectedQ;
2851 return SparseMatrix_pseudo_diameter_unweighted(A, 0, FALSE, &end1, &end2, &connectedQ);
2852}
2853
2854int SparseMatrix_connectedQ(SparseMatrix A0){
2855 int root = 0, nlevel, *levelset_ptr = NULL, *levelset = NULL, *mask = NULL, connected;
2856 SparseMatrix A = A0;
2857
2858 if (!SparseMatrix_is_symmetric(A, TRUE)){
2859 if (A->m != A->n) return FALSE;
2860 A = SparseMatrix_symmetrize(A, TRUE);
2861 }
2862
2863 SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
2864 connected = (levelset_ptr[nlevel] == A->m);
2865
2866 FREE(levelset_ptr);
2867 FREE(levelset);
2868 FREE(mask);
2869 if (A != A0) SparseMatrix_delete(A);
2870
2871 return connected;
2872}
2873
2874
2875void SparseMatrix_decompose_to_supervariables(SparseMatrix A, int *ncluster, int **cluster, int **clusterp){
2876 /* nodes for a super variable if they share exactly the same neighbors. This is know as modules in graph theory.
2877 We work on columns only and columns with the same pattern are grouped as a super variable
2878 */
2879 int *ia = A->ia, *ja = A->ja, n = A->n, m = A->m;
2880 int *super = NULL, *nsuper = NULL, i, j, *mask = NULL, isup, *newmap, isuper;
2881
2882 super = MALLOC(sizeof(int)*((size_t)(n)));
2883 nsuper = MALLOC(sizeof(int)*((size_t)(n+1)));
2884 mask = MALLOC(sizeof(int)*((size_t)n));
2885 newmap = MALLOC(sizeof(int)*((size_t)n));
2886 nsuper++;
2887
2888 isup = 0;
2889 for (i = 0; i < n; i++) super[i] = isup;/* every node belongs to super variable 0 by default */
2890 nsuper[0] = n;
2891 for (i = 0; i < n; i++) mask[i] = -1;
2892 isup++;
2893
2894 for (i = 0; i < m; i++){
2895#ifdef DEBUG_PRINT1
2896 printf("\n");
2897 printf("doing row %d-----\n",i+1);
2898#endif
2899 for (j = ia[i]; j < ia[i+1]; j++){
2900 isuper = super[ja[j]];
2901 nsuper[isuper]--;/* those entries will move to a different super vars*/
2902 }
2903 for (j = ia[i]; j < ia[i+1]; j++){
2904 isuper = super[ja[j]];
2905 if (mask[isuper] < i){
2906 mask[isuper] = i;
2907 if (nsuper[isuper] == 0){/* all nodes in the isuper group exist in this row */
2908#ifdef DEBUG_PRINT1
2909 printf("node %d keep super node id %d\n",ja[j]+1,isuper+1);
2910#endif
2911 nsuper[isuper] = 1;
2912 newmap[isuper] = isuper;
2913 } else {
2914 newmap[isuper] = isup;
2915 nsuper[isup] = 1;
2916#ifdef DEBUG_PRINT1
2917 printf("make node %d into supernode %d\n",ja[j]+1,isup+1);
2918#endif
2919 super[ja[j]] = isup++;
2920 }
2921 } else {
2922#ifdef DEBUG_PRINT1
2923 printf("node %d join super node %d\n",ja[j]+1,newmap[isuper]+1);
2924#endif
2925 super[ja[j]] = newmap[isuper];
2926 nsuper[newmap[isuper]]++;
2927 }
2928 }
2929#ifdef DEBUG_PRINT1
2930 printf("nsuper=");
2931 for (j = 0; j < isup; j++) printf("(%d,%d),",j+1,nsuper[j]);
2932 printf("\n");
2933#endif
2934 }
2935#ifdef DEBUG_PRINT1
2936 for (i = 0; i < n; i++){
2937 printf("node %d is in supernode %d\n",i, super[i]);
2938 }
2939#endif
2940#ifdef PRINT
2941 fprintf(stderr, "n = %d, nsup = %d\n",n,isup);
2942#endif
2943 /* now accumulate super nodes */
2944 nsuper--;
2945 nsuper[0] = 0;
2946 for (i = 0; i < isup; i++) nsuper[i+1] += nsuper[i];
2947
2948 *cluster = newmap;
2949 for (i = 0; i < n; i++){
2950 isuper = super[i];
2951 (*cluster)[nsuper[isuper]++] = i;
2952 }
2953 for (i = isup; i > 0; i--) nsuper[i] = nsuper[i-1];
2954 nsuper[0] = 0;
2955 *clusterp = nsuper;
2956 *ncluster = isup;
2957
2958#ifdef PRINT
2959 for (i = 0; i < *ncluster; i++){
2960 printf("{");
2961 for (j = (*clusterp)[i]; j < (*clusterp)[i+1]; j++){
2962 printf("%d, ",(*cluster)[j]);
2963 }
2964 printf("},");
2965 }
2966 printf("\n");
2967#endif
2968
2969 FREE(mask);
2970 FREE(super);
2971
2972}
2973
2974SparseMatrix SparseMatrix_get_augmented(SparseMatrix A){
2975 /* convert matrix A to an augmente dmatrix {{0,A},{A^T,0}} */
2976 int *irn = NULL, *jcn = NULL;
2977 void *val = NULL;
2978 int nz = A->nz, type = A->type;
2979 int m = A->m, n = A->n, i, j;
2980 SparseMatrix B = NULL;
2981 if (!A) return NULL;
2982 if (nz > 0){
2983 irn = MALLOC(sizeof(int)*((size_t)nz)*2);
2984 jcn = MALLOC(sizeof(int)*((size_t)nz)*2);
2985 }
2986
2987 if (A->a){
2988 assert(A->size != 0 && nz > 0);
2989 val = MALLOC(A->size*2*((size_t)nz));
2990 MEMCPY(val, A->a, A->size*((size_t)nz));
2991 MEMCPY((void*)(((char*) val) + ((size_t)nz)*A->size), A->a, A->size*((size_t)nz));
2992 }
2993
2994 nz = 0;
2995 for (i = 0; i < m; i++){
2996 for (j = (A->ia)[i]; j < (A->ia)[i+1]; j++){
2997 irn[nz] = i;
2998 jcn[nz++] = (A->ja)[j] + m;
2999 }
3000 }
3001 for (i = 0; i < m; i++){
3002 for (j = (A->ia)[i]; j < (A->ia)[i+1]; j++){
3003 jcn[nz] = i;
3004 irn[nz++] = (A->ja)[j] + m;
3005 }
3006 }
3007
3008 B = SparseMatrix_from_coordinate_arrays(nz, m + n, m + n, irn, jcn, val, type, A->size);
3009 SparseMatrix_set_symmetric(B);
3010 SparseMatrix_set_pattern_symmetric(B);
3011 if (irn) FREE(irn);
3012 if (jcn) FREE(jcn);
3013 if (val) FREE(val);
3014 return B;
3015
3016}
3017
3018SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A, int bipartite_options){
3019 SparseMatrix B;
3020 switch (bipartite_options){
3021 case BIPARTITE_RECT:
3022 if (A->m == A->n) return A;
3023 break;
3024 case BIPARTITE_PATTERN_UNSYM:
3025 if (A->m == A->n && SparseMatrix_is_symmetric(A, TRUE)) return A;
3026 break;
3027 case BIPARTITE_UNSYM:
3028 if (A->m == A->n && SparseMatrix_is_symmetric(A, FALSE)) return A;
3029 break;
3030 case BIPARTITE_ALWAYS:
3031 break;
3032 default:
3033 assert(0);
3034 }
3035 B = SparseMatrix_get_augmented(A);
3036 SparseMatrix_delete(A);
3037 return B;
3038}
3039
3040SparseMatrix SparseMatrix_get_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices){
3041 /* get the submatrix from row/columns indices[0,...,l-1].
3042 row rindices[i] will be the new row i
3043 column cindices[i] will be the new column i.
3044 if rindices = NULL, it is assume that 1 -- nrow is needed. Same for cindices/ncol.
3045 */
3046 int nz = 0, i, j, *irn, *jcn, *ia = A->ia, *ja = A->ja, m = A->m, n = A->n;
3047 int *cmask, *rmask;
3048 void *v = NULL;
3049 SparseMatrix B = NULL;
3050 int irow = 0, icol = 0;
3051
3052 if (nrow <= 0 || ncol <= 0) return NULL;
3053
3054
3055
3056 rmask = MALLOC(sizeof(int)*((size_t)m));
3057 cmask = MALLOC(sizeof(int)*((size_t)n));
3058 for (i = 0; i < m; i++) rmask[i] = -1;
3059 for (i = 0; i < n; i++) cmask[i] = -1;
3060
3061 if (rindices){
3062 for (i = 0; i < nrow; i++) {
3063 if (rindices[i] >= 0 && rindices[i] < m){
3064 rmask[rindices[i]] = irow++;
3065 }
3066 }
3067 } else {
3068 for (i = 0; i < nrow; i++) {
3069 rmask[i] = irow++;
3070 }
3071 }
3072
3073 if (cindices){
3074 for (i = 0; i < ncol; i++) {
3075 if (cindices[i] >= 0 && cindices[i] < n){
3076 cmask[cindices[i]] = icol++;
3077 }
3078 }
3079 } else {
3080 for (i = 0; i < ncol; i++) {
3081 cmask[i] = icol++;
3082 }
3083 }
3084
3085 for (i = 0; i < m; i++){
3086 if (rmask[i] < 0) continue;
3087 for (j = ia[i]; j < ia[i+1]; j++){
3088 if (cmask[ja[j]] < 0) continue;
3089 nz++;
3090 }
3091 }
3092
3093
3094 switch (A->type){
3095 case MATRIX_TYPE_REAL:{
3096 real *a = (real*) A->a;
3097 real *val;
3098 irn = MALLOC(sizeof(int)*((size_t)nz));
3099 jcn = MALLOC(sizeof(int)*((size_t)nz));
3100 val = MALLOC(sizeof(real)*((size_t)nz));
3101
3102 nz = 0;
3103 for (i = 0; i < m; i++){
3104 if (rmask[i] < 0) continue;
3105 for (j = ia[i]; j < ia[i+1]; j++){
3106 if (cmask[ja[j]] < 0) continue;
3107 irn[nz] = rmask[i];
3108 jcn[nz] = cmask[ja[j]];
3109 val[nz++] = a[j];
3110 }
3111 }
3112 v = (void*) val;
3113 break;
3114 }
3115 case MATRIX_TYPE_COMPLEX:{
3116 real *a = (real*) A->a;
3117 real *val;
3118
3119 irn = MALLOC(sizeof(int)*((size_t)nz));
3120 jcn = MALLOC(sizeof(int)*((size_t)nz));
3121 val = MALLOC(sizeof(real)*2*((size_t)nz));
3122
3123 nz = 0;
3124 for (i = 0; i < m; i++){
3125 if (rmask[i] < 0) continue;
3126 for (j = ia[i]; j < ia[i+1]; j++){
3127 if (cmask[ja[j]] < 0) continue;
3128 irn[nz] = rmask[i];
3129 jcn[nz] = cmask[ja[j]];
3130 val[2*nz] = a[2*j];
3131 val[2*nz+1] = a[2*j+1];
3132 nz++;
3133 }
3134 }
3135 v = (void*) val;
3136 break;
3137 }
3138 case MATRIX_TYPE_INTEGER:{
3139 int *a = (int*) A->a;
3140 int *val;
3141
3142 irn = MALLOC(sizeof(int)*((size_t)nz));
3143 jcn = MALLOC(sizeof(int)*((size_t)nz));
3144 val = MALLOC(sizeof(int)*((size_t)nz));
3145
3146 nz = 0;
3147 for (i = 0; i < m; i++){
3148 if (rmask[i] < 0) continue;
3149 for (j = ia[i]; j < ia[i+1]; j++){
3150 if (cmask[ja[j]] < 0) continue;
3151 irn[nz] = rmask[i];
3152 jcn[nz] = cmask[ja[j]];
3153 val[nz] = a[j];
3154 nz++;
3155 }
3156 }
3157 v = (void*) val;
3158 break;
3159 }
3160 case MATRIX_TYPE_PATTERN:
3161 irn = MALLOC(sizeof(int)*((size_t)nz));
3162 jcn = MALLOC(sizeof(int)*((size_t)nz));
3163 nz = 0;
3164 for (i = 0; i < m; i++){
3165 if (rmask[i] < 0) continue;
3166 for (j = ia[i]; j < ia[i+1]; j++){
3167 if (cmask[ja[j]] < 0) continue;
3168 irn[nz] = rmask[i];
3169 jcn[nz++] = cmask[ja[j]];
3170 }
3171 }
3172 break;
3173 case MATRIX_TYPE_UNKNOWN:
3174 FREE(rmask);
3175 FREE(cmask);
3176 return NULL;
3177 default:
3178 FREE(rmask);
3179 FREE(cmask);
3180 return NULL;
3181 }
3182
3183 B = SparseMatrix_from_coordinate_arrays(nz, nrow, ncol, irn, jcn, v, A->type, A->size);
3184 FREE(cmask);
3185 FREE(rmask);
3186 FREE(irn);
3187 FREE(jcn);
3188 if (v) FREE(v);
3189
3190
3191 return B;
3192
3193}
3194
3195SparseMatrix SparseMatrix_exclude_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices){
3196 /* get a submatrix by excluding rows and columns */
3197 int *r, *c, nr, nc, i;
3198 SparseMatrix B;
3199
3200 if (nrow <= 0 && ncol <= 0) return A;
3201
3202 r = MALLOC(sizeof(int)*((size_t)A->m));
3203 c = MALLOC(sizeof(int)*((size_t)A->n));
3204
3205 for (i = 0; i < A->m; i++) r[i] = i;
3206 for (i = 0; i < A->n; i++) c[i] = i;
3207 for (i = 0; i < nrow; i++) {
3208 if (rindices[i] >= 0 && rindices[i] < A->m){
3209 r[rindices[i]] = -1;
3210 }
3211 }
3212 for (i = 0; i < ncol; i++) {
3213 if (cindices[i] >= 0 && cindices[i] < A->n){
3214 c[cindices[i]] = -1;
3215 }
3216 }
3217
3218 nr = nc = 0;
3219 for (i = 0; i < A->m; i++) {
3220 if (r[i] > 0) r[nr++] = r[i];
3221 }
3222 for (i = 0; i < A->n; i++) {
3223 if (c[i] > 0) c[nc++] = c[i];
3224 }
3225
3226 B = SparseMatrix_get_submatrix(A, nr, nc, r, c);
3227
3228 FREE(r);
3229 FREE(c);
3230 return B;
3231
3232}
3233
3234SparseMatrix SparseMatrix_largest_component(SparseMatrix A){
3235 SparseMatrix B;
3236 int ncomp;
3237 int *comps = NULL;
3238 int *comps_ptr = NULL;
3239 int i;
3240 int nmax, imax = 0;
3241
3242 if (!A) return NULL;
3243 A = SparseMatrix_to_square_matrix(A, BIPARTITE_RECT);
3244 SparseMatrix_weakly_connected_components(A, &ncomp, &comps, &comps_ptr);
3245 if (ncomp == 1) {
3246 B = A;
3247 } else {
3248 nmax = 0;
3249 for (i = 0; i < ncomp; i++){
3250 if (nmax < comps_ptr[i+1] - comps_ptr[i]){
3251 nmax = comps_ptr[i+1] - comps_ptr[i];
3252 imax = i;
3253 }
3254 }
3255 B = SparseMatrix_get_submatrix(A, nmax, nmax, &comps[comps_ptr[imax]], &comps[comps_ptr[imax]]);
3256 }
3257 FREE(comps);
3258 FREE(comps_ptr);
3259 return B;
3260
3261
3262}
3263
3264SparseMatrix SparseMatrix_delete_sparse_columns(SparseMatrix A, int threshold, int **new2old, int *nnew, int inplace){
3265 /* delete sparse columns of threshold or less entries in A. After than number of columns will be nnew, and
3266 the mapping from new matrix column to old matrix column is new2old.
3267 On entry, if new2old is NULL, it is allocated.
3268 */
3269 SparseMatrix B;
3270 int *ia, *ja;
3271 int *old2new;
3272 int i;
3273 old2new = MALLOC(sizeof(int)*((size_t)A->n));
3274 for (i = 0; i < A->n; i++) old2new[i] = -1;
3275
3276 *nnew = 0;
3277 B = SparseMatrix_transpose(A);
3278 ia = B->ia; ja = B->ja;
3279 for (i = 0; i < B->m; i++){
3280 if (ia[i+1] > ia[i] + threshold){
3281 (*nnew)++;
3282 }
3283 }
3284 if (!(*new2old)) *new2old = MALLOC(sizeof(int)*((size_t)(*nnew)));
3285
3286 *nnew = 0;
3287 for (i = 0; i < B->m; i++){
3288 if (ia[i+1] > ia[i] + threshold){
3289 (*new2old)[*nnew] = i;
3290 old2new[i] = *nnew;
3291 (*nnew)++;
3292 }
3293 }
3294 SparseMatrix_delete(B);
3295
3296 if (inplace){
3297 B = A;
3298 } else {
3299 B = SparseMatrix_copy(A);
3300 }
3301 ia = B->ia; ja = B->ja;
3302 for (i = 0; i < ia[B->m]; i++){
3303 assert(old2new[ja[i]] >= 0);
3304 ja[i] = old2new[ja[i]];
3305 }
3306 B->n = *nnew;
3307
3308 FREE(old2new);
3309 return B;
3310
3311
3312}
3313
3314SparseMatrix SparseMatrix_delete_empty_columns(SparseMatrix A, int **new2old, int *nnew, int inplace){
3315 return SparseMatrix_delete_sparse_columns(A, 0, new2old, nnew, inplace);
3316}
3317
3318SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A){
3319 real *a;
3320 int i;
3321
3322 if (A->a) FREE(A->a);
3323 A->a = MALLOC(sizeof(real)*((size_t)A->nz));
3324 a = (real*) (A->a);
3325 for (i = 0; i < A->nz; i++) a[i] = 1.;
3326 A->type = MATRIX_TYPE_REAL;
3327 A->size = sizeof(real);
3328 return A;
3329
3330}
3331
3332SparseMatrix SparseMatrix_complement(SparseMatrix A, int undirected){
3333 /* find the complement graph A^c, such that {i,h}\in E(A_c) iff {i,j} \notin E(A). Only structural matrix is returned. */
3334 SparseMatrix B = A;
3335 int *ia, *ja;
3336 int m = A->m, n = A->n;
3337 int *mask, nz = 0;
3338 int *irn, *jcn;
3339 int i, j;
3340
3341 if (undirected) B = SparseMatrix_symmetrize(A, TRUE);
3342 assert(m == n);
3343
3344 ia = B->ia; ja = B->ja;
3345 mask = MALLOC(sizeof(int)*((size_t)n));
3346 irn = MALLOC(sizeof(int)*(((size_t)n)*((size_t)n) - ((size_t)A->nz)));
3347 jcn = MALLOC(sizeof(int)*(((size_t)n)*((size_t)n) - ((size_t)A->nz)));
3348
3349 for (i = 0; i < n; i++){
3350 mask[i] = -1;
3351 }
3352
3353 for (i = 0; i < n; i++){
3354 for (j = ia[i]; j < ia[i+1]; j++){
3355 mask[ja[j]] = i;
3356 }
3357 for (j = 0; j < n; j++){
3358 if (mask[j] != i){
3359 irn[nz] = i;
3360 jcn[nz++] = j;
3361 }
3362 }
3363 }
3364
3365 if (B != A) SparseMatrix_delete(B);
3366 B = SparseMatrix_from_coordinate_arrays(nz, m, n, irn, jcn, NULL, MATRIX_TYPE_PATTERN, 0);
3367 FREE(irn);
3368 FREE(jcn);
3369 return B;
3370}
3371
3372int SparseMatrix_k_centers(SparseMatrix D0, int weighted, int K, int root, int **centers, int centering, real **dist0){
3373 /*
3374 Input:
3375 D: the graph. If weighted, the entry values is used.
3376 weighted: whether to treat the graph as weighted
3377 K: the number of centers
3378 root: the start node to find the k center.
3379 centering: whether the distance should be centered so that sum_k dist[n*k+i] = 0
3380 Output:
3381 centers: the list of nodes that form the k-centers. If centers = NULL on input, it will be allocated.
3382 dist: of dimension k*n, dist[k*n: (k+1)*n) gives the distance of every node to the k-th center.
3383 return: flag. if not zero, the graph is not connected, or out of memory.
3384 */
3385 SparseMatrix D = D0;
3386 int m = D->m, n = D->n;
3387 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
3388 int aggressive = FALSE;
3389 int connectedQ, end1, end2;
3390 enum {K_CENTER_DISCONNECTED = 1, K_CENTER_MEM};
3391 real *dist_min = NULL, *dist_sum = NULL, dmax, dsum;
3392 real *dist = NULL;
3393 int nlist, *list = NULL;
3394 int flag = 0, i, j, k, nlevel;
3395 int check_connected = FALSE;
3396
3397 if (!SparseMatrix_is_symmetric(D, FALSE)){
3398 D = SparseMatrix_symmetrize(D, FALSE);
3399 }
3400
3401 assert(m == n);
3402
3403 dist_min = MALLOC(sizeof(real)*n);
3404 dist_sum = MALLOC(sizeof(real)*n);
3405 for (i = 0; i < n; i++) dist_min[i] = -1;
3406 for (i = 0; i < n; i++) dist_sum[i] = 0;
3407 if (!(*centers)) *centers = MALLOC(sizeof(int)*K);
3408 if (!(*dist0)) *dist0 = MALLOC(sizeof(real)*K*n);
3409 if (!weighted){
3410 dist = MALLOC(sizeof(real)*n);
3411 SparseMatrix_pseudo_diameter_unweighted(D, root, aggressive, &end1, &end2, &connectedQ);
3412 if (check_connected && !connectedQ) {
3413 flag = K_CENTER_DISCONNECTED;
3414 goto RETURN;
3415 }
3416 root = end1;
3417 for (k = 0; k < K; k++){
3418 (*centers)[k] = root;
3419 // fprintf(stderr,"k = %d, root = %d\n",k, root+1);
3420 SparseMatrix_level_sets(D, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
3421 if (check_connected) assert(levelset_ptr[nlevel] == n);
3422 for (i = 0; i < nlevel; i++) {
3423 for (j = levelset_ptr[i]; j < levelset_ptr[i+1]; j++){
3424 (*dist0)[k*n+levelset[j]] = i;
3425 if (k == 0){
3426 dist_min[levelset[j]] = i;
3427 } else {
3428 dist_min[levelset[j]] = MIN(dist_min[levelset[j]], i);
3429 }
3430 dist_sum[levelset[j]] += i;
3431 }
3432 }
3433
3434 /* root = argmax_i min_roots dist(i, roots) */
3435 dmax = dist_min[0];
3436 dsum = dist_sum[0];
3437 root = 0;
3438 for (i = 0; i < n; i++) {
3439 if (!check_connected && dist_min[i] < 0) continue;/* if the graph is disconnected, then we can not count on every node to be in level set.
3440 Usee dist_min<0 to identify those not in level set */
3441 if (dmax < dist_min[i] || (dmax == dist_min[i] && dsum < dist_sum[i])){/* tie break with avg dist */
3442 dmax = dist_min[i];
3443 dsum = dist_sum[i];
3444 root = i;
3445 }
3446 }
3447 }
3448 } else {
3449 SparseMatrix_pseudo_diameter_weighted(D, root, aggressive, &end1, &end2, &connectedQ);
3450 if (check_connected && !connectedQ) return K_CENTER_DISCONNECTED;
3451 root = end1;
3452 list = MALLOC(sizeof(int)*n);
3453
3454 for (k = 0; k < K; k++){
3455 //fprintf(stderr,"k = %d, root = %d\n",k, root+1);
3456 (*centers)[k] = root;
3457 dist = &((*dist0)[k*n]);
3458 flag = Dijkstra(D, root, dist, &nlist, list, &dmax);
3459 if (flag){
3460 flag = K_CENTER_MEM;
3461 goto RETURN;
3462 }
3463 if (check_connected) assert(nlist == n);
3464 for (i = 0; i < n; i++){
3465 if (k == 0){
3466 dist_min[i] = dist[i];
3467 } else {
3468 dist_min[i] = MIN(dist_min[i], dist[i]);
3469 }
3470 dist_sum[i] += dist[i];
3471 }
3472
3473 /* root = argmax_i min_roots dist(i, roots) */
3474 dmax = dist_min[0];
3475 dsum = dist_sum[0];
3476 root = 0;
3477 for (i = 0; i < n; i++) {
3478 if (!check_connected && dist_min[i] < 0) continue;/* if the graph is disconnected, then we can not count on every node to be in level set.
3479 Usee dist_min<0 to identify those not in level set */
3480 if (dmax < dist_min[i] || (dmax == dist_min[i] && dsum < dist_sum[i])){/* tie break with avg dist */
3481 dmax = dist_min[i];
3482 dsum = dist_sum[i];
3483 root = i;
3484 }
3485 }
3486 }
3487 dist = NULL;
3488 }
3489
3490 if (centering){
3491 for (i = 0; i < n; i++) dist_sum[i] /= k;
3492 for (k = 0; k < K; k++){
3493 for (i = 0; i < n; i++){
3494 (*dist0)[k*n+i] -= dist_sum[i];
3495 }
3496 }
3497 }
3498
3499 RETURN:
3500 if (levelset_ptr) FREE(levelset_ptr);
3501 if (levelset) FREE(levelset);
3502 if (mask) FREE(mask);
3503
3504 if (D != D0) SparseMatrix_delete(D);
3505 if (dist) FREE(dist);
3506 if (dist_min) FREE(dist_min);
3507 if (dist_sum) FREE(dist_sum);
3508 if (list) FREE(list);
3509 return flag;
3510
3511}
3512
3513
3514
3515int SparseMatrix_k_centers_user(SparseMatrix D0, int weighted, int K, int *centers_user, int centering, real **dist0){
3516 /*
3517 Input:
3518 D: the graph. If weighted, the entry values is used.
3519 weighted: whether to treat the graph as weighted
3520 K: the number of centers
3521 root: the start node to find the k center.
3522 centering: whether the distance should be centered so that sum_k dist[n*k+i] = 0
3523 centers_user: the list of nodes that form the k-centers, GIVEN BY THE USER
3524 Output:
3525 dist: of dimension k*n, dist[k*n: (k+1)*n) gives the distance of every node to the k-th center.
3526 return: flag. if not zero, the graph is not connected, or out of memory.
3527 */
3528 SparseMatrix D = D0;
3529 int m = D->m, n = D->n;
3530 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
3531 int aggressive = FALSE;
3532 int connectedQ, end1, end2;
3533 enum {K_CENTER_DISCONNECTED = 1, K_CENTER_MEM};
3534 real *dist_min = NULL, *dist_sum = NULL, dmax;
3535 real *dist = NULL;
3536 int nlist, *list = NULL;
3537 int flag = 0, i, j, k, nlevel;
3538 int root;
3539
3540 if (!SparseMatrix_is_symmetric(D, FALSE)){
3541 D = SparseMatrix_symmetrize(D, FALSE);
3542 }
3543
3544 assert(m == n);
3545
3546 dist_min = MALLOC(sizeof(real)*n);
3547 dist_sum = MALLOC(sizeof(real)*n);
3548 for (i = 0; i < n; i++) dist_sum[i] = 0;
3549 if (!(*dist0)) *dist0 = MALLOC(sizeof(real)*K*n);
3550 if (!weighted){
3551 dist = MALLOC(sizeof(real)*n);
3552 root = centers_user[0];
3553 SparseMatrix_pseudo_diameter_unweighted(D, root, aggressive, &end1, &end2, &connectedQ);
3554 if (!connectedQ) {
3555 flag = K_CENTER_DISCONNECTED;
3556 goto RETURN;
3557 }
3558 for (k = 0; k < K; k++){
3559 root = centers_user[k];
3560 SparseMatrix_level_sets(D, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
3561 assert(levelset_ptr[nlevel] == n);
3562 for (i = 0; i < nlevel; i++) {
3563 for (j = levelset_ptr[i]; j < levelset_ptr[i+1]; j++){
3564 (*dist0)[k*n+levelset[j]] = i;
3565 if (k == 0){
3566 dist_min[levelset[j]] = i;
3567 } else {
3568 dist_min[levelset[j]] = MIN(dist_min[levelset[j]], i);
3569 }
3570 dist_sum[levelset[j]] += i;
3571 }
3572 }
3573
3574 }
3575 } else {
3576 root = centers_user[0];
3577 SparseMatrix_pseudo_diameter_weighted(D, root, aggressive, &end1, &end2, &connectedQ);
3578 if (!connectedQ) return K_CENTER_DISCONNECTED;
3579 list = MALLOC(sizeof(int)*n);
3580
3581 for (k = 0; k < K; k++){
3582 root = centers_user[k];
3583 // fprintf(stderr,"k = %d, root = %d\n",k, root+1);
3584 dist = &((*dist0)[k*n]);
3585 flag = Dijkstra(D, root, dist, &nlist, list, &dmax);
3586 if (flag){
3587 flag = K_CENTER_MEM;
3588 dist = NULL;
3589 goto RETURN;
3590 }
3591 assert(nlist == n);
3592 for (i = 0; i < n; i++){
3593 if (k == 0){
3594 dist_min[i] = dist[i];
3595 } else {
3596 dist_min[i] = MIN(dist_min[i], dist[i]);
3597 }
3598 dist_sum[i] += dist[i];
3599 }
3600
3601 }
3602 dist = NULL;
3603 }
3604
3605 if (centering){
3606 for (i = 0; i < n; i++) dist_sum[i] /= k;
3607 for (k = 0; k < K; k++){
3608 for (i = 0; i < n; i++){
3609 (*dist0)[k*n+i] -= dist_sum[i];
3610 }
3611 }
3612 }
3613
3614 RETURN:
3615 if (levelset_ptr) FREE(levelset_ptr);
3616 if (levelset) FREE(levelset);
3617 if (mask) FREE(mask);
3618
3619 if (D != D0) SparseMatrix_delete(D);
3620 if (dist) FREE(dist);
3621 if (dist_min) FREE(dist_min);
3622 if (dist_sum) FREE(dist_sum);
3623 if (list) FREE(list);
3624 return flag;
3625
3626}
3627
3628
3629
3630
3631SparseMatrix SparseMatrix_from_dense(int m, int n, real *x){
3632 /* wrap a mxn matrix into a sparse matrix. the {i,j} entry of the matrix is in x[i*n+j], 0<=i<m; 0<=j<n */
3633 int i, j, *ja;
3634 real *a;
3635 SparseMatrix A = SparseMatrix_new(m, n, m*n, MATRIX_TYPE_REAL, FORMAT_CSR);
3636
3637 A->ia[0] = 0;
3638 for (i = 1; i <= m; i++) (A->ia)[i] = (A->ia)[i-1] + n;
3639
3640 ja = A->ja;
3641 a = (real*) A->a;
3642 for (i = 0; i < m; i++){
3643 for (j = 0; j < n; j++) {
3644 ja[j] = j;
3645 a[j] = x[i*n+j];
3646 }
3647 ja += n; a += j;
3648 }
3649 A->nz = m*n;
3650 return A;
3651
3652}
3653
3654
3655int SparseMatrix_distance_matrix(SparseMatrix D0, int weighted, real **dist0){
3656 /*
3657 Input:
3658 D: the graph. If weighted, the entry values is used.
3659 weighted: whether to treat the graph as weighted
3660 Output:
3661 dist: of dimension nxn, dist[i*n+j] gives the distance of node i to j.
3662 return: flag. if not zero, the graph is not connected, or out of memory.
3663 */
3664 SparseMatrix D = D0;
3665 int m = D->m, n = D->n;
3666 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
3667 real *dist = NULL;
3668 int nlist, *list = NULL;
3669 int flag = 0, i, j, k, nlevel;
3670 real dmax;
3671
3672 if (!SparseMatrix_is_symmetric(D, FALSE)){
3673 D = SparseMatrix_symmetrize(D, FALSE);
3674 }
3675
3676 assert(m == n);
3677
3678 if (!(*dist0)) *dist0 = MALLOC(sizeof(real)*n*n);
3679 for (i = 0; i < n*n; i++) (*dist0)[i] = -1;
3680
3681 if (!weighted){
3682 for (k = 0; k < n; k++){
3683 SparseMatrix_level_sets(D, k, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
3684 assert(levelset_ptr[nlevel] == n);
3685 for (i = 0; i < nlevel; i++) {
3686 for (j = levelset_ptr[i]; j < levelset_ptr[i+1]; j++){
3687 (*dist0)[k*n+levelset[j]] = i;
3688 }
3689 }
3690 }
3691 } else {
3692 list = MALLOC(sizeof(int)*n);
3693 for (k = 0; k < n; k++){
3694 dist = &((*dist0)[k*n]);
3695 flag = Dijkstra(D, k, dist, &nlist, list, &dmax);
3696 }
3697 }
3698
3699 if (levelset_ptr) FREE(levelset_ptr);
3700 if (levelset) FREE(levelset);
3701 if (mask) FREE(mask);
3702
3703 if (D != D0) SparseMatrix_delete(D);
3704 if (list) FREE(list);
3705 return flag;
3706
3707}
3708
3709SparseMatrix SparseMatrix_distance_matrix_k_centers(int K, SparseMatrix D, int weighted){
3710 /* return a sparse matrix whichj represent the k-center and distance from every node to them.
3711 The matrix will have k*n entries
3712 */
3713 int flag;
3714 real *dist = NULL;
3715 int m = D->m, n = D->n;
3716 int root = 0;
3717 int *centers = NULL;
3718 real d;
3719 int i, j, center;
3720 SparseMatrix B, C;
3721 int centering = FALSE;
3722
3723 assert(m == n);
3724
3725 B = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
3726
3727 flag = SparseMatrix_k_centers(D, weighted, K, root, &centers, centering, &dist);
3728 assert(!flag);
3729
3730 for (i = 0; i < K; i++){
3731 center = centers[i];
3732 for (j = 0; j < n; j++){
3733 d = dist[i*n + j];
3734 B = SparseMatrix_coordinate_form_add_entries(B, 1, &center, &j, &d);
3735 B = SparseMatrix_coordinate_form_add_entries(B, 1, &j, &center, &d);
3736 }
3737 }
3738
3739 C = SparseMatrix_from_coordinate_format(B);
3740 SparseMatrix_delete(B);
3741
3742 FREE(centers);
3743 FREE(dist);
3744 return C;
3745}
3746
3747SparseMatrix SparseMatrix_distance_matrix_khops(int khops, SparseMatrix D0, int weighted){
3748 /*
3749 Input:
3750 khops: number of hops allow. If khops < 0, this will give a dense distances. Otherwise it gives a sparse matrix that represent the k-neighborhood graph
3751 D: the graph. If weighted, the entry values is used.
3752 weighted: whether to treat the graph as weighted
3753 Output:
3754 DD: of dimension nxn. DD[i,j] gives the shortest path distance, subject to the fact that the short oath must be of <= khops.
3755 return: flag. if not zero, the graph is not connected, or out of memory.
3756 */
3757 SparseMatrix D = D0, B, C;
3758 int m = D->m, n = D->n;
3759 int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
3760 real *dist = NULL;
3761 int nlist, *list = NULL;
3762 int flag = 0, i, j, k, itmp, nlevel;
3763 real dmax, dtmp;
3764
3765 if (!SparseMatrix_is_symmetric(D, FALSE)){
3766 D = SparseMatrix_symmetrize(D, FALSE);
3767 }
3768
3769 assert(m == n);
3770
3771 B = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
3772
3773 if (!weighted){
3774 for (k = 0; k < n; k++){
3775 SparseMatrix_level_sets_khops(khops, D, k, &nlevel, &levelset_ptr, &levelset, &mask, TRUE);
3776 for (i = 0; i < nlevel; i++) {
3777 for (j = levelset_ptr[i]; j < levelset_ptr[i+1]; j++){
3778 itmp = levelset[j]; dtmp = i;
3779 if (k != itmp) B = SparseMatrix_coordinate_form_add_entries(B, 1, &k, &itmp, &dtmp);
3780 }
3781 }
3782 }
3783 } else {
3784 list = MALLOC(sizeof(int)*n);
3785 dist = MALLOC(sizeof(real)*n);
3786 /*
3787 Dijkstra_khops(khops, D, 60, dist, &nlist, list, &dmax);
3788 for (j = 0; j < nlist; j++){
3789 fprintf(stderr,"{%d,%d}=%f,",60,list[j],dist[list[j]]);
3790 }
3791 fprintf(stderr,"\n");
3792 Dijkstra_khops(khops, D, 94, dist, &nlist, list, &dmax);
3793 for (j = 0; j < nlist; j++){
3794 fprintf(stderr,"{%d,%d}=%f,",94,list[j],dist[list[j]]);
3795 }
3796 fprintf(stderr,"\n");
3797 exit(1);
3798
3799 */
3800
3801 for (k = 0; k < n; k++){
3802 SparseMatrix_level_sets_khops(khops, D, k, &nlevel, &levelset_ptr, &levelset, &mask, FALSE);
3803 assert(nlevel-1 <= khops);/* the first level is the root */
3804 flag = Dijkstra_masked(D, k, dist, &nlist, list, &dmax, mask);
3805 assert(!flag);
3806 for (i = 0; i < nlevel; i++) {
3807 for (j = levelset_ptr[i]; j < levelset_ptr[i+1]; j++){
3808 assert(mask[levelset[j]] == i+1);
3809 mask[levelset[j]] = -1;
3810 }
3811 }
3812 for (j = 0; j < nlist; j++){
3813 itmp = list[j]; dtmp = dist[itmp];
3814 if (k != itmp) B = SparseMatrix_coordinate_form_add_entries(B, 1, &k, &itmp, &dtmp);
3815 }
3816 }
3817 }
3818
3819 C = SparseMatrix_from_coordinate_format(B);
3820 SparseMatrix_delete(B);
3821
3822 if (levelset_ptr) FREE(levelset_ptr);
3823 if (levelset) FREE(levelset);
3824 if (mask) FREE(mask);
3825 if (dist) FREE(dist);
3826
3827 if (D != D0) SparseMatrix_delete(D);
3828 if (list) FREE(list);
3829 /* I can not find a reliable way to make the matrix symmetric. Right now I use a mask array to
3830 limit consider of only nodes with in k hops, but even this is not symmetric. e.g.,
3831 . 10 10 10 10
3832 .A---B---C----D----E
3833 . 2 | | 2
3834 . G----F
3835 . 2
3836 If we set hops = 4, and from A, it can not see F (which is 5 hops), hence distance(A,E) =40
3837 but from E, it can see all nodes (all within 4 hops), so distance(E, A)=36.
3838 .
3839 may be there is a better way to ensure symmetric, but for now we just symmetrize it
3840 */
3841 D = SparseMatrix_symmetrize(C, FALSE);
3842 SparseMatrix_delete(C);
3843 return D;
3844
3845}
3846
3847#if PQ
3848void SparseMatrix_kcore_decomposition(SparseMatrix A, int *coreness_max0, int **coreness_ptr0, int **coreness_list0){
3849 /* give an undirected graph A, find the k-coreness of each vertex
3850 A: a graph. Will be made undirected and self loop removed
3851 coreness_max: max core number.
3852 coreness_ptr: array of size (coreness_max + 2), element [corness_ptr[i], corness_ptr[i+1])
3853 . of array coreness_list gives the vertices with core i, i <= coreness_max
3854 coreness_list: array of size n = A->m
3855 */
3856 SparseMatrix B;
3857 int i, j, *ia, *ja, n = A->m, nz, istatus, neighb;
3858 PriorityQueue pq = NULL;
3859 int gain, deg, k, deg_max = 0, deg_old;
3860 int *coreness_ptr, *coreness_list, coreness_now;
3861 int *mask;
3862
3863 assert(A->m == A->n);
3864 B = SparseMatrix_symmetrize(A, FALSE);
3865 B = SparseMatrix_remove_diagonal(B);
3866 ia = B->ia;
3867 ja = B->ja;
3868
3869 mask = MALLOC(sizeof(int)*n);
3870 for (i = 0; i < n; i++) mask[i] = -1;
3871
3872 pq = PriorityQueue_new(n, n-1);
3873 for (i = 0; i < n; i++){
3874 deg = ia[i+1] - ia[i];
3875 deg_max = MAX(deg_max, deg);
3876 gain = n - 1 - deg;
3877 pq = PriorityQueue_push(pq, i, gain);
3878 //fprintf(stderr,"insert %d with gain %d\n",i, gain);
3879 }
3880
3881
3882 coreness_ptr = MALLOC(sizeof(int)*(deg_max+2));
3883 coreness_list = MALLOC(sizeof(int)*n);
3884 deg_old = 0;
3885 coreness_ptr[deg_old] = 0;
3886 coreness_now = 0;
3887
3888 nz = 0;
3889 while (PriorityQueue_pop(pq, &k, &gain)){
3890 deg = (n-1) - gain;
3891 if (deg > deg_old) {
3892 //fprintf(stderr,"deg = %d, cptr[%d--%d]=%d\n",deg, deg_old + 1, deg, nz);
3893 for (j = deg_old + 1; j <= deg; j++) coreness_ptr[j] = nz;
3894 coreness_now = deg;
3895 deg_old = deg;
3896 }
3897 coreness_list[nz++] = k;
3898 mask[k] = coreness_now;
3899 //fprintf(stderr,"=== \nremove node %d with gain %d, mask with %d, nelement=%d\n",k, gain, coreness_now, pq->count);
3900 for (j = ia[k]; j < ia[k+1]; j++){
3901 neighb = ja[j];
3902 if (mask[neighb] < 0){
3903 gain = PriorityQueue_get_gain(pq, neighb);
3904 //fprintf(stderr,"update node %d with gain %d, nelement=%d\n",neighb, gain+1, pq->count);
3905 istatus = PriorityQueue_remove(pq, neighb);
3906 assert(istatus != 0);
3907 pq = PriorityQueue_push(pq, neighb, gain + 1);
3908 }
3909 }
3910 }
3911 coreness_ptr[coreness_now + 1] = nz;
3912
3913 *coreness_max0 = coreness_now;
3914 *coreness_ptr0 = coreness_ptr;
3915 *coreness_list0 = coreness_list;
3916
3917 if (Verbose){
3918 for (i = 0; i <= coreness_now; i++){
3919 if (coreness_ptr[i+1] - coreness_ptr[i] > 0){
3920 fprintf(stderr,"num_in_core[%d] = %d: ",i, coreness_ptr[i+1] - coreness_ptr[i]);
3921#if 0
3922 for (j = coreness_ptr[i]; j < coreness_ptr[i+1]; j++){
3923 fprintf(stderr,"%d,",coreness_list[j]);
3924 }
3925#endif
3926 fprintf(stderr,"\n");
3927 }
3928 }
3929 }
3930 if (Verbose)
3931
3932
3933 if (B != A) SparseMatrix_delete(B);
3934 FREE(mask);
3935}
3936
3937void SparseMatrix_kcoreness(SparseMatrix A, int **coreness){
3938
3939 int coreness_max, *coreness_ptr = NULL, *coreness_list = NULL, i, j;
3940
3941 if (!(*coreness)) coreness = MALLOC(sizeof(int)*A->m);
3942
3943 SparseMatrix_kcore_decomposition(A, &coreness_max, &coreness_ptr, &coreness_list);
3944
3945 for (i = 0; i <= coreness_max; i++){
3946 for (j = coreness_ptr[i]; j < coreness_ptr[i+1]; j++){
3947 (*coreness)[coreness_list[j]] = i;
3948 }
3949 }
3950
3951 assert(coreness_ptr[coreness_ptr[coreness_max+1]] = A->m);
3952
3953}
3954
3955
3956
3957
3958void SparseMatrix_khair_decomposition(SparseMatrix A, int *hairness_max0, int **hairness_ptr0, int **hairness_list0){
3959 /* define k-hair as the largest subgraph of the graph such that the degree of each node is <=k.
3960 Give an undirected graph A, find the k-hairness of each vertex
3961 A: a graph. Will be made undirected and self loop removed
3962 hairness_max: max hair number.
3963 hairness_ptr: array of size (hairness_max + 2), element [corness_ptr[i], corness_ptr[i+1])
3964 . of array hairness_list gives the vertices with hair i, i <= hairness_max
3965 hairness_list: array of size n = A->m
3966 */
3967 SparseMatrix B;
3968 int i, j, jj, *ia, *ja, n = A->m, nz, istatus, neighb;
3969 PriorityQueue pq = NULL;
3970 int gain, deg = 0, k, deg_max = 0, deg_old;
3971 int *hairness_ptr, *hairness_list, l;
3972 int *mask;
3973
3974 assert(A->m == A->n);
3975 B = SparseMatrix_symmetrize(A, FALSE);
3976 B = SparseMatrix_remove_diagonal(B);
3977 ia = B->ia;
3978 ja = B->ja;
3979
3980 mask = MALLOC(sizeof(int)*n);
3981 for (i = 0; i < n; i++) mask[i] = -1;
3982
3983 pq = PriorityQueue_new(n, n-1);
3984 for (i = 0; i < n; i++){
3985 deg = ia[i+1] - ia[i];
3986 deg_max = MAX(deg_max, deg);
3987 gain = deg;
3988 pq = PriorityQueue_push(pq, i, gain);
3989 }
3990
3991
3992 hairness_ptr = MALLOC(sizeof(int)*(deg_max+2));
3993 hairness_list = MALLOC(sizeof(int)*n);
3994 deg_old = deg_max;
3995 hairness_ptr[deg_old + 1] = n;
3996
3997 nz = n - 1;
3998 while (PriorityQueue_pop(pq, &k, &gain)){
3999 deg = gain;
4000 mask[k] = deg;
4001
4002 if (deg < deg_old) {
4003 //fprintf(stderr,"cptr[%d--%d]=%d\n",deg, deg_old + 1, nz);
4004 for (j = deg_old; j >= deg; j--) hairness_ptr[j] = nz + 1;
4005
4006 for (jj = hairness_ptr[deg_old]; jj < hairness_ptr [deg_old+1]; jj++){
4007 l = hairness_list[jj];
4008 //fprintf(stderr,"=== \nremove node hairness_list[%d]= %d, mask with %d, nelement=%d\n",jj, l, deg_old, pq->count);
4009 for (j = ia[l]; j < ia[l+1]; j++){
4010 neighb = ja[j];
4011 if (neighb == k) deg--;/* k was masked. But we do need to update ts degree */
4012 if (mask[neighb] < 0){
4013 gain = PriorityQueue_get_gain(pq, neighb);
4014 //fprintf(stderr,"update node %d with deg %d, nelement=%d\n",neighb, gain-1, pq->count);
4015 istatus = PriorityQueue_remove(pq, neighb);
4016 assert(istatus != 0);
4017 pq = PriorityQueue_push(pq, neighb, gain - 1);
4018 }
4019 }
4020 }
4021 mask[k] = 0;/* because a bunch of nodes are removed, k may not be the best node! Unmask */
4022 pq = PriorityQueue_push(pq, k, deg);
4023 PriorityQueue_pop(pq, &k, &gain);
4024 deg = gain;
4025 mask[k] = deg;
4026 deg_old = deg;
4027 }
4028 //fprintf(stderr,"-------- node with highes deg is %d, deg = %d\n",k,deg);
4029 //fprintf(stderr,"hairness_lisrt[%d]=%d, mask[%d] = %d\n",nz,k, k, deg);
4030 assert(deg == deg_old);
4031 hairness_list[nz--] = k;
4032 }
4033 hairness_ptr[deg] = nz + 1;
4034 assert(nz + 1 == 0);
4035 for (i = 0; i < deg; i++) hairness_ptr[i] = 0;
4036
4037 *hairness_max0 = deg_max;
4038 *hairness_ptr0 = hairness_ptr;
4039 *hairness_list0 = hairness_list;
4040
4041 if (Verbose){
4042 for (i = 0; i <= deg_max; i++){
4043 if (hairness_ptr[i+1] - hairness_ptr[i] > 0){
4044 fprintf(stderr,"num_in_hair[%d] = %d: ",i, hairness_ptr[i+1] - hairness_ptr[i]);
4045#if 0
4046 for (j = hairness_ptr[i]; j < hairness_ptr[i+1]; j++){
4047 fprintf(stderr,"%d,",hairness_list[j]);
4048 }
4049#endif
4050 fprintf(stderr,"\n");
4051 }
4052 }
4053 }
4054 if (Verbose)
4055
4056
4057 if (B != A) SparseMatrix_delete(B);
4058 FREE(mask);
4059}
4060
4061
4062void SparseMatrix_khairness(SparseMatrix A, int **hairness){
4063
4064 int hairness_max, *hairness_ptr = NULL, *hairness_list = NULL, i, j;
4065
4066 if (!(*hairness)) hairness = MALLOC(sizeof(int)*A->m);
4067
4068 SparseMatrix_khair_decomposition(A, &hairness_max, &hairness_ptr, &hairness_list);
4069
4070 for (i = 0; i <= hairness_max; i++){
4071 for (j = hairness_ptr[i]; j < hairness_ptr[i+1]; j++){
4072 (*hairness)[hairness_list[j]] = i;
4073 }
4074 }
4075
4076 assert(hairness_ptr[hairness_ptr[hairness_max+1]] = A->m);
4077
4078}
4079#endif
4080
4081void SparseMatrix_page_rank(SparseMatrix A, real teleport_probablity, int weighted, real epsilon, real **page_rank){
4082 /* A(i,j)/Sum_k A(i,k) gives the probablity of the random surfer walking from i to j
4083 A: n x n square matrix
4084 weighted: whether to use the wedge weights (matrix entries)
4085 page_rank: array of length n. If *page_rank was null on entry, will be assigned.
4086
4087 */
4088 int n = A->n;
4089 int i, j;
4090 int *ia = A->ia, *ja = A->ja;
4091 real *x, *y, *diag, res;
4092 real *a = NULL;
4093 int iter = 0;
4094
4095 assert(A->m == n);
4096 assert(teleport_probablity >= 0);
4097
4098 if (weighted){
4099 switch (A->type){
4100 case MATRIX_TYPE_REAL:
4101 a = (real*) A->a;
4102 break;
4103 case MATRIX_TYPE_COMPLEX:/* take real part */
4104 a = MALLOC(sizeof(real)*n);
4105 for (i = 0; i < n; i++) a[i] = ((real*) A->a)[2*i];
4106 break;
4107 case MATRIX_TYPE_INTEGER:
4108 a = MALLOC(sizeof(real)*n);
4109 for (i = 0; i < n; i++) a[i] = ((int*) A->a)[i];
4110 break;
4111 case MATRIX_TYPE_PATTERN:
4112 case MATRIX_TYPE_UNKNOWN:
4113 default:
4114 weighted = FALSE;
4115 break;
4116 }
4117 }
4118
4119
4120 if (!(*page_rank)) *page_rank = MALLOC(sizeof(real)*n);
4121 x = *page_rank;
4122
4123 diag = MALLOC(sizeof(real)*n);
4124 for (i = 0; i < n; i++) diag[i] = 0;
4125 y = MALLOC(sizeof(real)*n);
4126
4127 for (i = 0; i < n; i++) x[i] = 1./n;
4128
4129 /* find the column sum */
4130 if (weighted){
4131 for (i = 0; i < n; i++){
4132 for (j = ia[i]; j < ia[i+1]; j++){
4133 if (ja[j] == i) continue;
4134 diag[i] += ABS(a[j]);
4135 }
4136 }
4137 } else {
4138 for (i = 0; i < n; i++){
4139 for (j = ia[i]; j < ia[i+1]; j++){
4140 if (ja[j] == i) continue;
4141 diag[i]++;
4142 }
4143 }
4144 }
4145 for (i = 0; i < n; i++) diag[i] = 1./MAX(diag[i], MACHINEACC);
4146
4147 /* iterate */
4148 do {
4149 iter++;
4150 for (i = 0; i < n; i++) y[i] = 0;
4151 if (weighted){
4152 for (i = 0; i < n; i++){
4153 for (j = ia[i]; j < ia[i+1]; j++){
4154 if (ja[j] == i) continue;
4155 y[ja[j]] += a[j]*x[i]*diag[i];
4156 }
4157 }
4158 } else {
4159 for (i = 0; i < n; i++){
4160 for (j = ia[i]; j < ia[i+1]; j++){
4161 if (ja[j] == i) continue;
4162 y[ja[j]] += x[i]*diag[i];
4163 }
4164 }
4165 }
4166 for (i = 0; i < n; i++){
4167 y[i] = (1-teleport_probablity)*y[i] + teleport_probablity/n;
4168 }
4169
4170 /*
4171 fprintf(stderr,"\n============\nx=");
4172 for (i = 0; i < n; i++) fprintf(stderr,"%f,",x[i]);
4173 fprintf(stderr,"\nx=");
4174 for (i = 0; i < n; i++) fprintf(stderr,"%f,",y[i]);
4175 fprintf(stderr,"\n");
4176 */
4177
4178 res = 0;
4179 for (i = 0; i < n; i++) res += ABS(x[i] - y[i]);
4180 if (Verbose) fprintf(stderr,"page rank iter -- %d, res = %f\n",iter, res);
4181 MEMCPY(x, y, sizeof(real)*n);
4182 } while (res > epsilon);
4183
4184 FREE(y);
4185 FREE(diag);
4186 if (a && a != A->a) FREE(a);
4187}
4188