Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_MMWrite: write a matrix to a Matrix Market file
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-License-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
10 : // Contributors.txt for a full list of contributors. Created, in part, with
11 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
12 : // DM22-0790
13 :
14 : // Contributed by Timothy A. Davis, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // LAGraph_MMWrite: write a matrix to a Matrix Market file.
19 :
20 : // Writes a matrix to a file in the Matrix Market format. See LAGraph_MMRead
21 : // for a description of the format.
22 :
23 : // The Matrix Market format is described at:
24 : // https://math.nist.gov/MatrixMarket/formats.html
25 :
26 : // Parts of this code are from SuiteSparse/CHOLMOD/Check/cholmod_write.c, and
27 : // are used here by permission of the author of CHOLMOD/Check (T. A. Davis).
28 :
29 : #include "LG_internal.h"
30 :
31 : #undef LG_FREE_WORK
32 : #define LG_FREE_WORK \
33 : { \
34 : LAGraph_Free ((void **) &I, NULL) ; \
35 : LAGraph_Free ((void **) &J, NULL) ; \
36 : LAGraph_Free ((void **) &K, NULL) ; \
37 : LAGraph_Free ((void **) &X, NULL) ; \
38 : GrB_free (&AT) ; \
39 : GrB_free (&M) ; \
40 : GrB_free (&C) ; \
41 : }
42 :
43 : #undef LG_FREE_ALL
44 : #define LG_FREE_ALL LG_FREE_WORK
45 :
46 : //------------------------------------------------------------------------------
47 : // print_double
48 : //------------------------------------------------------------------------------
49 :
50 : // Print a double value to the file, using the shortest format that ensures the
51 : // value is written precisely. Returns true if successful, false if an I/O
52 : // error occurred.
53 :
54 120798 : static bool print_double
55 : (
56 : FILE *f, // file to print to
57 : double x // value to print
58 : )
59 : {
60 :
61 : char s [MAXLINE], *p ;
62 120798 : int64_t i, dest = 0, src = 0 ;
63 : int width, ok ;
64 :
65 : //--------------------------------------------------------------------------
66 : // handle Inf and NaN
67 : //--------------------------------------------------------------------------
68 :
69 120798 : if (isnan (x))
70 : {
71 1 : return (fprintf (f, "nan") > 0) ;
72 : }
73 120797 : if (isinf (x))
74 : {
75 12 : return (fprintf (f, (x < 0) ? "-inf" : "inf") > 0) ;
76 : }
77 :
78 : //--------------------------------------------------------------------------
79 : // find the smallest acceptable precision
80 : //--------------------------------------------------------------------------
81 :
82 869437 : for (width = 6 ; width < 20 ; width++)
83 : {
84 : double y ;
85 869437 : sprintf (s, "%.*g", width, x) ;
86 869437 : sscanf (s, "%lg", &y) ;
87 869437 : if (x == y) break ;
88 : }
89 :
90 : //--------------------------------------------------------------------------
91 : // shorten the string
92 : //--------------------------------------------------------------------------
93 :
94 : // change "e+0" to "e", change "e+" to "e", and change "e-0" to "e-"
95 1774378 : for (i = 0 ; i < MAXLINE && s [i] != '\0' ; i++)
96 : {
97 1671099 : if (s [i] == 'e')
98 : {
99 17506 : if (s [i+1] == '+')
100 : {
101 3532 : dest = i+1 ;
102 3532 : if (s [i+2] == '0')
103 : {
104 : // delete characters s[i+1] and s[i+2]
105 3324 : src = i+3 ;
106 : }
107 : else
108 : {
109 : // delete characters s[i+1]
110 208 : src = i+2 ;
111 : }
112 : }
113 13974 : else if (s [i+1] == '-')
114 : {
115 13974 : dest = i+2 ;
116 13974 : if (s [i+2] == '0')
117 : {
118 : // delete character s[i+2]
119 12258 : src = i+3 ;
120 : }
121 : else
122 : {
123 : // no change
124 1716 : break ;
125 : }
126 : }
127 31788 : while (s [src] != '\0')
128 : {
129 15998 : s [dest++] = s [src++] ;
130 : }
131 15790 : s [dest] = '\0' ;
132 15790 : break ;
133 : }
134 : }
135 :
136 : // delete the leading "0" if present and not necessary
137 120785 : p = s ;
138 120785 : s [MAXLINE-1] = '\0' ;
139 120785 : i = strlen (s) ;
140 120785 : if (i > 2 && s [0] == '0' && s [1] == '.')
141 : {
142 : // change "0.x" to ".x"
143 9359 : p = s + 1 ;
144 : }
145 111426 : else if (i > 3 && s [0] == '-' && s [1] == '0' && s [2] == '.')
146 : {
147 : // change "-0.x" to "-.x"
148 4528 : s [1] = '-' ;
149 4528 : p = s + 1 ;
150 : }
151 :
152 : #if 0
153 : // double-check
154 : i = sscanf (p, "%lg", &z) ;
155 : if (i != 1 || y != z)
156 : {
157 : // oops! something went wrong in the "e+0" edit, above.
158 : // this "cannot" happen
159 : sprintf (s, "%.*g", width, x) ;
160 : p = s ;
161 : }
162 : #endif
163 :
164 : //--------------------------------------------------------------------------
165 : // print the value to the file
166 : //--------------------------------------------------------------------------
167 :
168 120785 : return (fprintf (f, "%s", p) > 0) ;
169 : }
170 :
171 : //------------------------------------------------------------------------------
172 : // LAGraph_MMWrite: write a matrix to a MatrixMarket file
173 : //------------------------------------------------------------------------------
174 :
175 1367 : int LAGraph_MMWrite
176 : (
177 : // input:
178 : GrB_Matrix A, // matrix to write to the file
179 : FILE *f, // file to write it to, must be already open
180 : FILE *fcomments, // optional file with extra comments, may be NULL
181 : char *msg
182 : )
183 : {
184 :
185 : //--------------------------------------------------------------------------
186 : // check inputs
187 : //--------------------------------------------------------------------------
188 :
189 1367 : LG_CLEAR_MSG ;
190 1367 : void *X = NULL ;
191 1367 : GrB_Index *I = NULL, *J = NULL, *K = NULL ;
192 1367 : GrB_Matrix M = NULL, AT = NULL, C = NULL ;
193 1367 : LG_ASSERT (A != NULL, GrB_NULL_POINTER) ;
194 1366 : LG_ASSERT (f != NULL, GrB_NULL_POINTER) ;
195 :
196 : //--------------------------------------------------------------------------
197 : // determine the basic matrix properties
198 : //--------------------------------------------------------------------------
199 :
200 : GrB_Index nrows, ncols, nvals ;
201 1366 : GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ;
202 1366 : GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
203 1366 : GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ;
204 1366 : GrB_Index n = nrows ;
205 :
206 : //--------------------------------------------------------------------------
207 : // determine if the matrix is dense
208 : //--------------------------------------------------------------------------
209 :
210 1366 : MM_fmt_enum MM_fmt = MM_coordinate ;
211 :
212 : // guard against integer overflow
213 1366 : if (((double) nrows * (double) ncols < (double) INT64_MAX) &&
214 1366 : (nvals == nrows * ncols))
215 : {
216 62 : MM_fmt = MM_array ;
217 : }
218 :
219 : //--------------------------------------------------------------------------
220 : // determine the entry type
221 : //--------------------------------------------------------------------------
222 :
223 1366 : GrB_Type type = NULL ;
224 : char atype_name [LAGRAPH_MAX_NAME_LEN] ;
225 1366 : atype_name [0] = '\0' ;
226 1366 : LAGraph_Matrix_TypeName (atype_name, A, msg) ;
227 1366 : LAGraph_TypeFromName (&type, atype_name, msg) ;
228 :
229 1366 : MM_type_enum MM_type = MM_integer ;
230 :
231 1366 : if (type == GrB_BOOL || type == GrB_INT8 || type == GrB_INT16 ||
232 909 : type == GrB_INT32 || type == GrB_INT64 || type == GrB_UINT8 ||
233 525 : type == GrB_UINT16 || type == GrB_UINT32 || type == GrB_UINT64)
234 : {
235 907 : MM_type = MM_integer ;
236 : }
237 459 : else if (type == GrB_FP32 || type == GrB_FP64)
238 : {
239 458 : MM_type = MM_real ;
240 : }
241 : #if 0
242 : else if (type == GxB_FC32 || type == GxB_FC64)
243 : {
244 : MM_type = MM_complex ;
245 : }
246 : #endif
247 : else
248 : {
249 1 : LG_ASSERT_MSG (false, GrB_NOT_IMPLEMENTED, "type not supported") ;
250 : }
251 :
252 : //--------------------------------------------------------------------------
253 : // determine symmetry
254 : //--------------------------------------------------------------------------
255 :
256 1365 : MM_storage_enum MM_storage = MM_general ;
257 :
258 1365 : if (nrows == ncols)
259 : {
260 : // AT = A'
261 1841 : GRB_TRY (GrB_Matrix_new (&AT, type, n, n)) ;
262 1190 : GRB_TRY (GrB_transpose (AT, NULL, NULL, A, NULL)) ;
263 :
264 : //----------------------------------------------------------------------
265 : // check for symmetry
266 : //----------------------------------------------------------------------
267 :
268 1084 : bool isequal = false ;
269 1084 : LG_TRY (LAGraph_Matrix_IsEqual (&isequal, A, AT, msg)) ;
270 816 : if (isequal)
271 : {
272 263 : MM_storage = MM_symmetric ;
273 : }
274 :
275 : //----------------------------------------------------------------------
276 : // check for skew-symmetry
277 : //----------------------------------------------------------------------
278 :
279 : // for signed types only
280 816 : if (MM_storage == MM_general)
281 : {
282 : // select the operator
283 553 : GrB_UnaryOp op = NULL ;
284 553 : if (type == GrB_INT8 ) op = GrB_AINV_INT8 ;
285 512 : else if (type == GrB_INT16) op = GrB_AINV_INT16 ;
286 471 : else if (type == GrB_INT32) op = GrB_AINV_INT32 ;
287 394 : else if (type == GrB_INT64) op = GrB_AINV_INT64 ;
288 292 : else if (type == GrB_FP32 ) op = GrB_AINV_FP32 ;
289 251 : else if (type == GrB_FP64 ) op = GrB_AINV_FP64 ;
290 : #if 0
291 : else if (type == GxB_FC32 ) op = GxB_AINV_FC32 ;
292 : else if (type == GxB_FC64 ) op = GxB_AINV_FC64 ;
293 : #endif
294 553 : if (op != NULL)
295 : {
296 456 : GRB_TRY (GrB_apply (AT, NULL, NULL, op, AT, NULL)) ;
297 455 : LG_TRY (LAGraph_Matrix_IsEqual (&isequal, A, AT, msg)) ;
298 320 : if (isequal)
299 : {
300 114 : MM_storage = MM_skew_symmetric ;
301 : }
302 : }
303 : }
304 :
305 : //----------------------------------------------------------------------
306 : // check for Hermitian (not yet supported)
307 : //----------------------------------------------------------------------
308 :
309 : #if 0
310 : if (MM_type == MM_complex && MM_storage == MM_general)
311 : {
312 : LG_TRY (LAGraph_Matrix_IsEqualOp (&isequal, A, AT,
313 : LAGraph_HERMITIAN_ComplexFP64, msg)) ;
314 : if (isequal)
315 : {
316 : MM_storage = MM_hermitian ;
317 : }
318 : }
319 : #endif
320 :
321 680 : GrB_free (&AT) ;
322 : }
323 :
324 : //--------------------------------------------------------------------------
325 : // determine if the matrix is structural-only
326 : //--------------------------------------------------------------------------
327 :
328 714 : bool is_structural = false ;
329 714 : if (! (MM_storage == MM_skew_symmetric || MM_storage == MM_hermitian))
330 : {
331 600 : if (type == GrB_BOOL)
332 : {
333 179 : GRB_TRY (GrB_reduce (&is_structural, NULL, GrB_LAND_MONOID_BOOL,
334 : A, NULL)) ;
335 : }
336 : else
337 : {
338 421 : GRB_TRY (GrB_Matrix_new (&C, GrB_BOOL, nrows, ncols)) ;
339 332 : GrB_BinaryOp op = NULL ;
340 332 : if (type == GrB_INT8 ) op = GrB_EQ_INT8 ;
341 323 : else if (type == GrB_INT16 ) op = GrB_EQ_INT16 ;
342 314 : else if (type == GrB_INT32 ) op = GrB_EQ_INT32 ;
343 282 : else if (type == GrB_INT64 ) op = GrB_EQ_INT64 ;
344 198 : else if (type == GrB_UINT8 ) op = GrB_EQ_UINT8 ;
345 189 : else if (type == GrB_UINT16) op = GrB_EQ_UINT16 ;
346 180 : else if (type == GrB_UINT32) op = GrB_EQ_UINT32 ;
347 171 : else if (type == GrB_UINT64) op = GrB_EQ_UINT64 ;
348 162 : else if (type == GrB_FP32 ) op = GrB_EQ_FP32 ;
349 153 : else if (type == GrB_FP64 ) op = GrB_EQ_FP64 ;
350 : #if 0
351 : else if (type == GxB_FC32 ) op = GrB_EQ_FC32 ;
352 : else if (type == GxB_FC64 ) op = GrB_EQ_FC64 ;
353 : #endif
354 332 : GRB_TRY (GrB_apply (C, NULL, NULL, op, A, 1, NULL)) ;
355 269 : GRB_TRY (GrB_reduce (&is_structural, NULL, GrB_LAND_MONOID_BOOL,
356 : C, NULL)) ;
357 269 : GrB_free (&C) ;
358 : }
359 448 : if (is_structural)
360 : {
361 225 : MM_type = MM_pattern ;
362 225 : MM_fmt = MM_coordinate ;
363 : }
364 : }
365 :
366 : //--------------------------------------------------------------------------
367 : // write the Matrix Market header
368 : //--------------------------------------------------------------------------
369 :
370 562 : FPRINTF (f, "%%%%MatrixMarket matrix") ;
371 :
372 562 : switch (MM_fmt)
373 : {
374 532 : default :
375 532 : case MM_coordinate : FPRINTF (f, " coordinate") ; break ;
376 30 : case MM_array : FPRINTF (f, " array") ; break ;
377 : }
378 :
379 562 : switch (MM_type)
380 : {
381 170 : default :
382 170 : case MM_real : FPRINTF (f, " real") ; break ;
383 167 : case MM_integer : FPRINTF (f, " integer") ; break ;
384 : // case MM_complex : FPRINTF (f, " complex") ; break ;
385 225 : case MM_pattern : FPRINTF (f, " pattern") ; break ;
386 : }
387 :
388 562 : switch (MM_storage)
389 : {
390 215 : default :
391 215 : case MM_general : FPRINTF (f, " general\n") ; break ;
392 233 : case MM_symmetric : FPRINTF (f, " symmetric\n") ; break ;
393 114 : case MM_skew_symmetric : FPRINTF (f, " skew-symmetric\n") ; break ;
394 : // case MM_hermitian : FPRINTF (f, " Hermitian\n") ; break ;
395 : }
396 :
397 562 : FPRINTF (f, "%%%%GraphBLAS type ") ;
398 562 : if (type == GrB_BOOL ) { FPRINTF (f, "bool\n") ; }
399 383 : else if (type == GrB_INT8 ) { FPRINTF (f, "int8_t\n") ; }
400 357 : else if (type == GrB_INT16 ) { FPRINTF (f, "int16_t\n") ; }
401 331 : else if (type == GrB_INT32 ) { FPRINTF (f, "int32_t\n") ; }
402 288 : else if (type == GrB_INT64 ) { FPRINTF (f, "int64_t\n") ; }
403 198 : else if (type == GrB_UINT8 ) { FPRINTF (f, "uint8_t\n") ; }
404 191 : else if (type == GrB_UINT16) { FPRINTF (f, "uint16_t\n") ; }
405 184 : else if (type == GrB_UINT32) { FPRINTF (f, "uint32_t\n") ; }
406 177 : else if (type == GrB_UINT64) { FPRINTF (f, "uint64_t\n") ; }
407 170 : else if (type == GrB_FP32 ) { FPRINTF (f, "float\n") ; }
408 144 : else if (type == GrB_FP64 ) { FPRINTF (f, "double\n") ; }
409 : #if 0
410 : else if (type == GxB_FC32 ) { FPRINTF (f, "float complex\n") ; }
411 : else if (type == GxB_FC64 ) { FPRINTF (f, "double complex\n") ; }
412 : #endif
413 :
414 : //--------------------------------------------------------------------------
415 : // include any additional comments
416 : //--------------------------------------------------------------------------
417 :
418 562 : if (fcomments != NULL)
419 : {
420 : char buffer [MAXLINE] ;
421 12 : while (fgets (buffer, MAXLINE-1, fcomments) != NULL)
422 : {
423 8 : FPRINTF (f, "%%%s", buffer) ;
424 : }
425 : }
426 :
427 : //--------------------------------------------------------------------------
428 : // print the first line
429 : //--------------------------------------------------------------------------
430 :
431 562 : bool is_general = (MM_storage == MM_general) ;
432 562 : GrB_Index nvals_to_print = nvals ;
433 :
434 562 : if (!is_general)
435 : {
436 : // count the entries on the diagonal
437 347 : int64_t nself_edges = 0 ;
438 347 : LG_TRY (LG_nself_edges (&nself_edges, A, msg)) ;
439 : // nvals_to_print = # of entries in tril(A), including diagonal
440 130 : nvals_to_print = nself_edges + (nvals - nself_edges) / 2 ;
441 : }
442 :
443 345 : if (MM_fmt == MM_array)
444 : {
445 : // write `nrows ncols` if the array format is used
446 19 : FPRINTF (f, "%" PRIu64 " %" PRIu64 "\n",
447 : nrows, ncols) ;
448 : }
449 : else
450 : {
451 : // otherwise write `nrows ncols nvals` for the coordinate format
452 326 : FPRINTF (f, "%" PRIu64 " %" PRIu64 " %" PRIu64 "\n",
453 : nrows, ncols, nvals_to_print) ;
454 : }
455 :
456 345 : if (nvals_to_print == 0)
457 : {
458 : // quick return if nothing more to do
459 2 : LG_FREE_ALL ;
460 2 : return (GrB_SUCCESS) ;
461 : }
462 :
463 : //--------------------------------------------------------------------------
464 : // extract and print tuples
465 : //--------------------------------------------------------------------------
466 :
467 343 : LG_TRY (LAGraph_Malloc ((void **) &I, nvals, sizeof (GrB_Index), msg)) ;
468 293 : LG_TRY (LAGraph_Malloc ((void **) &J, nvals, sizeof (GrB_Index), msg)) ;
469 243 : LG_TRY (LAGraph_Malloc ((void **) &K, nvals, sizeof (GrB_Index), msg)) ;
470 329549 : for (int64_t k = 0 ; k < nvals ; k++)
471 : {
472 329356 : K [k] = k ;
473 : }
474 :
475 193 : GrB_Index nvals_printed = 0 ;
476 193 : bool coord = (MM_fmt == MM_coordinate) ;
477 :
478 : #define WRITE_TUPLES(ctype,is_unsigned,is_signed,is_real,is_complex) \
479 : { \
480 : ctype *X = NULL ; \
481 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (ctype), msg)) ;\
482 : GRB_TRY (GrB_Matrix_extractTuples (I, J, X, &nvals, A)) ; \
483 : LG_TRY (LG_msort3 ((int64_t *) J, (int64_t *) I, \
484 : (int64_t *) K, nvals, msg)) ; \
485 : for (int64_t k = 0 ; k < nvals ; k++) \
486 : { \
487 : /* convert the row and column index to 1-based */ \
488 : GrB_Index i = I [k] + 1 ; \
489 : GrB_Index j = J [k] + 1 ; \
490 : ctype x = X [K [k]] ; \
491 : if (is_general || i >= j) \
492 : { \
493 : /* print the row and column index of the tuple */ \
494 : if (coord) FPRINTF (f, "%" PRIu64 " %" PRIu64 " ", i, j) ; \
495 : /* print the value of the tuple */ \
496 : if (is_structural) \
497 : { \
498 : /* print nothing */ ; \
499 : } \
500 : else if (is_unsigned) \
501 : { \
502 : FPRINTF (f, "%" PRIu64, (uint64_t) x) ; \
503 : } \
504 : else if (is_signed) \
505 : { \
506 : FPRINTF (f, "%" PRId64, (int64_t) x) ; \
507 : } \
508 : else if (is_real) \
509 : { \
510 : LG_ASSERT_MSG (print_double (f, (double) x), \
511 : LAGRAPH_IO_ERROR, "Unable to write to file") ; \
512 : } \
513 : /* else if (is_complex) */ \
514 : /* { */ \
515 : /* LG_ASSERT_MSG (print_double (f, creal (x)), */ \
516 : /* LAGRAPH_IO_ERROR, */ \
517 : /* "Unable to write to file") ; */ \
518 : /* FPRINTF (f, " ") ; */ \
519 : /* LG_ASSERT_MSG (print_double (f, cimag (x)), */ \
520 : /* LAGRAPH_IO_ERROR, */ \
521 : /* "Unable to write to file") ; */ \
522 : /* } */ \
523 : FPRINTF (f, "\n") ; \
524 : } \
525 : nvals_printed++ ; \
526 : } \
527 : LG_TRY (LAGraph_Free ((void **) &X, NULL)) ; \
528 : }
529 :
530 16069 : if (type == GrB_BOOL ) WRITE_TUPLES (bool , 1, 0, 0, 0)
531 200 : else if (type == GrB_INT8 ) WRITE_TUPLES (int8_t , 0, 1, 0, 0)
532 192 : else if (type == GrB_INT16 ) WRITE_TUPLES (int16_t , 0, 1, 0, 0)
533 244 : else if (type == GrB_INT32 ) WRITE_TUPLES (int32_t , 0, 1, 0, 0)
534 438 : else if (type == GrB_INT64 ) WRITE_TUPLES (int64_t , 0, 1, 0, 0)
535 95 : else if (type == GrB_UINT8 ) WRITE_TUPLES (uint8_t , 1, 0, 0, 0)
536 91 : else if (type == GrB_UINT16 ) WRITE_TUPLES (uint16_t, 1, 0, 0, 0)
537 87 : else if (type == GrB_UINT32 ) WRITE_TUPLES (uint32_t, 1, 0, 0, 0)
538 83 : else if (type == GrB_UINT64 ) WRITE_TUPLES (uint64_t, 1, 0, 0, 0)
539 119 : else if (type == GrB_FP32 ) WRITE_TUPLES (float , 0, 0, 1, 0)
540 202769 : else if (type == GrB_FP64 ) WRITE_TUPLES (double , 0, 0, 1, 0)
541 : #if 0
542 : else if (type == GxB_FC32 ) WRITE_TUPLES (GxB_FC32_t, 0, 0, 0, 1) ;
543 : else if (type == GxB_FC64 ) WRITE_TUPLES (GxB_FC64_t, 0, 0, 0, 1) ;
544 : #endif
545 :
546 : ASSERT (nvals_to_print == nvals_printed) ;
547 :
548 : //--------------------------------------------------------------------------
549 : // free workspace and return
550 : //--------------------------------------------------------------------------
551 :
552 105 : LG_FREE_ALL ;
553 105 : return (GrB_SUCCESS) ;
554 : }
|
