Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_MMRead: read a matrix from 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_MMRead: read a matrix from a Matrix Market file
19 :
20 : // Parts of this code are from SuiteSparse/CHOLMOD/Check/cholmod_read.c, and
21 : // are used here by permission of the author of CHOLMOD/Check (T. A. Davis).
22 :
23 : // The Matrix Market format is described at:
24 : // https://math.nist.gov/MatrixMarket/formats.html
25 :
26 : // Return values:
27 : // GrB_SUCCESS: input file and output matrix are valid
28 : // LAGRAPH_IO_ERROR: the input file cannot be read or has invalid content
29 : // GrB_NULL_POINTER: A or f are NULL on input
30 : // GrB_NOT_IMPLEMENTED: complex types not yet supported
31 : // other: return values directly from GrB_* methods
32 :
33 : #define LG_FREE_WORK \
34 : { \
35 : LAGraph_Free ((void **) &I, NULL) ; \
36 : LAGraph_Free ((void **) &J, NULL) ; \
37 : LAGraph_Free ((void **) &X, NULL) ; \
38 : }
39 :
40 : #define LG_FREE_ALL \
41 : { \
42 : LG_FREE_WORK ; \
43 : GrB_free (A) ; \
44 : }
45 :
46 : #include "LG_internal.h"
47 :
48 : //------------------------------------------------------------------------------
49 : // get_line
50 : //------------------------------------------------------------------------------
51 :
52 : // Read one line of the file, return true if successful, false if EOF.
53 : // The string is returned in buf, converted to lower case.
54 :
55 4065091 : static inline bool get_line
56 : (
57 : FILE *f, // file open for reading
58 : char *buf // size MAXLINE+1
59 : )
60 : {
61 :
62 : // check inputs
63 : ASSERT (f != NULL) ;
64 : ASSERT (buf != NULL) ;
65 :
66 : // read the line from the file
67 4065091 : buf [0] = '\0' ;
68 4065091 : buf [1] = '\0' ;
69 4065091 : if (fgets (buf, MAXLINE, f) == NULL)
70 : {
71 : // EOF or other I/O error
72 1 : return (false) ;
73 : }
74 4065090 : buf [MAXLINE] = '\0' ;
75 :
76 : // convert the string to lower case
77 75762915 : for (int k = 0 ; k < MAXLINE && buf [k] != '\0' ; k++)
78 : {
79 71697825 : buf [k] = tolower (buf [k]) ;
80 : }
81 4065090 : return (true) ;
82 : }
83 :
84 : //------------------------------------------------------------------------------
85 : // is_blank_line
86 : //------------------------------------------------------------------------------
87 :
88 : // returns true if buf is a blank line or comment, false otherwise.
89 :
90 4061578 : static inline bool is_blank_line
91 : (
92 : char *buf // size MAXLINE+1, never NULL
93 : )
94 : {
95 :
96 : // check inputs
97 : ASSERT (buf != NULL) ;
98 :
99 : // check if comment line
100 4061578 : if (buf [0] == '%')
101 : {
102 : // line is a comment
103 5256 : return (true) ;
104 : }
105 :
106 : // check if blank line
107 4056352 : for (int k = 0 ; k <= MAXLINE ; k++)
108 : {
109 4056352 : int c = buf [k] ;
110 4056352 : if (c == '\0')
111 : {
112 : // end of line
113 30 : break ;
114 : }
115 4056322 : if (!isspace (c))
116 : {
117 : // non-space character; this is not an error
118 4056292 : return (false) ;
119 : }
120 : }
121 :
122 : // line is blank
123 30 : return (true) ;
124 : }
125 :
126 : //------------------------------------------------------------------------------
127 : // read_double
128 : //------------------------------------------------------------------------------
129 :
130 : // Read a single double value from a string. The string may be any string
131 : // recognized by sscanf, or inf, -inf, +inf, or nan. The token infinity is
132 : // also OK instead of inf (only the first 3 letters of inf* or nan* are
133 : // significant, and the rest are ignored).
134 :
135 3157602 : static inline bool read_double // true if successful, false if failure
136 : (
137 : char *p, // string containing the value
138 : double *rval // value to read in
139 : )
140 : {
141 3157602 : while (*p && isspace (*p)) p++ ; // skip any spaces
142 :
143 3157602 : if (MATCH (p, "inf", 3) || MATCH (p, "+inf", 4))
144 : {
145 96 : (*rval) = INFINITY ;
146 : }
147 3157506 : else if (MATCH (p, "-inf", 4))
148 : {
149 46 : (*rval) = -INFINITY ;
150 : }
151 3157460 : else if (MATCH (p, "nan", 3))
152 : {
153 4 : (*rval) = NAN ;
154 : }
155 : else
156 : {
157 3157456 : if (sscanf (p, "%lg", rval) != 1)
158 : {
159 : // invalid file format, EOF, or other I/O error
160 1 : return (false) ;
161 : }
162 : }
163 3157601 : return (true) ;
164 : }
165 :
166 : //------------------------------------------------------------------------------
167 : // read_entry: read a numerical value and typecast to the given type
168 : //------------------------------------------------------------------------------
169 :
170 4054330 : static inline bool read_entry // returns true if successful, false if failure
171 : (
172 : char *p, // string containing the value
173 : GrB_Type type, // type of value to read
174 : bool structural, // if true, then the value is 1
175 : uint8_t *x // value read in, a pointer to space of size of the type
176 : )
177 : {
178 :
179 4054330 : int64_t ival = 1 ;
180 4054330 : double rval = 1, zval = 0 ;
181 :
182 4054330 : while (*p && isspace (*p)) p++ ; // skip any spaces
183 :
184 4054330 : if (type == GrB_BOOL)
185 : {
186 767170 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
187 767170 : if (ival < 0 || ival > 1)
188 : {
189 : // entry out of range
190 1 : return (false) ;
191 : }
192 767169 : bool *result = (bool *) x ;
193 767169 : result [0] = (bool) ival ;
194 : }
195 3287160 : else if (type == GrB_INT8)
196 : {
197 511 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
198 511 : if (ival < INT8_MIN || ival > INT8_MAX)
199 : {
200 : // entry out of range
201 1 : return (false) ;
202 : }
203 510 : int8_t *result = (int8_t *) x ;
204 510 : result [0] = (int8_t) ival ;
205 : }
206 3286649 : else if (type == GrB_INT16)
207 : {
208 523 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
209 523 : if (ival < INT16_MIN || ival > INT16_MAX)
210 : {
211 : // entry out of range
212 1 : return (false) ;
213 : }
214 522 : int16_t *result = (int16_t *) x ;
215 522 : result [0] = (int16_t) ival ;
216 : }
217 3286126 : else if (type == GrB_INT32)
218 : {
219 1603 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
220 1603 : if (ival < INT32_MIN || ival > INT32_MAX)
221 : {
222 : // entry out of range
223 1 : return (false) ;
224 : }
225 1602 : int32_t *result = (int32_t *) x ;
226 1602 : result [0] = (int32_t) ival ;
227 : }
228 3284523 : else if (type == GrB_INT64)
229 : {
230 3026 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
231 3026 : int64_t *result = (int64_t *) x ;
232 3026 : result [0] = (int64_t) ival ;
233 : }
234 3281497 : else if (type == GrB_UINT8)
235 : {
236 290 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
237 290 : if (ival < 0 || ival > UINT8_MAX)
238 : {
239 : // entry out of range
240 1 : return (false) ;
241 : }
242 289 : uint8_t *result = (uint8_t *) x ;
243 289 : result [0] = (uint8_t) ival ;
244 : }
245 3281207 : else if (type == GrB_UINT16)
246 : {
247 290 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
248 290 : if (ival < 0 || ival > UINT16_MAX)
249 : {
250 : // entry out of range
251 1 : return (false) ;
252 : }
253 289 : uint16_t *result = (uint16_t *) x ;
254 289 : result [0] = (uint16_t) ival ;
255 : }
256 3280917 : else if (type == GrB_UINT32)
257 : {
258 290 : if (!structural && sscanf (p, "%" SCNd64, &ival) != 1) return (false) ;
259 290 : if (ival < 0 || ival > UINT32_MAX)
260 : {
261 : // entry out of range
262 1 : return (false) ;
263 : }
264 289 : uint32_t *result = (uint32_t *) x ;
265 289 : result [0] = (uint32_t) ival ;
266 : }
267 3280627 : else if (type == GrB_UINT64)
268 : {
269 312 : uint64_t uval = 1 ;
270 312 : if (!structural && sscanf (p, "%" SCNu64, &uval) != 1) return (false) ;
271 312 : uint64_t *result = (uint64_t *) x ;
272 312 : result [0] = (uint64_t) uval ;
273 : }
274 3280315 : else if (type == GrB_FP32)
275 : {
276 1090459 : if (!structural && !read_double (p, &rval)) return (false) ;
277 1090459 : float *result = (float *) x ;
278 1090459 : result [0] = (float) rval ;
279 : }
280 2189856 : else if (type == GrB_FP64)
281 : {
282 2189856 : if (!structural && !read_double (p, &rval)) return (false) ;
283 2189855 : double *result = (double *) x ;
284 2189855 : result [0] = rval ;
285 : }
286 : #if 0
287 : else if (type == GxB_FC32)
288 : {
289 : if (!structural && !read_double (p, &rval)) return (false) ;
290 : while (*p && !isspace (*p)) p++ ; // skip real part
291 : if (!structural && !read_double (p, &zval)) return (false) ;
292 : float *result = (float *) x ;
293 : result [0] = (float) rval ; // real part
294 : result [1] = (float) zval ; // imaginary part
295 : }
296 : else if (type == GxB_FC64)
297 : {
298 : if (!structural && !read_double (p, &rval)) return (false) ;
299 : while (*p && !isspace (*p)) p++ ; // skip real part
300 : if (!structural && !read_double (p, &zval)) return (false) ;
301 : double *result = (double *) x ;
302 : result [0] = rval ; // real part
303 : result [1] = zval ; // imaginary part
304 : }
305 : #endif
306 :
307 4054322 : return (true) ;
308 : }
309 :
310 : //------------------------------------------------------------------------------
311 : // negate_scalar: negate a scalar value
312 : //------------------------------------------------------------------------------
313 :
314 : // negate the scalar x. Do nothing for bool or uint*.
315 :
316 1340 : static inline void negate_scalar
317 : (
318 : GrB_Type type,
319 : uint8_t *x
320 : )
321 : {
322 :
323 1340 : if (type == GrB_INT8)
324 : {
325 210 : int8_t *value = (int8_t *) x ;
326 210 : (*value) = - (*value) ;
327 : }
328 1130 : else if (type == GrB_INT16)
329 : {
330 210 : int16_t *value = (int16_t *) x ;
331 210 : (*value) = - (*value) ;
332 : }
333 920 : else if (type == GrB_INT32)
334 : {
335 210 : int32_t *value = (int32_t *) x ;
336 210 : (*value) = - (*value) ;
337 : }
338 710 : else if (type == GrB_INT64)
339 : {
340 210 : int64_t *value = (int64_t *) x ;
341 210 : (*value) = - (*value) ;
342 : }
343 500 : else if (type == GrB_FP32)
344 : {
345 290 : float *value = (float *) x ;
346 290 : (*value) = - (*value) ;
347 : }
348 210 : else if (type == GrB_FP64)
349 : {
350 210 : double *value = (double *) x ;
351 210 : (*value) = - (*value) ;
352 : }
353 : #if 0
354 : else if (type == GxB_FC32)
355 : {
356 : float complex *value = (float complex *) x ;
357 : (*value) = - (*value) ;
358 : }
359 : else if (type == GxB_FC64)
360 : {
361 : double complex *value = (double complex *) x ;
362 : (*value) = - (*value) ;
363 : }
364 : #endif
365 1340 : }
366 :
367 : //------------------------------------------------------------------------------
368 : // set_value
369 : //------------------------------------------------------------------------------
370 :
371 : // Add the (i,j,x) triplet to the I,J,X arrays as the kth triplet, and
372 : // increment k. No typecasting is done.
373 :
374 6291344 : static inline void set_value
375 : (
376 : size_t typesize, // size of the numerical type, in bytes
377 : GrB_Index i,
378 : GrB_Index j,
379 : uint8_t *x, // scalar, an array of size at least typesize
380 : GrB_Index *I,
381 : GrB_Index *J,
382 : uint8_t *X,
383 : GrB_Index *k // # of triplets
384 : )
385 : {
386 6291344 : I [*k] = i ;
387 6291344 : J [*k] = j ;
388 6291344 : memcpy (X + ((*k) * typesize), x, typesize) ;
389 6291344 : (*k)++ ;
390 6291344 : }
391 :
392 : //------------------------------------------------------------------------------
393 : // LAGraph_MMRead
394 : //------------------------------------------------------------------------------
395 :
396 1974 : int LAGraph_MMRead
397 : (
398 : // output:
399 : GrB_Matrix *A, // handle of matrix to create
400 : // input:
401 : FILE *f, // file to read from, already open
402 : char *msg
403 : )
404 : {
405 :
406 : //--------------------------------------------------------------------------
407 : // check inputs
408 : //--------------------------------------------------------------------------
409 :
410 1974 : GrB_Index *I = NULL, *J = NULL ;
411 1974 : uint8_t *X = NULL ;
412 1974 : LG_CLEAR_MSG ;
413 1974 : LG_ASSERT (A != NULL, GrB_NULL_POINTER) ;
414 1973 : LG_ASSERT (f != NULL, GrB_NULL_POINTER) ;
415 1972 : (*A) = NULL ;
416 :
417 : //--------------------------------------------------------------------------
418 : // set the default properties
419 : //--------------------------------------------------------------------------
420 :
421 1972 : MM_fmt_enum MM_fmt = MM_coordinate ;
422 1972 : MM_type_enum MM_type = MM_real ;
423 1972 : MM_storage_enum MM_storage = MM_general ;
424 1972 : GrB_Type type = GrB_FP64 ;
425 1972 : size_t typesize = sizeof (double) ;
426 1972 : GrB_Index nrows = 0 ;
427 1972 : GrB_Index ncols = 0 ;
428 1972 : GrB_Index nvals = 0 ;
429 :
430 : //--------------------------------------------------------------------------
431 : // read the Matrix Market header
432 : //--------------------------------------------------------------------------
433 :
434 : // Read the header. This consists of zero or more comment lines (blank, or
435 : // starting with a "%" in the first column), followed by a single data line
436 : // containing two or three numerical values. The first line is normally:
437 : //
438 : // %%MatrixMarket matrix <fmt> <type> <storage>
439 : //
440 : // but this is optional. The 2nd line is also optional (the %%MatrixMarket
441 : // line is required for this 2nd line to be recognized):
442 : //
443 : // %%GraphBLAS type <Ctype>
444 : //
445 : // where the Ctype is one of: bool, int8_t, int16_t, int32_t, int64_t,
446 : // uint8_t, uint16_t, uint32_t, uint64_t, float, or double.
447 : //
448 : // If the %%MatrixMarket line is not present, then the <fmt> <type> and
449 : // <storage> are implicit. If the first data line contains 3 items,
450 : // then the implicit header is:
451 : //
452 : // %%MatrixMarket matrix coordinate real general
453 : // %%GraphBLAS type double
454 : //
455 : // If the first data line contains 2 items (nrows ncols), then the implicit
456 : // header is:
457 : //
458 : // %%MatrixMarket matrix array real general
459 : // %%GraphBLAS type double
460 : //
461 : // The implicit header is an extension of the Matrix Market format.
462 :
463 : char buf [MAXLINE+1] ;
464 :
465 1972 : bool got_mm_header = false ;
466 1972 : bool got_first_data_line = false ;
467 : int64_t line ;
468 :
469 10713 : for (line = 1 ; get_line (f, buf) ; line++)
470 : {
471 :
472 : //----------------------------------------------------------------------
473 : // parse the line
474 : //----------------------------------------------------------------------
475 :
476 10713 : if ((line == 1) && MATCH (buf, "%%matrixmarket", 14))
477 1956 : {
478 :
479 : //------------------------------------------------------------------
480 : // read a Matrix Market header
481 : //------------------------------------------------------------------
482 :
483 : // %%MatrixMarket matrix <fmt> <type> <storage>
484 : // if present, it must be the first line in the file.
485 :
486 1964 : got_mm_header = true ;
487 1964 : char *p = buf + 14 ;
488 :
489 : //------------------------------------------------------------------
490 : // get "matrix" token and discard it
491 : //------------------------------------------------------------------
492 :
493 3928 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
494 :
495 1964 : if (!MATCH (p, "matrix", 6))
496 : {
497 : // invalid Matrix Market object
498 1 : LG_ASSERT_MSG (false,
499 : LAGRAPH_IO_ERROR, "invalid MatrixMarket header"
500 : " ('matrix' token missing)") ;
501 : }
502 1963 : p += 6 ; // skip past token "matrix"
503 :
504 : //------------------------------------------------------------------
505 : // get the fmt token
506 : //------------------------------------------------------------------
507 :
508 3926 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
509 :
510 1963 : if (MATCH (p, "coordinate", 10))
511 : {
512 1876 : MM_fmt = MM_coordinate ;
513 1876 : p += 10 ;
514 : }
515 87 : else if (MATCH (p, "array", 5))
516 : {
517 86 : MM_fmt = MM_array ;
518 86 : p += 5 ;
519 : }
520 : else
521 : {
522 : // invalid Matrix Market format
523 1 : LG_ASSERT_MSG (false,
524 : LAGRAPH_IO_ERROR, "invalid format in MatrixMarket header"
525 : " (format must be 'coordinate' or 'array')") ;
526 : }
527 :
528 : //------------------------------------------------------------------
529 : // get the Matrix Market type token
530 : //------------------------------------------------------------------
531 :
532 3924 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
533 :
534 1962 : if (MATCH (p, "real", 4))
535 : {
536 621 : MM_type = MM_real ;
537 621 : type = GrB_FP64 ;
538 621 : typesize = sizeof (double) ;
539 621 : p += 4 ;
540 : }
541 1341 : else if (MATCH (p, "integer", 7))
542 : {
543 799 : MM_type = MM_integer ;
544 799 : type = GrB_INT64 ;
545 799 : typesize = sizeof (int64_t) ;
546 799 : p += 7 ;
547 : }
548 542 : else if (MATCH (p, "complex", 7))
549 : {
550 1 : MM_type = MM_complex ;
551 : #if 0
552 : type = GxB_FC64 ;
553 : typesize = sizeof (GxB_FC64_t) ;
554 : p += 7 ;
555 : #endif
556 1 : LG_ASSERT_MSG (false,
557 : GrB_NOT_IMPLEMENTED, "complex types not supported") ;
558 : }
559 541 : else if (MATCH (p, "pattern", 7))
560 : {
561 540 : MM_type = MM_pattern ;
562 540 : type = GrB_BOOL ;
563 540 : typesize = sizeof (bool) ;
564 540 : p += 7 ;
565 : }
566 : else
567 : {
568 : // invalid Matrix Market type
569 1 : LG_ASSERT_MSG (false,
570 : LAGRAPH_IO_ERROR, "invalid MatrixMarket type") ;
571 : }
572 :
573 : //------------------------------------------------------------------
574 : // get the storage token
575 : //------------------------------------------------------------------
576 :
577 3920 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
578 :
579 1960 : if (MATCH (p, "general", 7))
580 : {
581 1170 : MM_storage = MM_general ;
582 : }
583 790 : else if (MATCH (p, "symmetric", 9))
584 : {
585 616 : MM_storage = MM_symmetric ;
586 : }
587 174 : else if (MATCH (p, "skew-symmetric", 14))
588 : {
589 171 : MM_storage = MM_skew_symmetric ;
590 : }
591 3 : else if (MATCH (p, "hermitian", 9))
592 : {
593 2 : MM_storage = MM_hermitian ;
594 : }
595 : else
596 : {
597 : // invalid Matrix Market storage
598 1 : LG_ASSERT_MSG (false,
599 : LAGRAPH_IO_ERROR, "invalid MatrixMarket storage") ;
600 : }
601 :
602 : //------------------------------------------------------------------
603 : // ensure the combinations are valid
604 : //------------------------------------------------------------------
605 :
606 1959 : if (MM_type == MM_pattern)
607 : {
608 : // (coodinate) x (pattern) x (general or symmetric)
609 540 : LG_ASSERT_MSG (
610 : (MM_fmt == MM_coordinate &&
611 : (MM_storage == MM_general || MM_storage == MM_symmetric)),
612 : LAGRAPH_IO_ERROR,
613 : "invalid MatrixMarket pattern combination") ;
614 : }
615 :
616 1957 : if (MM_storage == MM_hermitian)
617 : {
618 : // (coordinate or array) x (complex) x (Hermitian)
619 1 : LG_ASSERT_MSG (MM_type == MM_complex,
620 : LAGRAPH_IO_ERROR,
621 : "invalid MatrixMarket complex combination") ;
622 : }
623 :
624 : }
625 8749 : else if (got_mm_header && MATCH (buf, "%%graphblas", 11))
626 1498 : {
627 :
628 : //------------------------------------------------------------------
629 : // %%GraphBLAS structured comment
630 : //------------------------------------------------------------------
631 :
632 1548 : char *p = buf + 11 ;
633 3096 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
634 :
635 1548 : if (MATCH (p, "type", 4) && !got_first_data_line)
636 : {
637 :
638 : //--------------------------------------------------------------
639 : // %%GraphBLAS type <Ctype>
640 : //--------------------------------------------------------------
641 :
642 : // This must appear after the %%MatrixMarket header and before
643 : // the first data line. Otherwise the %%GraphBLAS line is
644 : // treated as a pure comment.
645 :
646 1502 : p += 4 ;
647 3004 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
648 :
649 : // Ctype is one of: bool, int8_t, int16_t, int32_t, int64_t,
650 : // uint8_t, uint16_t, uint32_t, uint64_t, float, or double.
651 : // The complex types "float complex", or "double complex" are
652 : // not yet supported.
653 :
654 1502 : if (MATCH (p, "bool", 4))
655 : {
656 519 : type = GrB_BOOL ;
657 519 : typesize = sizeof (bool) ;
658 : }
659 983 : else if (MATCH (p, "int8_t", 6))
660 : {
661 59 : type = GrB_INT8 ;
662 59 : typesize = sizeof (int8_t) ;
663 : }
664 924 : else if (MATCH (p, "int16_t", 7))
665 : {
666 60 : type = GrB_INT16 ;
667 60 : typesize = sizeof (int16_t) ;
668 : }
669 864 : else if (MATCH (p, "int32_t", 7))
670 : {
671 176 : type = GrB_INT32 ;
672 176 : typesize = sizeof (int32_t) ;
673 : }
674 688 : else if (MATCH (p, "int64_t", 7))
675 : {
676 160 : type = GrB_INT64 ;
677 160 : typesize = sizeof (int64_t) ;
678 : }
679 528 : else if (MATCH (p, "uint8_t", 7))
680 : {
681 32 : type = GrB_UINT8 ;
682 32 : typesize = sizeof (uint8_t) ;
683 : }
684 496 : else if (MATCH (p, "uint16_t", 8))
685 : {
686 31 : type = GrB_UINT16 ;
687 31 : typesize = sizeof (uint16_t) ;
688 : }
689 465 : else if (MATCH (p, "uint32_t", 8))
690 : {
691 31 : type = GrB_UINT32 ;
692 31 : typesize = sizeof (uint32_t) ;
693 : }
694 434 : else if (MATCH (p, "uint64_t", 8))
695 : {
696 32 : type = GrB_UINT64 ;
697 32 : typesize = sizeof (uint64_t) ;
698 : }
699 402 : else if (MATCH (p, "float complex", 13))
700 : {
701 : #if 0
702 : type = GxB_FC32 ;
703 : typesize = sizeof (GxB_FC32_t) ;
704 : #endif
705 1 : LG_ASSERT_MSG (false,
706 : GrB_NOT_IMPLEMENTED, "complex types not supported") ;
707 : }
708 401 : else if (MATCH (p, "double complex", 14))
709 : {
710 : #if 0
711 : type = GxB_FC64 ;
712 : typesize = sizeof (GxB_FC64_t) ;
713 : #endif
714 1 : LG_ASSERT_MSG (false,
715 : GrB_NOT_IMPLEMENTED, "complex types not supported") ;
716 : }
717 400 : else if (MATCH (p, "float", 5))
718 : {
719 96 : type = GrB_FP32 ;
720 96 : typesize = sizeof (float) ;
721 : }
722 304 : else if (MATCH (p, "double", 6))
723 : {
724 303 : type = GrB_FP64 ;
725 303 : typesize = sizeof (double) ;
726 : }
727 : else
728 : {
729 : // unknown type
730 1 : LG_ASSERT_MSG (false,
731 : LAGRAPH_IO_ERROR, "unknown type") ;
732 : }
733 :
734 1499 : if (MM_storage == MM_skew_symmetric && (type == GrB_BOOL ||
735 160 : type == GrB_UINT8 || type == GrB_UINT16 ||
736 159 : type == GrB_UINT32 || type == GrB_UINT64))
737 : {
738 : // matrices with unsigned types cannot be skew-symmetric
739 1 : LG_ASSERT_MSG (false, LAGRAPH_IO_ERROR,
740 : "skew-symmetric matrices cannot have an unsigned type");
741 : }
742 : }
743 : else
744 : {
745 : // %%GraphBLAS line but no "type" as the 2nd token; ignore it
746 46 : continue ;
747 : }
748 :
749 : }
750 7201 : else if (is_blank_line (buf))
751 : {
752 :
753 : // -----------------------------------------------------------------
754 : // blank line or comment line
755 : // -----------------------------------------------------------------
756 :
757 5241 : continue ;
758 :
759 : }
760 : else
761 : {
762 :
763 : // -----------------------------------------------------------------
764 : // read the first data line
765 : // -----------------------------------------------------------------
766 :
767 : // format: [nrows ncols nvals] or just [nrows ncols]
768 :
769 1960 : got_first_data_line = true ;
770 1960 : int nitems = sscanf (buf, "%" SCNu64 " %" SCNu64 " %" SCNu64,
771 : &nrows, &ncols, &nvals) ;
772 :
773 1960 : if (nitems == 2)
774 : {
775 : // a dense matrix
776 43 : if (!got_mm_header)
777 : {
778 : // if no header, treat it as if it were
779 : // %%MatrixMarket matrix array real general
780 4 : MM_fmt = MM_array ;
781 4 : MM_type = MM_real ;
782 4 : MM_storage = MM_general ;
783 4 : type = GrB_FP64 ;
784 4 : typesize = sizeof (double) ;
785 : }
786 43 : if (MM_storage == MM_general)
787 : {
788 : // dense general matrix
789 25 : nvals = nrows * ncols ;
790 : }
791 : else
792 : {
793 : // dense symmetric, skew-symmetric, or hermitian matrix
794 18 : nvals = nrows + ((nrows * nrows - nrows) / 2) ;
795 : }
796 : }
797 1917 : else if (nitems == 3)
798 : {
799 : // a sparse matrix
800 1916 : if (!got_mm_header)
801 : {
802 : // if no header, treat it as if it were
803 : // %%MatrixMarket matrix coordinate real general
804 4 : MM_fmt = MM_coordinate ;
805 4 : MM_type = MM_real ;
806 4 : MM_storage = MM_general ;
807 4 : type = GrB_FP64 ;
808 4 : typesize = sizeof (double) ;
809 : }
810 : }
811 : else
812 : {
813 : // wrong number of items in first data line
814 1 : LG_ASSERT_MSGF (false,
815 : LAGRAPH_IO_ERROR, "invalid 1st data line"
816 : " (line %" PRId64 " of input file)", line) ;
817 : }
818 :
819 1959 : if (nrows != ncols)
820 : {
821 : // a rectangular matrix must be in the general storage
822 90 : LG_ASSERT_MSG (MM_storage == MM_general,
823 : LAGRAPH_IO_ERROR, "invalid rectangular storage") ;
824 : }
825 :
826 : //------------------------------------------------------------------
827 : // header has been read in
828 : //------------------------------------------------------------------
829 :
830 1958 : break ;
831 : }
832 : }
833 :
834 : //--------------------------------------------------------------------------
835 : // create the matrix
836 : //--------------------------------------------------------------------------
837 :
838 1958 : GRB_TRY (GrB_Matrix_new (A, type, nrows, ncols)) ;
839 :
840 : //--------------------------------------------------------------------------
841 : // quick return for empty matrix
842 : //--------------------------------------------------------------------------
843 :
844 1806 : if (nrows == 0 || ncols == 0 || nvals == 0)
845 : {
846 : // success: return an empty matrix. This is not an error.
847 11 : return (GrB_SUCCESS) ;
848 : }
849 :
850 : //--------------------------------------------------------------------------
851 : // allocate space for the triplets
852 : //--------------------------------------------------------------------------
853 :
854 1795 : GrB_Index nvals3 = ((MM_storage == MM_general) ? 1 : 2) * (nvals + 1) ;
855 1795 : LG_TRY (LAGraph_Malloc ((void **) &I, nvals3, sizeof (GrB_Index), msg)) ;
856 1745 : LG_TRY (LAGraph_Malloc ((void **) &J, nvals3, sizeof (GrB_Index), msg)) ;
857 1695 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals3, typesize, msg)) ;
858 :
859 : //--------------------------------------------------------------------------
860 : // read in the triplets
861 : //--------------------------------------------------------------------------
862 :
863 1645 : GrB_Index i = -1, j = 0 ;
864 1645 : GrB_Index nvals2 = 0 ;
865 4055967 : for (int64_t k = 0 ; k < nvals ; k++)
866 : {
867 :
868 : //----------------------------------------------------------------------
869 : // get the next triplet, skipping blank lines and comment lines
870 : //----------------------------------------------------------------------
871 :
872 : uint8_t x [MAXLINE] ; // scalar value
873 :
874 : while (true)
875 45 : {
876 :
877 : //------------------------------------------------------------------
878 : // read the file until finding the next triplet
879 : //------------------------------------------------------------------
880 :
881 4054378 : bool ok = get_line (f, buf) ;
882 4054378 : line++ ;
883 4054388 : LG_ASSERT_MSG (ok, LAGRAPH_IO_ERROR, "premature EOF") ;
884 4054377 : if (is_blank_line (buf))
885 : {
886 : // blank line or comment
887 45 : continue ;
888 : }
889 :
890 : //------------------------------------------------------------------
891 : // get the row and column index
892 : //------------------------------------------------------------------
893 :
894 4054332 : char *p = buf ;
895 4054332 : if (MM_fmt == MM_array)
896 : {
897 : // array format, column major order
898 977 : i++ ;
899 977 : if (i == nrows)
900 : {
901 152 : j++ ;
902 152 : if (MM_storage == MM_general)
903 : {
904 : // dense matrix in column major order
905 74 : i = 0 ;
906 : }
907 : else
908 : {
909 : // dense matrix in column major order, only the lower
910 : // triangular form is present, including the diagonal
911 78 : i = j ;
912 : }
913 : }
914 : }
915 : else
916 : {
917 : // coordinate format; read the row index and column index
918 4053355 : int inputs = sscanf (p, "%" SCNu64 " %" SCNu64, &i, &j) ;
919 4053355 : LG_ASSERT_MSGF (inputs == 2, LAGRAPH_IO_ERROR,
920 : "line %" PRId64 " of input file: indices invalid", line) ;
921 : // check the indices (they are 1-based in the MM file format)
922 4053354 : LG_ASSERT_MSGF (i >= 1 && i <= nrows, GrB_INDEX_OUT_OF_BOUNDS,
923 : "line %" PRId64 " of input file: row index %" PRIu64
924 : " out of range (must be in range 1 to %" PRIu64")",
925 : line, i, nrows) ;
926 4053354 : LG_ASSERT_MSGF (j >= 1 && j <= ncols, GrB_INDEX_OUT_OF_BOUNDS,
927 : "line %" PRId64 " of input file: column index %" PRIu64
928 : " out of range (must be in range 1 to %" PRIu64")",
929 : line, j, ncols) ;
930 : // convert from 1-based to 0-based.
931 4053353 : i-- ;
932 4053353 : j-- ;
933 : // advance p to the 3rd token to get the value of the entry
934 4053353 : while (*p && isspace (*p)) p++ ; // skip any leading spaces
935 17619616 : while (*p && !isspace (*p)) p++ ; // skip the row index
936 8106706 : while (*p && isspace (*p)) p++ ; // skip any spaces
937 17350928 : while (*p && !isspace (*p)) p++ ; // skip the column index
938 : }
939 :
940 : //------------------------------------------------------------------
941 : // read the value of the entry
942 : //------------------------------------------------------------------
943 :
944 8155058 : while (*p && isspace (*p)) p++ ; // skip any spaces
945 :
946 4054330 : ok = read_entry (p, type, MM_type == MM_pattern, x) ;
947 4054330 : LG_ASSERT_MSGF (ok, LAGRAPH_IO_ERROR, "entry value invalid on line"
948 : " %" PRId64 " of input file", line) ;
949 :
950 : //------------------------------------------------------------------
951 : // set the value in the matrix
952 : //------------------------------------------------------------------
953 :
954 4054322 : set_value (typesize, i, j, x, I, J, X, &nvals2) ;
955 :
956 : //------------------------------------------------------------------
957 : // also set the A(j,i) entry, if symmetric
958 : //------------------------------------------------------------------
959 :
960 4054322 : if (i != j && MM_storage != MM_general)
961 : {
962 2237022 : if (MM_storage == MM_symmetric)
963 : {
964 2235682 : set_value (typesize, j, i, x, I, J, X, &nvals2) ;
965 : }
966 1340 : else if (MM_storage == MM_skew_symmetric)
967 : {
968 1340 : negate_scalar (type, x) ;
969 1340 : set_value (typesize, j, i, x, I, J, X, &nvals2) ;
970 : }
971 : #if 0
972 : else if (MM_storage == MM_hermitian)
973 : {
974 : double complex *value = (double complex *) x ;
975 : (*value) = conj (*value) ;
976 : set_value (typesize, j, i, x, I, J, X, &nvals2) ;
977 : }
978 : #endif
979 : }
980 :
981 : // one more entry has been read in
982 4054322 : break ;
983 : }
984 : }
985 :
986 : //--------------------------------------------------------------------------
987 : // build the final matrix
988 : //--------------------------------------------------------------------------
989 :
990 1634 : if (type == GrB_BOOL)
991 : {
992 528 : GRB_TRY (GrB_Matrix_build_BOOL (*A, I, J, (bool *) X, nvals2, NULL)) ;
993 : }
994 1106 : else if (type == GrB_INT8)
995 : {
996 46 : GRB_TRY (GrB_Matrix_build_INT8 (*A, I, J, (int8_t *) X, nvals2, NULL)) ;
997 : }
998 1060 : else if (type == GrB_INT16)
999 : {
1000 47 : GRB_TRY (GrB_Matrix_build_INT16 (*A, I, J, (int16_t *) X, nvals2, NULL)) ;
1001 : }
1002 1013 : else if (type == GrB_INT32)
1003 : {
1004 137 : GRB_TRY (GrB_Matrix_build_INT32 (*A, I, J, (int32_t *) X, nvals2, NULL)) ;
1005 : }
1006 876 : else if (type == GrB_INT64)
1007 : {
1008 232 : GRB_TRY (GrB_Matrix_build_INT64 (*A, I, J, (int64_t *) X, nvals2, NULL)) ;
1009 : }
1010 644 : else if (type == GrB_UINT8)
1011 : {
1012 24 : GRB_TRY (GrB_Matrix_build_UINT8 (*A, I, J, (uint8_t *) X, nvals2, NULL)) ;
1013 : }
1014 620 : else if (type == GrB_UINT16)
1015 : {
1016 24 : GRB_TRY (GrB_Matrix_build_UINT16 (*A, I, J, (uint16_t *) X, nvals2, NULL)) ;
1017 : }
1018 596 : else if (type == GrB_UINT32)
1019 : {
1020 24 : GRB_TRY (GrB_Matrix_build_UINT32 (*A, I, J, (uint32_t *) X, nvals2, NULL)) ;
1021 : }
1022 572 : else if (type == GrB_UINT64)
1023 : {
1024 26 : GRB_TRY (GrB_Matrix_build_UINT64 (*A, I, J, (uint64_t *) X, nvals2, NULL)) ;
1025 : }
1026 546 : else if (type == GrB_FP32)
1027 : {
1028 84 : GRB_TRY (GrB_Matrix_build_FP32 (*A, I, J, (float *) X, nvals2, NULL)) ;
1029 : }
1030 462 : else if (type == GrB_FP64)
1031 : {
1032 462 : GRB_TRY (GrB_Matrix_build_FP64 (*A, I, J, (double *) X, nvals2, NULL)) ;
1033 : }
1034 : #if 0
1035 : else if (type == GxB_FC32)
1036 : {
1037 : GRB_TRY (GxB_Matrix_build_FC32 (*A, I, J, (GxB_FC32_t *) X, nvals2, NULL)) ;
1038 : }
1039 : else if (type == GxB_FC64)
1040 : {
1041 : GRB_TRY (GxB_Matrix_build_FC64 (*A, I, J, (GxB_FC64_t *) X, nvals2, NULL)) ;
1042 : }
1043 : #endif
1044 :
1045 : //--------------------------------------------------------------------------
1046 : // free workspace and return result
1047 : //--------------------------------------------------------------------------
1048 :
1049 1247 : LG_FREE_WORK ;
1050 1247 : return (GrB_SUCCESS) ;
1051 : }
|
