Graphviz  2.41.20170921.2350
SparseMatrix.c
Go to the documentation of this file.
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 
28 static 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;
35  size = 2*sizeof(real);
36  break;
38  size = sizeof(int);
39  break;
41  size = 0;
42  break;
44  size = 0;
45  break;
46  default:
47  size = 0;
48  break;
49  }
50 
51  return size;
52 }
53 
55  SparseMatrix B;
60  return A;
61 }
63  /* make it strictly low diag only, and set flag to undirected */
64  SparseMatrix B;
67  return SparseMatrix_remove_upper(B);
68 }
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:
137  return NULL;
138  default:
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 
151 SparseMatrix 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);
160  return A;
161 }
162 
164  SparseMatrix B;
165  if (SparseMatrix_is_symmetric(A, pattern_symmetric_only)) {
166  B = SparseMatrix_copy(A);
168  }
169  B = SparseMatrix_transpose(A);
170  if (!B) return NULL;
171  A = SparseMatrix_add(A, B);
176 }
177 
178 int 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){
285  } else {
288  }
289  RETURN:
290  FREE(mask);
291 
293  return res;
294 }
295 
296 static 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;
324  return A;
325 }
326 
327 static 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 
351 static 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 
384 SparseMatrix 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 }
397 SparseMatrix 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 
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 }
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 
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 
537  switch (A->format){
538  case FORMAT_CSR:
540  break;
541  case FORMAT_CSC:
542  assert(0);/* not implemented yet...
543  SparseMatrix_print_csc(c, A);*/
544  break;
545  case FORMAT_COORD:
547  break;
548  default:
549  assert(0);
550  }
551 }
552 
553 
554 
555 
556 
557 static 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 
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 
646 void 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  }
656  fclose(f);
657 
658 }
659 
660 
661 
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 
708  SparseMatrix A = NULL;
709  FILE *f;
710  f = fopen(name, "rb");
711 
713  return A;
714 }
715 
716 static 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 
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 
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 }
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 
828 static 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 
957 SparseMatrix 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 
962 SparseMatrix 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 
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 
1100 static 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 
1116 static 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 
1156 static 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 
1192 void 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 
1217 void 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);
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 
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 }
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 
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:
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 
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:
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 */
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 
1968 SparseMatrix 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 
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 
2073 SparseMatrix 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 
2153  return A;
2154 }
2155 
2156 
2157 
2158 
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 
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 
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);
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 
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 
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 
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  }
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  }
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 
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 
2378 SparseMatrix 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 
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 
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 
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 
2495 void 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 
2545 void 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 }
2552 void 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 
2558 void 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 
2590  real dist;/* distance to root */
2591  int id;/*node id */
2592 };
2593 typedef struct nodedata_struct* nodedata;
2594 
2595 static 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 
2609 static 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 
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 
2720 static 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 
2724 static 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 
2732 real 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 
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 
2796 real 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 
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 
2850  int end1, end2, connectedQ;
2851  return SparseMatrix_pseudo_diameter_unweighted(A, 0, FALSE, &end1, &end2, &connectedQ);
2852 }
2853 
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 
2875 void 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 
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);
3011  if (irn) FREE(irn);
3012  if (jcn) FREE(jcn);
3013  if (val) FREE(val);
3014  return B;
3015 
3016 }
3017 
3019  SparseMatrix B;
3020  switch (bipartite_options){
3021  case BIPARTITE_RECT:
3022  if (A->m == A->n) return A;
3023  break;
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  }
3037  return B;
3038 }
3039 
3040 SparseMatrix 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 
3195 SparseMatrix 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 
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;
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 
3264 SparseMatrix 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  }
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 
3314 SparseMatrix 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 
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 
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 
3372 int 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)){
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 
3515 int 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)){
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 
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;
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 
3655 int 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)){
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 
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 
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 
3741 
3742  FREE(centers);
3743  FREE(dist);
3744  return C;
3745 }
3746 
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)){
3767  }
3768 
3769  assert(m == n);
3770 
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 
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  */
3843  return D;
3844 
3845 }
3846 
3847 #if PQ
3848 void 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);
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 
3937 void 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 
3958 void 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);
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 
4062 void 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 
4081 void 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 }
SparseMatrix SparseMatrix_get_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices)
void SparseMatrix_export_binary_fp(FILE *f, SparseMatrix A)
Definition: SparseMatrix.c:626
#define MAX(a, b)
Definition: agerror.c:17
void SparseMatrix_khairness(SparseMatrix A, int **hairness)
SparseMatrix SparseMatrix_multiply3(SparseMatrix A, SparseMatrix B, SparseMatrix C)
#define ABS(a)
Definition: arith.h:45
void SparseMatrix_level_sets_internal(int khops, SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask)
SparseMatrix SparseMatrix_import_binary_fp(FILE *f)
Definition: SparseMatrix.c:662
double xmax
Definition: geometry.c:20
SparseMatrix SparseMatrix_multiply_by_scaler(SparseMatrix A, real s)
SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum)
int BinaryHeap_reset(BinaryHeap h, int id, void *item)
Definition: BinaryHeap.c:208
SparseMatrix SparseMatrix_symmetrize_nodiag(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:163
real SparseMatrix_pseudo_diameter_only(SparseMatrix A)
#define MIN(a, b)
Definition: arith.h:35
SparseMatrix SparseMatrix_normalize_to_rowsum1(SparseMatrix A)
void SparseMatrix_page_rank(SparseMatrix A, real teleport_probablity, int weighted, real epsilon, real **page_rank)
SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A, int bipartite_options)
SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A)
SparseMatrix SparseMatrix_copy(SparseMatrix A)
SparseMatrix SparseMatrix_scaled_by_vector(SparseMatrix A, real *v, int apply_to_row)
#define C
Definition: pack.c:29
SparseMatrix SparseMatrix_remove_upper(SparseMatrix A)
#define FREE
Definition: PriorityQueue.c:23
SparseMatrix SparseMatrix_import_binary(char *name)
Definition: SparseMatrix.c:707
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
Definition: SparseMatrix.c:384
PriorityQueue PriorityQueue_push(PriorityQueue q, int i, int gain)
Definition: PriorityQueue.c:65
#define assert(x)
Definition: cghdr.h:47
SparseMatrix SparseMatrix_exclude_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices)
BinaryHeap BinaryHeap_new(int(*cmp)(void *item1, void *item2))
Definition: BinaryHeap.c:16
int PriorityQueue_remove(PriorityQueue q, int i)
int SparseMatrix_connectedQ(SparseMatrix A0)
double xmin
Definition: geometry.c:20
#define SparseMatrix_known_symmetric(A)
Definition: SparseMatrix.h:176
SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A)
void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag)
Definition: SparseMatrix.c:646
#define clear_flag(a, flag)
Definition: general.h:38
double ymax
Definition: geometry.c:20
void SparseMatrix_multiply_dense(SparseMatrix A, int ATransposed, real *v, int vTransposed, real **res, int res_transposed, int dim)
double ymin
Definition: geometry.c:20
void SparseMatrix_decompose_to_supervariables(SparseMatrix A, int *ncluster, int **cluster, int **clusterp)
SparseMatrix 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)
Definition: SparseMatrix.c:962
int SparseMatrix_has_diagonal(SparseMatrix A)
int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only)
Definition: SparseMatrix.c:178
int SparseMatrix_k_centers_user(SparseMatrix D0, int weighted, int K, int *centers_user, int centering, real **dist0)
void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int transposed)
int BinaryHeap_insert(BinaryHeap h, void *item)
Definition: BinaryHeap.c:137
SparseMatrix SparseMatrix_largest_component(SparseMatrix A)
void SparseMatrix_kcore_decomposition(SparseMatrix A, int *coreness_max0, int **coreness_ptr0, int **coreness_list0)
SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix A)
int
Definition: grammar.c:1264
#define max(x, y)
Definition: stress.c:794
SparseMatrix SparseMatrix_delete_sparse_columns(SparseMatrix A, int threshold, int **new2old, int *nnew, int inplace)
SparseMatrix SparseMatrix_general_new(int m, int n, int nz, int type, size_t sz, int format)
Definition: SparseMatrix.c:397
Definition: circular.h:30
void SparseMatrix_print(char *c, SparseMatrix A)
Definition: SparseMatrix.c:536
SparseMatrix SparseMatrix_get_augmented(SparseMatrix A)
void SparseMatrix_kcoreness(SparseMatrix A, int **coreness)
#define SparseMatrix_set_pattern_symmetric(A)
Definition: SparseMatrix.h:163
SparseMatrix SparseMatrix_from_coordinate_format_not_compacted(SparseMatrix A, int what_to_sum)
Definition: SparseMatrix.c:812
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A)
void SparseMatrix_export(FILE *f, SparseMatrix A)
Definition: SparseMatrix.c:778
SparseMatrix SparseMatrix_apply_fun_general(SparseMatrix A, void(*fun)(int i, int j, int n, double *x))
int SparseMatrix_distance_matrix(SparseMatrix D0, int weighted, real **dist0)
void SparseMatrix_print_csr(char *c, SparseMatrix A)
Definition: SparseMatrix.c:420
PriorityQueue PriorityQueue_new(int n, int ngain)
Definition: PriorityQueue.c:27
#define MALLOC
Definition: PriorityQueue.c:21
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:151
void SparseMatrix_level_sets(SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask)
void SparseMatrix_level_sets_khops(int khops, SparseMatrix A, int root, int *nlevel, int **levelset_ptr, int **levelset, int **mask, int reinitialize_mask)
void * BinaryHeap_extract_min(BinaryHeap h)
Definition: BinaryHeap.c:169
void SparseMatrix_delete(SparseMatrix A)
Definition: SparseMatrix.c:411
Definition: grammar.c:79
#define REALLOC
Definition: PriorityQueue.c:22
SparseMatrix SparseMatrix_from_dense(int m, int n, real *x)
real SparseMatrix_pseudo_diameter_unweighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ)
#define MEMCPY
Definition: PriorityQueue.c:24
if(aagss+aagstacksize-1<=aagssp)
Definition: grammar.c:1332
SparseMatrix SparseMatrix_sort(SparseMatrix A)
Definition: SparseMatrix.c:54
SparseMatrix SparseMatrix_delete_empty_columns(SparseMatrix A, int **new2old, int *nnew, int inplace)
struct nodedata_struct * nodedata
#define NULL
Definition: logic.h:39
real SparseMatrix_pseudo_diameter_weighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ)
int SparseMatrix_k_centers(SparseMatrix D0, int weighted, int K, int root, int **centers, int centering, real **dist0)
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
Definition: SparseMatrix.c:957
SparseMatrix SparseMatrix_crop(SparseMatrix A, real epsilon)
SparseMatrix SparseMatrix_distance_matrix_k_centers(int K, SparseMatrix D, int weighted)
EXTERN unsigned char Verbose
Definition: globals.h:64
for(;;)
Definition: grammar.c:1846
void * BinaryHeap_get_item(BinaryHeap h, int id)
Definition: BinaryHeap.c:225
SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B)
void SparseMatrix_khair_decomposition(SparseMatrix A, int *hairness_max0, int **hairness_ptr0, int **hairness_list0)
SparseMatrix SparseMatrix_make_undirected(SparseMatrix A)
Definition: SparseMatrix.c:62
int PriorityQueue_pop(PriorityQueue q, int *i, int *gain)
SparseMatrix SparseMatrix_distance_matrix_khops(int khops, SparseMatrix D0, int weighted)
#define SparseMatrix_known_strucural_symmetric(A)
Definition: SparseMatrix.h:177
SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A)
SparseMatrix SparseMatrix_to_complex(SparseMatrix A)
void SparseMatrix_weakly_connected_components(SparseMatrix A0, int *ncomp, int **comps, int **comps_ptr)
int PriorityQueue_get_gain(PriorityQueue q, int i)
double dist(Site *s, Site *t)
Definition: site.c:41
SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double(*fun)(double x))
void BinaryHeap_delete(BinaryHeap h, void(*del)(void *item))
Definition: BinaryHeap.c:34
SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B)
Definition: SparseMatrix.c:966
SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, int nentries, int *irn, int *jcn, void *val)
#define MACHINEACC
Definition: general.h:120
#define SparseMatrix_set_symmetric(A)
Definition: SparseMatrix.h:162
void SparseMatrix_print_coord(char *c, SparseMatrix A)
Definition: SparseMatrix.c:481
#define FALSE
Definition: cgraph.h:35
#define SYMMETRY_EPSILON
Definition: SparseMatrix.h:19
SparseMatrix SparseMatrix_from_coordinate_format(SparseMatrix A)
Definition: SparseMatrix.c:797
SparseMatrix SparseMatrix_transpose(SparseMatrix A)
Definition: SparseMatrix.c:69
#define SparseMatrix_set_undirected(A)
Definition: SparseMatrix.h:161
SparseMatrix SparseMatrix_complement(SparseMatrix A, int undirected)
#define TRUE
Definition: cgraph.h:38
#define real
Definition: general.h:34