Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
matrix_io.c
Go to the documentation of this file.
1 
9 #include "my_util.h"
10 
11 
12 void *my_malloc (int sz)
13 {
14  void *ptr = malloc(sz);
15 
16  if (!ptr){
17  printf ("Error: can't allocate %d bytes\n", sz);
18  exit (0);
19  }
20 #ifdef DEBUG
21  printf("allocating %d bytes to %p\n", sz, ptr);
22 #endif
23 
24  return ptr;
25 }
26 
27 /* generates random integer in [low, high) */
28 int random_integer (int low, int high)
29 {
30  int r;
31  r = ((high-low)*drand48()) + low;
32  if (r==high) r--;
33 
34  return r;
35 }
36 
37 /* generates random double in [low, high) */
38 double random_double (double low, double high)
39 {
40  return ((high-low)*drand48()) + low;
41 }
42 /*
43 MT-LEVEL
44  Safe
45 
46 DESCRIPTION
47  This family of functions generates pseudo-random numbers
48  using the well-known linear congruential algorithm and 48-
49  bit integer arithmetic.
50 
51  Functions drand48() and erand48() return non-negative
52  double-precision floating-point values uniformly distributed
53  over the interval [0.0, 1.0).
54 */
55 
56 
57 int str_to_mem_unit(char *str){
58  int v = 0;
59 
60  while (*str >= '0' && *str <='9'){
61  v = v*10 + *str - '0'; str++;
62  }
63 
64  if (*str == 'm' || *str == 'M'){
65  v *= 1024*1024;
66  }
67  else if (*str == 'k' || *str == 'K'){
68  v *= 1024;
69  }
70 
71  return v;
72 }
73 
74 
75 #include "mmio.h"
76 
77 
79 {
80  if (!mm_is_matrix(matcode)) return 0;
81  if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
82  if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
83  if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
84  mm_is_skew(matcode))) return 0;
85  return 1;
86 }
87 
88 int mm_read_banner(FILE *f, MM_typecode *matcode)
89 {
90  char line[MM_MAX_LINE_LENGTH];
91  char banner[MM_MAX_TOKEN_LENGTH];
92  char mtx[MM_MAX_TOKEN_LENGTH];
93  char crd[MM_MAX_TOKEN_LENGTH];
94  char data_type[MM_MAX_TOKEN_LENGTH];
95  char storage_scheme[MM_MAX_TOKEN_LENGTH];
96  char *p;
97 
98 
99  mm_clear_typecode(matcode);
100 
101  if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
102  return MM_PREMATURE_EOF;
103 
104  if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
105  storage_scheme) != 5)
106  return MM_PREMATURE_EOF;
107 
108  for (p=mtx; *p!='\0'; *p=tolower(*p),p++); /* convert to lower case */
109  for (p=crd; *p!='\0'; *p=tolower(*p),p++);
110  for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
111  for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
112 
113  /* check for banner */
114  if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
115  return MM_NO_HEADER;
116 
117  /* first field should be "mtx" */
118  if (strcmp(mtx, MM_MTX_STR) != 0)
119  return MM_UNSUPPORTED_TYPE;
120  mm_set_matrix(matcode);
121 
122 
123  /* second field describes whether this is a sparse matrix (in coordinate
124  storgae) or a dense array */
125 
126 
127  if (strcmp(crd, MM_SPARSEROW_STR) == 0)
128  mm_set_sparserow(matcode);
129  else
130  if (strcmp(crd, MM_COORDINATE_STR) == 0)
131  mm_set_coordinate(matcode);
132  else
133  if (strcmp(crd, MM_DENSE_STR) == 0)
134  mm_set_dense(matcode);
135  else
136  return MM_UNSUPPORTED_TYPE;
137 
138 
139  /* third field */
140 
141  if (strcmp(data_type, MM_REAL_STR) == 0)
142  mm_set_real(matcode);
143  else
144  if (strcmp(data_type, MM_COMPLEX_STR) == 0)
145  mm_set_complex(matcode);
146  else
147  if (strcmp(data_type, MM_PATTERN_STR) == 0)
148  mm_set_pattern(matcode);
149  else
150  if (strcmp(data_type, MM_INT_STR) == 0)
151  mm_set_integer(matcode);
152  else
153  return MM_UNSUPPORTED_TYPE;
154 
155 
156  /* fourth field */
157 
158  if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
159  mm_set_general(matcode);
160  else
161  if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
162  mm_set_symmetric(matcode);
163  else
164  if (strcmp(storage_scheme, MM_HERM_STR) == 0)
165  mm_set_hermitian(matcode);
166  else
167  if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
168  mm_set_skew(matcode);
169  else
170  return MM_UNSUPPORTED_TYPE;
171 
172 
173  return 0;
174 }
175 
176 int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
177 {
178  if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
180  else
181  return 0;
182 }
183 
184 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
185 {
186  char line[MM_MAX_LINE_LENGTH];
187  int num_items_read;
188 
189  /* set return null parameter values, in case we exit with errors */
190  *M = *N = *nz = 0;
191 
192  /* now continue scanning until you reach the end-of-comments */
193  do
194  {
195  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
196  return MM_PREMATURE_EOF;
197  }while (line[0] == '%');
198 
199  /* line[] is either blank or has M,N, nz */
200  if (sscanf(line, "%d %d %d", M, N, nz) == 3)
201  return 0;
202 
203  else
204  do
205  {
206  num_items_read = fscanf(f, "%d %d %d", M, N, nz);
207  if (num_items_read == EOF) return MM_PREMATURE_EOF;
208  }
209  while (num_items_read != 3);
210 
211  return 0;
212 }
213 
214 
215 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
216 {
217  char line[MM_MAX_LINE_LENGTH];
218  int num_items_read;
219  int nz;
220 
221  /* set return null parameter values, in case we exit with errors */
222  *M = *N = 0;
223 
224  /* now continue scanning until you reach the end-of-comments */
225  do
226  {
227  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
228  return MM_PREMATURE_EOF;
229  }while (line[0] == '%');
230 
231  /* line[] is either blank or has M,N, nz */
232  if (sscanf(line, "%d %d %*d", M, N) == 2)
233  return 0;
234 
235  else /* we have a blank line */
236  do
237  {
238  num_items_read = fscanf(f, "%d %d %*d", M, N);
239  if (num_items_read == EOF) return MM_PREMATURE_EOF;
240  }
241  while (num_items_read != 2);
242 
243  return 0;
244 }
245 
246 int mm_write_mtx_array_size(FILE *f, int M, int N)
247 {
248  if (fprintf(f, "%d %d\n", M, N) != 2)
250  else
251  return 0;
252 }
253 
254 
255 
256 /*-------------------------------------------------------------------------*/
257 
258 /******************************************************************/
259 /* use when I[], J[], and val[]J, and val[] are already allocated */
260 /******************************************************************/
261 
262 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
263  double val[], MM_typecode matcode)
264 {
265  int i;
266  if (mm_is_complex(matcode))
267  {
268  for (i=0; i<nz; i++)
269  if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
270  != 4) return MM_PREMATURE_EOF;
271  }
272  else if (mm_is_real(matcode))
273  {
274  for (i=0; i<nz; i++)
275  {
276  if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
277  != 3) return MM_PREMATURE_EOF;
278 
279  }
280  }
281 
282  else if (mm_is_pattern(matcode))
283  {
284  for (i=0; i<nz; i++)
285  if (fscanf(f, "%d %d", &I[i], &J[i])
286  != 2) return MM_PREMATURE_EOF;
287  }
288  else
289  return MM_UNSUPPORTED_TYPE;
290 
291  return 0;
292 
293 }
294 
295 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
296  double *real, double *imag, MM_typecode matcode)
297 {
298  if (mm_is_complex(matcode))
299  {
300  if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
301  != 4) return MM_PREMATURE_EOF;
302  }
303  else if (mm_is_real(matcode))
304  {
305  if (fscanf(f, "%d %d %lg\n", I, J, real)
306  != 3) return MM_PREMATURE_EOF;
307 
308  }
309 
310  else if (mm_is_pattern(matcode))
311  {
312  if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
313  }
314  else
315  return MM_UNSUPPORTED_TYPE;
316 
317  return 0;
318 
319 }
320 
321 
322 /************************************************************************
323  mm_read_mtx_crd() fills M, N, nz, array of values, and return
324  type code, e.g. 'MCRS'
325 
326  if matrix is complex, values[] is of size 2*nz,
327  (nz pairs of real/imaginary values)
328 ************************************************************************/
329 
330 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
331  double **val, MM_typecode *matcode)
332 {
333  int ret_code;
334  FILE *f;
335 
336  if (strcmp(fname, "stdin") == 0) f=stdin;
337  else
338  if ((f = fopen(fname, "r")) == NULL)
339  return MM_COULD_NOT_READ_FILE;
340 
341 
342  if ((ret_code = mm_read_banner(f, matcode)) != 0)
343  return ret_code;
344 
345  if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
346  mm_is_matrix(*matcode)))
347  return MM_UNSUPPORTED_TYPE;
348 
349  if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
350  return ret_code;
351 
352 
353  *I = (int *) malloc(*nz * sizeof(int));
354  *J = (int *) malloc(*nz * sizeof(int));
355  *val = NULL;
356 
357  if (mm_is_complex(*matcode))
358  {
359  *val = (double *) malloc(*nz * 2 * sizeof(double));
360  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
361  *matcode);
362  if (ret_code != 0) return ret_code;
363  }
364  else if (mm_is_real(*matcode))
365  {
366  *val = (double *) malloc(*nz * sizeof(double));
367  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
368  *matcode);
369  if (ret_code != 0) return ret_code;
370  }
371 
372  else if (mm_is_pattern(*matcode))
373  {
374  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
375  *matcode);
376  if (ret_code != 0) return ret_code;
377  }
378 
379  if (f != stdin) fclose(f);
380  return 0;
381 }
382 
383 int mm_write_banner(FILE *f, MM_typecode matcode)
384 {
385  char *str = mm_typecode_to_str(matcode);
386  int ret_code;
387 
388  ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
389  free(str);
390  if (ret_code !=2 )
392  else
393  return 0;
394 }
395 
396 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
397  double val[], MM_typecode matcode)
398 {
399  FILE *f;
400  int i;
401 
402  if (strcmp(fname, "stdout") == 0)
403  f = stdout;
404  else
405  if ((f = fopen(fname, "w")) == NULL)
407 
408  /* print banner followed by typecode */
409  fprintf(f, "%s ", MatrixMarketBanner);
410  fprintf(f, "%s\n", mm_typecode_to_str(matcode));
411 
412  /* print matrix sizes and nonzeros */
413  fprintf(f, "%d %d %d\n", M, N, nz);
414 
415  /* print values */
416  if (mm_is_pattern(matcode))
417  for (i=0; i<nz; i++)
418  fprintf(f, "%d %d\n", I[i], J[i]);
419  else
420  if (mm_is_real(matcode))
421  for (i=0; i<nz; i++)
422  fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
423  else
424  if (mm_is_complex(matcode))
425  for (i=0; i<nz; i++)
426  fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i],
427  val[2*i+1]);
428  else
429  {
430  if (f != stdout) fclose(f);
431  return MM_UNSUPPORTED_TYPE;
432  }
433 
434  if (f !=stdout) fclose(f);
435 
436  return 0;
437 }
438 
439 
441 {
443  char *types[4];
444  int error =0;
445 
446  /* check for MTX type */
447  if (mm_is_matrix(matcode))
448  types[0] = MM_MTX_STR;
449  else
450  error=1;
451 
452  /* check for CRD or ARR matrix */
453  if (mm_is_sparserow(matcode))
454  types[1] = MM_SPARSEROW_STR;
455  else
456  if (mm_is_coordinate(matcode))
457  types[1] = MM_COORDINATE_STR;
458  else
459  if (mm_is_dense(matcode))
460  types[1] = MM_DENSE_STR;
461  else
462  return NULL;
463 
464  /* check for element data type */
465  if (mm_is_real(matcode))
466  types[2] = MM_REAL_STR;
467  else
468  if (mm_is_complex(matcode))
469  types[2] = MM_COMPLEX_STR;
470  else
471  if (mm_is_pattern(matcode))
472  types[2] = MM_PATTERN_STR;
473  else
474  if (mm_is_integer(matcode))
475  types[2] = MM_INT_STR;
476  else
477  return NULL;
478 
479 
480  /* check for symmetry type */
481  if (mm_is_general(matcode))
482  types[3] = MM_GENERAL_STR;
483  else
484  if (mm_is_symmetric(matcode))
485  types[3] = MM_SYMM_STR;
486  else
487  if (mm_is_hermitian(matcode))
488  types[3] = MM_HERM_STR;
489  else
490  if (mm_is_skew(matcode))
491  types[3] = MM_SKEW_STR;
492  else
493  return NULL;
494 
495  sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
496  return strdup(buffer);
497 
498 }
499 
500 
501 #include "matrix_io.h"
502 
503 void csr2csc(int n, int m, int nz, double *a, int *col_idx, int *row_start,
504  double *csc_a, int *row_idx, int *col_start);
505 void coo2csr_in(int n, int nz, double *a, int *i_idx, int *j_idx);
506 void coo2csr(int n, int nz, double *a, int *i_idx, int *j_idx,
507  double *csr_a, int *col_idx, int *row_start);
508 
509 
510 /* write CSR format */
511 /* 1st line : % number_of_rows number_of_columns number_of_nonzeros
512  2nd line : % base of index
513  3rd line : row_number nz_r(=number_of_nonzeros_in_the_row)
514  next nz_r lines : column_index value(when a != NULL)
515  next line : row_number nz_r(=number_of_nonzeros_in_the_row)
516  next nz_r lines : column_index value(when a != NULL)
517  ...
518  */
519 
520 void write_csr (char *fn, int m, int n, int nz,
521  int *row_start, int *col_idx, double *a)
522 {
523  FILE *f;
524  int i, j;
525 
526  if ((f = fopen(fn, "w")) == NULL){
527  printf ("can't open file <%s> \n", fn);
528  exit(1);
529  }
530 
531  fprintf (f, "%s %d %d %d\n", "%", m, n, nz);
532  fprintf (f, "%s 0\n", "%");
533 
534  for (i=0; i<m; i++){
535  fprintf(f, "%d %d\n", i, row_start[i+1]-row_start[i]);
536 
537  for (j=row_start[i]; j<row_start[i+1]; j++){
538  if (a)
539  fprintf(f, "%d %20.19g\n", col_idx[j], a[j]);
540  else
541  fprintf(f, "%d\n", col_idx[j]);
542  }
543  }
544 
545  fclose (f);
546 
547 }
548 
549 
550 /* reads matrix market format (coordinate) and returns
551  csr format */
552 
553 void read_mm_matrix (char *fn, int *m, int *n, int *nz,
554  int **i_idx, int **j_idx, double **a)
555 {
556  MM_typecode matcode;
557  FILE *f;
558  int i,k;
559  int base=1, base2;
560 
561  if ((f = fopen(fn, "r")) == NULL) {
562  printf ("can't open file <%s> \n", fn);
563  exit(1);
564  }
565  if (mm_read_banner(f, &matcode) != 0){
566  printf("Could not process Matrix Market banner.\n");
567  exit(1);
568  }
569 
570  /* This is how one can screen matrix types if their application */
571  /* only supports a subset of the Matrix Market data types. */
572 
573  if (! (mm_is_matrix(matcode) && mm_is_sparse(matcode)) ){
574  printf("Sorry, this application does not support ");
575  printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
576  exit(1);
577  }
578 
579  /* find out size of sparse matrix .... */
580 
581  fscanf(f, "%*s %d %d %d", m, n, nz);
582 
583 // fscanf(f, "%d %d",&base2, &base); /* indentifies starting index */
584 
585  /* reserve memory for matrices */
586  /* if (mm_is_symmetric(matcode)){
587  *i_idx = (int *) my_malloc(*nz *2 * sizeof(int));
588  *j_idx = (int *) my_malloc(*nz *2 * sizeof(int));
589  *a = (double *) my_malloc(*nz *2 * sizeof(double));
590  }
591  else {*/
592  *i_idx = (int *) my_malloc(*nz * sizeof(int));
593  *j_idx = (int *) my_malloc(*nz * sizeof(int));
594  *a = (double *) my_malloc(*nz * sizeof(double));
595 
596  for (i=0; i<*nz; i++) {
597  (*i_idx)[i]=-1;
598  (*j_idx)[i]=-1;
599  (*a)[i] =-1;
600  }
601  // }
602 
603  if (!(*i_idx) || !(*j_idx) || !(*a)){
604  printf ("cannot allocate memory for %d, %d, %d sparse matrix\n", *m, *n, *nz);
605  exit(1);
606  }
607 
608 //printf("\n n:%d m: %d nz: %d base:%d\n", *m,*n,*nz, base);
609  k=0;
610 
611 //(*i_idx)[0]=base2-base;
612 //(*j_idx)[0]=base-base;
613  for (i=0; i<*nz; i++) {
614  //if (mm_is_pattern(matcode)){
615  fscanf(f, "%d %d", &(*i_idx)[i], &(*j_idx)[i]);
616  (*i_idx)[i] -= base; /* adjust from 1-based to 0-based */
617  (*j_idx)[i] -= base;
618  //printf("\n %d %d", (*i_idx)[i],(*j_idx)[i]);
619  (*a)[i] = random_double(-1, 1);
620  // }
621  /* else if (mm_is_real(matcode)){
622  fscanf(f, "%d %d %lg", &(*i_idx)[i], &(*j_idx)[i], &(*a)[i]);
623  (*i_idx)[i] -= base; // adjust from 1-based to 0-based
624  (*j_idx)[i] -= base;
625  }*/
626 
627  /* if (mm_is_symmetric(matcode)){
628  if ( (*i_idx)[i] != (*j_idx)[i] ){
629  (*i_idx)[*nz+k] = (*j_idx)[i];
630  (*j_idx)[*nz+k] = (*i_idx)[i];
631  (*a)[*nz+k] = (*a)[i];
632  k++;
633  }
634  }*/
635  }
636  *nz += k;
637 
638  fclose(f);
639  for (i=0; i<*nz; i++) {
640  //printf("\n %d %d ",(*i_idx)[i], (*j_idx)[i] );
641  }
642  coo2csr_in (*m, *nz, *a, *i_idx, *j_idx);
643 }
644 
645 /* reads Harwell-Boeing format matrix and returns CSR format */
646 
647 #define LINE_LEN 90
648 #define INPUT_WIDTH 80
649 
650 void read_hb_matrix (char *fn, int *m, int *n, int *nz,
651  int **row_start, int **col_idx, double **a)
652 {
653  MM_typecode matcode;
654  FILE *f;
655  int i, field, j, k;
656  int row, rhs;
657  char buffer[LINE_LEN];
658  char mat_format[4];
659  int *i_idx, *j_idx;
660  int *col_start, *row_idx;
661  double *coo_a;
662  char start_format[INPUT_WIDTH];
663  char idx_format[INPUT_WIDTH];
664  int start_input_repeat, start_input_width;
665  int idx_input_repeat, idx_input_width;
666  int start_lines, idx_lines, l;
667 
668  if ((f = fopen(fn, "r")) == NULL) {
669  printf ("can't open file <%s> \n", fn);
670  exit(1);
671  }
672  fgets (buffer, LINE_LEN, f); /* title line */
673  fgets (buffer, LINE_LEN, f); /* line counts */
674  sscanf(buffer, "%*d %d %d %*d %d", &start_lines, &idx_lines, &rhs);
675  fgets (buffer, LINE_LEN, f); /* elements counts */
676  sscanf(buffer, "%s %d %d %d", mat_format, m, n, nz);
677  printf("type=%s m=%d n=%d nz=%d rhs=%d\n", mat_format, *m, *n, *nz, rhs);
678  fgets (buffer, LINE_LEN, f); /* formats */
679  sscanf(buffer, "%s %s", start_format, idx_format);
680  sscanf(start_format, "(%d%*c%d)", &start_input_repeat, &start_input_width);
681  sscanf(idx_format, "(%d%*c%d)", &idx_input_repeat, &idx_input_width);
682  printf("%d start_input lines %d repeats %d: %d col_input lines %d repeats %d\n",
683  start_lines, start_input_width, start_input_repeat,
684  idx_lines, idx_input_width, idx_input_repeat);
685 
686  if (rhs)
687  fgets (buffer, LINE_LEN, f); /* rhs header */
688 
689  mm_clear_typecode(&matcode);
690  mm_set_matrix(&matcode);
691  mm_set_sparserow(&matcode);
692 
693  switch (mat_format[0]){
694  case 'r':
695  case 'R':
696  mm_set_real(&matcode);
697  break;
698  case 'c':
699  case 'C':
700  mm_set_complex(&matcode);
701  break;
702  case 'p':
703  case 'P':
704  mm_set_pattern(&matcode);
705  break;
706  default:
707  printf ("<%s> is not Harwell Boeing matrix format %s\n", fn, mat_format);
708  exit (0);
709  }
710 
711  switch (mat_format[1]){
712  case 'u':
713  case 'U':
714  case 'r':
715  case 'R':
716  mm_set_general(&matcode);
717  break;
718  case 's':
719  case 'S':
720  mm_set_symmetric(&matcode);
721  break;
722  case 'h':
723  case 'H':
724  mm_set_hermitian(&matcode);
725  break;
726  case 'z':
727  case 'Z':
728  mm_set_skew(&matcode);
729  break;
730  default:
731  printf ("<%s> is not Harwell Boeing matrix format %s\n", fn, mat_format);
732  exit (0);
733  }
734 
735  switch (mat_format[2]){
736  case 'a':
737  case 'A':
738  break;
739  case 'u':
740  case 'U':
741  printf ("<%s> is not unassembled Harwell Boeing matrix format %s\n",
742  fn, mat_format);
743  exit (0);
744  default:
745  printf ("<%s> is not Harwell Boeing matrix format %s\n", fn, mat_format);
746  exit (0);
747  }
748 
749 
750  *row_start = (int *) my_malloc((*m+1) * sizeof(int));
751 
752  /* reserve memory for matrices */
753  if (mm_is_symmetric(matcode)){
754  i_idx = (int *) my_malloc(*nz *2 * sizeof(int));
755  j_idx = (int *) my_malloc(*nz *2 * sizeof(int));
756  coo_a = (double *) my_malloc(*nz *2 * sizeof(double));
757  }
758  else {
759  col_start = (int *) my_malloc((*n+1) * sizeof(int));
760  row_idx = (int *) my_malloc(*nz * sizeof(int));
761 
762  *col_idx = (int *) my_malloc(*nz * sizeof(int));
763  *a = (double *) my_malloc(*nz * sizeof(double));
764  }
765 
766 
767  if (mm_is_symmetric(matcode)){
768  for (i=l=0; l<start_lines; l++){
769  fgets (buffer, LINE_LEN, f);
770  for (field=0; (field<start_input_repeat) && (i<=*m); field++){
771  for(j=0; j<start_input_width; j++)
772  start_format[j] = buffer[field*start_input_width+j];
773  start_format[j] = 0;
774 
775  (*row_start)[i++] = atoi(start_format)-1;
776  }
777  }
778 
779  for (i=l=k=row=0; l<idx_lines; l++){
780  fgets (buffer, LINE_LEN, f);
781  for (field=0; (field<idx_input_repeat) && (i<*nz); field++, i++){
782  for(j=0; j<idx_input_width; j++)
783  idx_format[j] = buffer[field*idx_input_width+j];
784  idx_format[j] = 0;
785 
786  if (i==(*row_start)[row+1]) row++;
787  i_idx[i] = row;
788  j_idx[i] = atoi(idx_format)-1;
789 
790  coo_a[i] = random_double(-1, 1);
791 
792  if ( i_idx[i] != j_idx[i] ){
793  i_idx[*nz+k] = j_idx[i];
794  j_idx[*nz+k] = i_idx[i];
795  coo_a[*nz+k] = coo_a[i];
796  k++;
797  }
798  }
799  }
800  *nz += k;
801 
802  *col_idx = (int *) my_malloc(*nz * sizeof(int));
803  *a = (double *) my_malloc(*nz * sizeof(double));
804  coo2csr (*m, *nz, coo_a, i_idx, j_idx, *a, *col_idx, *row_start);
805 
806  free (i_idx);
807  free (j_idx);
808  free (coo_a);
809 
810  }
811  else {
812  for (i=l=0; l<start_lines; l++){
813  fgets (buffer, LINE_LEN, f);
814  for (field=0; (field<start_input_repeat) && (i<=*n); field++){
815  for(j=0; j<start_input_width; j++)
816  start_format[j] = buffer[field*start_input_width+j];
817  start_format[j] = 0;
818 
819  col_start[i++] = atoi(start_format)-1;
820  }
821  }
822 
823 
824  for (i=l=0; l<idx_lines; l++){
825  fgets (buffer, LINE_LEN, f);
826  for (field=0; (field<idx_input_repeat) && (i<*nz); field++, i++){
827  for(j=0; j<idx_input_width; j++)
828  idx_format[j] = buffer[field*idx_input_width+j];
829  idx_format[j] = 0;
830 
831  row_idx[i] = atoi(idx_format)-1;
832 
833  (*a)[i] = random_double(-1, 1);
834  }
835  }
836 
837  csr2csc(*n, *m, *nz, NULL, row_idx, col_start, NULL, *col_idx, *row_start);
838 
839 
840  free (row_idx);
841  free (col_start);
842  }
843 
844  fclose(f);
845 }
846 
847 
848 void sort(int *col_idx, double *a, int start, int end)
849 {
850  int i, j, it;
851  double dt;
852 
853  for (i=end-1; i>start; i--)
854  for(j=start; j<i; j++)
855  if (col_idx[j] > col_idx[j+1]){
856 
857  if (a){
858  dt=a[j];
859  a[j]=a[j+1];
860  a[j+1]=dt;
861  }
862  it=col_idx[j];
863  col_idx[j]=col_idx[j+1];
864  col_idx[j+1]=it;
865 
866  }
867 }
868 
869 
870 
871 /* converts COO format to CSR format, in-place,
872  if SORT_IN_ROW is defined, each row is sorted in column index.
873  On return, i_idx contains row_start position */
874 
875 void coo2csr_in(int n, int nz, double *a, int *i_idx, int *j_idx)
876 {
877  int *row_start;
878  int i, j;
879  int init, i_next, j_next, i_pos;
880  double dt, a_next;
881 
882  row_start = (int *)malloc((n+1)*sizeof(int));
883  if (!row_start){
884  printf ("coo2csr_in: cannot allocate temporary memory\n");
885  exit (1);
886  }
887  for (i=0; i<=n; i++) row_start[i] = 0;
888 
889 #ifdef DEBUG
890  printf("\n idx \n" );
891  for (i=0; i<nz; i++){
892  printf("\n %d", i_idx[i]);
893  }
894 
895 #endif
896  /* determine row lengths */
897  for (i=0; i<nz; i++) {
898  if(i_idx[i]!=-1){
899  row_start[i_idx[i]+1]++;
900  }
901  }
902 
903 
904 #ifdef DEBUG
905  printf("\n row Start \n" );
906  for (i=0; i<n; i++){
907  printf("\n %d", row_start[i]);
908  }
909 #endif
910 
911  for (i=0; i<n; i++) row_start[i+1] += row_start[i];
912 
913 
914 #ifdef DEBUG
915  printf("\n row Start \n" );
916  for (i=0; i<n; i++){
917 
918 printf("\n %d", row_start[i]);
919 }
920 #endif
921 
922 
923  for (init=0; init<nz; ){
924  dt = a[init];
925  i = i_idx[init];
926  j = j_idx[init];
927  i_idx[init] = -1;
928 
929  while (1){
930  i_pos = row_start[i];
931  a_next = a[i_pos];
932  i_next = i_idx[i_pos];
933  j_next = j_idx[i_pos];
934 
935  a[i_pos] = dt;
936  j_idx[i_pos] = j;
937  i_idx[i_pos] = -1;
938  row_start[i]++;
939  if (i_next < 0) break;
940  dt = a_next;
941  i = i_next;
942  j = j_next;
943 
944  }
945  init++;
946  while ((i_idx[init] < 0) && (init < nz)) init++;
947  }
948 
949 
950  /* shift back row_start */
951  for (i=0; i<n; i++) i_idx[i+1] = row_start[i];
952  i_idx[0] = 0;
953 
954 #ifdef DEBUG
955  printf("\n row Start \n" );
956  for (i=0; i<n; i++){
957 
958 printf("\n %d", i_idx[i]);
959 }
960 #endif
961  for (i=0; i<n; i++){
962  sort (j_idx, a, i_idx[i], i_idx[i+1]);
963  }
964 
965 
966 }
967 
968 /* converts COO format to CSR format, not in-place,
969  if SORT_IN_ROW is defined, each row is sorted in column index */
970 
971 
972 void coo2csr(int n, int nz, double *a, int *i_idx, int *j_idx,
973  double *csr_a, int *col_idx, int *row_start)
974 {
975  int i, l;
976 
977  for (i=0; i<=n; i++) row_start[i] = 0;
978 
979  /* determine row lengths */
980  for (i=0; i<nz; i++) row_start[i_idx[i]+1]++;
981 
982 
983  for (i=0; i<n; i++) row_start[i+1] += row_start[i];
984 
985 
986  /* go through the structure once more. Fill in output matrix. */
987  for (l=0; l<nz; l++){
988  i = row_start[i_idx[l]];
989  csr_a[i] = a[l];
990  col_idx[i] = j_idx[l];
991  row_start[i_idx[l]]++;
992  }
993 
994  /* shift back row_start */
995  for (i=n; i>0; i--) row_start[i] = row_start[i-1];
996 
997  row_start[0] = 0;
998 
999  for (i=0; i<n; i++){
1000  sort (col_idx, csr_a, row_start[i], row_start[i+1]);
1001  }
1002 
1003 }
1004 
1005 
1006 /*
1007  converts CSR format to CSC format, not in-place,
1008  if a == NULL, only pattern is reorganized.
1009  the size of matrix is n x m.
1010  */
1011 
1012 void csr2csc(int n, int m, int nz, double *a, int *col_idx, int *row_start,
1013  double *csc_a, int *row_idx, int *col_start)
1014 {
1015  int i, j, k, l;
1016  int *ptr;
1017 
1018  for (i=0; i<=m; i++) col_start[i] = 0;
1019 
1020  /* determine column lengths */
1021  for (i=0; i<nz; i++) col_start[col_idx[i]+1]++;
1022 
1023  for (i=0; i<m; i++) col_start[i+1] += col_start[i];
1024 
1025 
1026  /* go through the structure once more. Fill in output matrix. */
1027 
1028  for (i=0, ptr=row_start; i<n; i++, ptr++)
1029  for (j=*ptr; j<*(ptr+1); j++){
1030  k = col_idx[j];
1031  l = col_start[k]++;
1032  row_idx[l] = i;
1033  if (a) csc_a[l] = a[j];
1034  }
1035 
1036  /* shift back col_start */
1037  for (i=m; i>0; i--) col_start[i] = col_start[i-1];
1038 
1039  col_start[0] = 0;
1040 }
#define MM_PATTERN_STR
Definition: mmio.h:121
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
Definition: matrix_io.c:176
#define mm_is_pattern(typecode)
Definition: mmio.h:45
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
Definition: matrix_io.c:396
#define MM_MAX_LINE_LENGTH
Definition: mmio.h:16
#define MM_SYMM_STR
Definition: mmio.h:118
#define mm_clear_typecode(typecode)
Definition: mmio.h:75
#define mm_set_real(typecode)
Definition: mmio.h:65
#define mm_set_pattern(typecode)
Definition: mmio.h:66
#define MM_GENERAL_STR
Definition: mmio.h:117
#define mm_set_coordinate(typecode)
Definition: mmio.h:59
#define MM_UNSUPPORTED_TYPE
Definition: mmio.h:88
#define LINE_LEN
Definition: matrix_io.c:647
#define mm_set_sparserow(typecode)
Definition: mmio.h:60
int random_integer(int low, int high)
Definition: matrix_io.c:28
#define mm_set_complex(typecode)
Definition: mmio.h:64
int mm_write_banner(FILE *f, MM_typecode matcode)
Definition: matrix_io.c:383
#define MM_PREMATURE_EOF
Definition: mmio.h:85
#define mm_set_matrix(typecode)
Definition: mmio.h:58
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
Definition: matrix_io.c:215
#define mm_set_dense(typecode)
Definition: mmio.h:62
#define mm_is_dense(typecode)
Definition: mmio.h:40
#define mm_is_sparserow(typecode)
Definition: mmio.h:38
#define MM_COORDINATE_STR
Definition: mmio.h:112
void read_hb_matrix(char *fn, int *m, int *n, int *nz, int **row_start, int **col_idx, double **a)
Definition: matrix_io.c:650
#define mm_is_skew(typecode)
Definition: mmio.h:50
#define r(a, b, c)
#define MatrixMarketBanner
Definition: mmio.h:17
MPI_Datatype * types
Definition: refMPIluBack.c:109
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
Definition: matrix_io.c:262
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)
Definition: matrix_io.c:295
void coo2csr(int n, int nz, double *a, int *i_idx, int *j_idx, double *csr_a, int *col_idx, int *row_start)
Definition: matrix_io.c:972
#define mm_set_hermitian(typecode)
Definition: mmio.h:73
#define MM_DENSE_STR
Definition: mmio.h:111
#define MM_MAX_TOKEN_LENGTH
Definition: mmio.h:18
void write_csr(char *fn, int m, int n, int nz, int *row_start, int *col_idx, double *a)
Definition: matrix_io.c:520
#define MM_INT_STR
Definition: mmio.h:116
int mm_write_mtx_array_size(FILE *f, int M, int N)
Definition: matrix_io.c:246
double random_double(double low, double high)
Definition: matrix_io.c:38
#define mm_set_skew(typecode)
Definition: mmio.h:72
#define mm_set_general(typecode)
Definition: mmio.h:71
void coo2csr_in(int n, int nz, double *a, int *i_idx, int *j_idx)
Definition: matrix_io.c:875
char * mm_typecode_to_str(MM_typecode matcode)
Definition: matrix_io.c:440
#define mm_is_general(typecode)
Definition: mmio.h:49
#define mm_is_real(typecode)
Definition: mmio.h:44
int mm_read_banner(FILE *f, MM_typecode *matcode)
Definition: matrix_io.c:88
#define mm_is_coordinate(typecode)
Definition: mmio.h:39
#define mm_is_sparse(typecode)
Definition: mmio.h:37
#define mm_is_hermitian(typecode)
Definition: mmio.h:51
#define MM_MTX_STR
Definition: mmio.h:109
void sort(int *col_idx, double *a, int start, int end)
Definition: matrix_io.c:848
#define MM_COULD_NOT_READ_FILE
Definition: mmio.h:84
void read_mm_matrix(char *fn, int *m, int *n, int *nz, int **i_idx, int **j_idx, double **a)
Definition: matrix_io.c:553
#define MM_COMPLEX_STR
Definition: mmio.h:114
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, double **val, MM_typecode *matcode)
Definition: matrix_io.c:330
int str_to_mem_unit(char *str)
Definition: matrix_io.c:57
#define MM_SKEW_STR
Definition: mmio.h:120
#define m(a, b, c)
#define mm_is_integer(typecode)
Definition: mmio.h:46
#define MM_REAL_STR
Definition: mmio.h:115
#define MM_HERM_STR
Definition: mmio.h:119
char MM_typecode[4]
Definition: mmio.h:20
#define MM_NO_HEADER
Definition: mmio.h:87
#define mm_is_matrix(typecode)
Definition: mmio.h:35
int mm_is_valid(MM_typecode matcode)
Definition: matrix_io.c:78
#define MM_SPARSEROW_STR
Definition: mmio.h:113
#define v(a, b, c)
#define mm_set_symmetric(typecode)
Definition: mmio.h:70
#define MM_COULD_NOT_WRITE_FILE
Definition: mmio.h:90
#define mm_is_complex(typecode)
Definition: mmio.h:43
void csr2csc(int n, int m, int nz, double *a, int *col_idx, int *row_start, double *csc_a, int *row_idx, int *col_start)
Definition: matrix_io.c:1012
void * my_malloc(int sz)
Definition: matrix_io.c:12
#define mm_is_symmetric(typecode)
Definition: mmio.h:48
double drand48()
double ** buffer
Definition: refMPIluBack.c:103
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
Definition: matrix_io.c:184
#define INPUT_WIDTH
Definition: matrix_io.c:648
#define mm_set_integer(typecode)
Definition: mmio.h:67