1 /* ========================================================================== */
2 /* === RBio/Source/RBio.c: C-callable RBio functions ======================== */
3 /* ========================================================================== */
4
5 /* Copyright 2009, Timothy A. Davis, All Rights Reserved.
6 Refer to RBio/Doc/license.txt for the RBio license. */
7
8 /* This file contains functions for writing/reading a sparse matrix to/from a
9 file in Rutherford-Boeing format. User-callable functions are declared
10 as PUBLIC. PRIVATE functions are available only within this file.
11 */
12
13 /* ========================================================================== */
14 /* === definitions ========================================================== */
15 /* ========================================================================== */
16
17 #include "RBio.h"
18
19 #ifdef INT
20 /* int version */
21 #define Int int
22 #define IDD "d"
23 #define RB(name) RB ## name ## _i
24 #else
25 /* Default: long (except for Windows, which is __int64) */
26 #define Int SuiteSparse_long
27 #define IDD SuiteSparse_long_idd
28 #define RB(name) RB ## name
29 #endif
30 #define ID "%" IDD
31
32 #define TRUE (1)
33 #define FALSE (0)
34 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
35 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
36 #define ABS(a) (((a) > 0) ? (a) : -(a))
37 #define ISNAN(a) ((a) != (a))
38 #define PRIVATE static
39 #define PUBLIC
40
41 #define SLEN 4096
42 #define FREE_WORK { SuiteSparse_free (w) ; \
43 SuiteSparse_free (cp) ; }
44
45 #define FREE_ALL { FREE_WORK ; \
46 SuiteSparse_free (Ap) ; \
47 SuiteSparse_free (Ai) ; \
48 SuiteSparse_free (Ax) ; \
49 SuiteSparse_free (Az) ; \
50 SuiteSparse_free (Zp) ; \
51 SuiteSparse_free (Zi) ; }
52
53 #define FREE_RAW { SuiteSparse_free (Ap) ; \
54 SuiteSparse_free (Ai) ; \
55 SuiteSparse_free (Ax) ; }
56
57
58 /* ========================================================================== */
59 /* === internal prototypes ================================================== */
60 /* ========================================================================== */
61
62 PRIVATE Int RB(format) /* return format to use (index in F_, C_format) */
63 (
64 /* input */
65 Int nnz, /* number of nonzeros */
66 double *x, /* of size nnz */
67 Int is_int, /* true if integer format is to be used */
68 double xmin, /* minimum value of x */
69 double xmax, /* maximum value of x */
70 Int fmt, /* initial format to use (index into F_format, ...) */
71
72 /* output */
73 char valfmt [21], /* Fortran format to use */
74 char valcfm [21], /* C format to use */
75 Int *valn /* number of entries per line */
76 ) ;
77
78 PRIVATE void RB(iformat)
79 (
80 /* input */
81 double xmin, /* smallest integer to print */
82 double xmax, /* largest integer to print */
83
84 /* output */
85 char indfmt [21], /* Fortran format to use */
86 char indcfm [21], /* C format to use */
87 Int *indn /* number of entries per line */
88 ) ;
89
90 PRIVATE Int RB(cards)
91 (
92 Int nitems, /* number of items to print */
93 Int nperline /* number of items per line */
94 ) ;
95
96 PRIVATE Int RB(iprint) /* returns TRUE if OK, FALSE otherwise */
97 (
98 /* input */
99 FILE *file, /* which file to write to */
100 char *indcfm, /* C format to use */
101 Int i, /* value to write */
102 Int indn, /* number of entries to write per line */
103
104 /* input/output */
105 Int *nbuf /* number of entries written to current line */
106 ) ;
107
108 PRIVATE Int RB(xprint) /* returns TRUE if OK, FALSE otherwise */
109 (
110 /* input */
111 FILE *file, /* which file to write to */
112 char *valcfm, /* C format to use */
113 double x, /* value to write */
114 Int valn, /* number of entries to write per line */
115 Int mkind, /* 0:real, 1:pattern, 2:complex, 3:integer */
116
117 /* input/output */
118 Int *nbuf /* number of entries written to current line */
119 ) ;
120
121 PRIVATE void RB(fill)
122 (
123 char *s, /* string to fill */
124 Int len, /* length of s (including trailing '\0') */
125 char c /* character to fill s with */
126 ) ;
127
128 PRIVATE Int RB(fix_mkind_in) /* return revised mkind */
129 (
130 Int mkind_in, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
131 double *Ax,
132 double *Az
133 ) ;
134
135 PRIVATE Int RB(writeTask) /* returns TRUE if OK, FALSE on failure */
136 (
137 /* input */
138 FILE *file, /* file to print to (already open) */
139 Int task, /* 0 to 3 (see above) */
140 Int nrow, /* A is nrow-by-ncol */
141 Int ncol,
142 Int mkind, /* 0:real, 1:pattern, 2:complex, 3:integer */
143 Int skind, /* -1:rect, 0:unsym, 1:sym, 2:hermitian, 3:skew */
144 Int *Ap, /* size ncol+1, column pointers */
145 Int *Ai, /* size anz=Ap[ncol], row indices */
146 double *Ax, /* size anz, real values */
147 double *Az, /* size anz, imaginary part (may be NULL) */
148 Int *Zp, /* size ncol+1, column pointers for Z (may be NULL) */
149 Int *Zi, /* size Zp[ncol], row indices for Z */
150 char *indcfm, /* C format for indices */
151 Int indn, /* # of indices per line */
152 char *valcfm, /* C format for values */
153 Int valn, /* # of values per line */
154
155 /* output */
156 Int *nnz, /* number of entries that will be printed to the file */
157
158 /* workspace */
159 Int *w, /* size MAX(nrow,ncol)+1 */
160 Int *cp /* size MAX(nrow,ncol)+1 */
161 ) ;
162
163 PRIVATE Int RB(read2) /* 0: OK, < 0: error, > 0: warning */
164 (
165 /* input */
166 FILE *file, /* must be already open for reading */
167 Int nrow, /* A is nrow-by-ncol */
168 Int ncol,
169 Int nnz, /* number of entries in A, from the header */
170 Int mkind, /* R: 0, P: 1, C: 2, I: 3 */
171 Int skind, /* R: -1, U: 0, S: 1, H: 2, Z: 3 */
172 Int build_upper, /* if TRUE and skind>0, then create upper part. Ai, Ax,
173 and Az must be twice the size as given below */
174
175 /* output */
176 Int *Ap, /* size ncol+1, column pointers for A */
177 Int *Ai, /* size nnz, row indices for A */
178 double *Ax, /* size nnz or 2*nnz if complex and Az NULL */
179 double *Az, /* size nnz, or NULL, for complex matrices only */
180
181 /* workspace */
182 Int *w, /* size MAX(nrow,ncol)+1 */
183 Int *cp, /* size ncol+1 */
184
185 /* input/workspace */
186 char *s, /* first line of column pointers on input */
187 Int slen
188 ) ;
189
190 PRIVATE Int RB(zcount) /* return number of explicit zeros in A */
191 (
192 Int nnz, /* number of entries to check */
193 Int mkind, /* R: 0, P: 1, C: 2, I: 3 */
194 double *Ax, /* NULL, size nnz, or 2*nnz */
195 double *Az /* NULL or size nnz */
196 ) ;
197
198 PRIVATE Int RB(extract) /* return # of explicit zero entries */
199 (
200 /* input */
201 Int ncol,
202 Int mkind, /* R: 0, P: 1, C: 2, I: 3 */
203
204 /* input/output */
205 Int *Ap, /* size ncol+1, column pointers A A */
206 Int *Ai, /* size nnz=Ap[ncol], row indices of A */
207 double *Ax, /* NULL, size nnz, or 2*nnz */
208 double *Az, /* NULL or size nnz */
209
210 /* output */
211 Int *Zp, /* size ncol+1, column pointers for Z */
212 Int *Zi /* size znz = Zp [ncol] = # of zeros in A on input */
213 ) ;
214
215 PRIVATE Int RB(readline) /* return string length, or -1 on error or EOF */
216 (
217 char *s, /* buffer to read into */
218 Int slen, /* s [0..slen] is of length slen+1 */
219 FILE *file /* file to read from */
220 ) ;
221
222 PRIVATE void RB(substring)
223 (
224 /* input */
225 char *s, /* input string */
226 Int len, /* length of string s [0..len] */
227 Int start, /* extract t [0:len-1] from s [start:start+len-1] */
228 Int tlen, /* length of substring t, excluding null terminator */
229
230 /* output */
231 char *t /* size tlen+1 */
232 ) ;
233
234 PRIVATE Int RB(xtoken) /* TRUE if token found, FALSE othewise */
235 (
236 /* input/output */
237 char *s, /* parse the next token in s [k..len] and update k */
238 Int len, /* length of s (input only) */
239 Int *k, /* start parsing at s [k] */
240 /* output */
241 double *x /* value of the token, or 0 if not found */
242 ) ;
243
244 PRIVATE Int RB(itoken)
245 (
246 /* input/output */
247 char *s, /* parse the next token in s [k..len] and update k */
248 Int len, /* length of s (input only) */
249 Int *k, /* start parsing at s [k] */
250 /* output */
251 Int *i /* value of the token, or 0 if not found */
252 ) ;
253
254 PRIVATE void RB(prune_space)
255 (
256 /* input/output */
257 char *s
258 ) ;
259
260 PRIVATE Int RB(header) /* 0: success, < 0: error, > 0: warning */
261 (
262 /* input */
263 FILE *file, /* must be already open for reading */
264
265 /* output */
266 char title [73], /* title, from first line of header */
267 char key [9], /* 8-character key, from first header line */
268 char mtype [4], /* RUA, RSA, PUA, PSA, RRA, etc */
269 Int *nrow, /* A is nrow-by-ncol */
270 Int *ncol,
271 Int *nnz, /* number of entries in A (tril(A) if symmetric) */
272 Int *nelnz, /* number of finite-elements */
273 char ptrfmt [21], /* Fortran format for column pointers */
274 char indfmt [21], /* Fortran format for row indices */
275 char valfmt [21], /* Fortran format for numerical values */
276 Int *mkind, /* R__: 0, P__: 1, C__: 2, I__: 3 */
277 Int *skind, /* _R_: -1, _U_: 0, _S_: 1, _H_: 2, _Z_: 3 */
278 Int *fem, /* __A: false, __E: true */
279 char *s, /* first line of data after the header */
280 Int slen /* s is of length slen+1 */
281 ) ;
282
283 PRIVATE Int RB(read1i) /* TRUE if OK, false otherwise */
284 (
285 FILE *file, /* file to read from (must be already open) */
286 char *s, /* buffer to use */
287 Int *len, /* strlen(s) */
288 Int slen, /* size of s */
289 Int *k, /* read position in s */
290 Int *x /* value read from the file */
291 ) ;
292
293 PRIVATE Int RB(read1x) /* TRUE if OK, false otherwise */
294 (
295 FILE *file, /* file to read from (must be already open) */
296 char *s, /* buffer to use */
297 Int *len, /* strlen(s) */
298 Int slen, /* size of s */
299 Int *k, /* read position in s */
300 double *x /* value read from the file */
301 ) ;
302
303 PRIVATE Int RB(iread) /* TRUE if OK, false otherwise */
304 (
305 /* input */
306 FILE *file, /* file to read from (must be already open) */
307 Int n, /* number of integers to read */
308 Int offset, /* offset to subtract from each integer */
309
310 /* output */
311 Int *A, /* read integers from file into A [0..n-1] */
312
313 /* input/workspace */
314 char *s, /* first line of input may be present in s */
315 Int slen /* s is of size slen+1 */
316 ) ;
317
318 PRIVATE Int RB(xread) /* TRUE if OK, false otherwise */
319 (
320 /* input */
321 FILE *file, /* file to read from (must be already open) */
322 Int n, /* number of values to read (2*n if complex) */
323 Int mkind, /* R: 0, P: 1, C: 2, I: 3 */
324
325 /* output */
326 double *Ax, /* read reals from file into Ax [0..n-1] */
327 double *Az, /* read reals from file into Az [0..n-1], may be NULL */
328
329 /* input/workspace */
330 char *s, /* first line of input may be present in s */
331 Int slen /* s is of size slen+1 */
332 ) ;
333
334 PRIVATE void RB(skipheader)
335 (
336 char *s,
337 Int slen,
338 FILE *file
339 ) ;
340
341
342 /* ========================================================================== */
343 /* === functions ============================================================ */
344 /* ========================================================================== */
345
346
347 /* -------------------------------------------------------------------------- */
348 /* RBget_entry: get numerical entry in the matrix at position p */
349 /* -------------------------------------------------------------------------- */
350
RB(get_entry)351 PUBLIC void RB(get_entry)
352 (
353 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
354 double *Ax, /* real part, or both if merged-complex */
355 double *Az, /* imaginary part if split-complex */
356 Int p, /* index of the entry */
357 double *xr, /* real part */
358 double *xz /* imaginary part */
359 )
360 {
361 if (mkind == 0 || mkind == 3)
362 {
363 /* A is real or integer */
364 *xr = Ax ? Ax [p] : 1 ;
365 *xz = 0 ;
366 }
367 else if (mkind == 2)
368 {
369 /* A is split-complex */
370 *xr = Ax ? Ax [p] : 1 ;
371 *xz = Ax ? Az [p] : 0 ;
372 }
373 else if (mkind == 4)
374 {
375 /* A is merged-complex */
376 *xr = Ax ? Ax [2*p ] : 1 ;
377 *xz = Ax ? Ax [2*p+1] : 0 ;
378 }
379 else /* if (mkind == 1) */
380 {
381 /* A is pattern-only */
382 *xr = 1 ;
383 *xz = 0 ;
384 }
385 }
386
387
388 /* -------------------------------------------------------------------------- */
389 /* RBput_entry: put numerical entry in the matrix in position p */
390 /* -------------------------------------------------------------------------- */
391
RB(put_entry)392 PUBLIC void RB(put_entry)
393 (
394 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
395 double *Ax, /* real part, or both if merged-complex */
396 double *Az, /* imaginary part if split-complex */
397 Int p, /* index of the entry */
398 double xr, /* real part */
399 double xz /* imaginary part */
400 )
401 {
402 if (mkind == 0 || mkind == 3)
403 {
404 /* A is real or integer; ignore xz */
405 if (Ax) Ax [p] = xr ;
406 if (Az) Az [p] = 0 ;
407 }
408 else if (mkind == 2)
409 {
410 /* A is split-complex */
411 if (Ax) Ax [p] = xr ;
412 if (Az) Az [p] = xz ;
413 }
414 else if (mkind == 4)
415 {
416 /* A is merged-complex */
417 if (Ax) Ax [2*p ] = xr ;
418 if (Ax) Ax [2*p+1] = xz ;
419 }
420 else /* if (mkind == 1) */
421 {
422 /* A is pattern-only; ignore xr and xz */
423 if (Ax) Ax [p] = 1 ;
424 if (Az) Az [p] = 0 ;
425 }
426 }
427
428
429 /* -------------------------------------------------------------------------- */
430 /* RBskipheader: skip past the header */
431 /* -------------------------------------------------------------------------- */
432
RB(skipheader)433 PRIVATE void RB(skipheader)
434 (
435 char *s, /* size slen+1 */
436 Int slen,
437 FILE *file /* file to read from */
438 )
439 {
440 s [0] = '\0' ;
441 RB(readline) (s, slen, file) ;
442 RB(readline) (s, slen, file) ;
443 RB(readline) (s, slen, file) ;
444 RB(readline) (s, slen, file) ;
445 RB(readline) (s, slen, file) ;
446 if ((s [0] == 'F') || (s [0] == 'f') || (s [0] == 'M') || (s[0] == 'm'))
447 {
448 RB(readline) (s, slen, file) ;
449 }
450 /* s now contains the first non-header line */
451 }
452
453
454 /* -------------------------------------------------------------------------- */
455 /* RBread2: read all but the header and construct the matrix */
456 /* -------------------------------------------------------------------------- */
457
RB(read2)458 PRIVATE Int RB(read2) /* 0: OK, < 0: error, > 0: warning */
459 (
460 /* input */
461 FILE *file, /* must be already open for reading */
462 Int nrow, /* A is nrow-by-ncol */
463 Int ncol,
464 Int nnz, /* number of entries in A, from the header */
465 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
466 Int skind, /* R: -1, U: 0, S: 1, H: 2, Z: 3 */
467 Int build_upper, /* if TRUE and skind>0, then create upper part. Ai, Ax,
468 and Az must be twice the size as given below */
469
470 /* output */
471 Int *Ap, /* size ncol+1, column pointers for A */
472 Int *Ai, /* size nnz, row indices for A */
473 double *Ax, /* size nnz or 2*nnz if complex and Az NULL */
474 double *Az, /* size nnz, or NULL, for complex matrices only */
475
476 /* workspace */
477 Int *w, /* size MAX(nrow,ncol)+1 */
478 Int *cp, /* size ncol+1 */
479
480 /* input/workspace */
481 char *s, /* first line of column pointers on input */
482 Int slen
483 )
484 {
485 double xr = 0, xz = 0 ;
486 Int p, i, j, k, ilast, alen, llen, psrc, pdst ;
487
488 /* ---------------------------------------------------------------------- */
489 /* skip past the header, if reading from a file */
490 /* ---------------------------------------------------------------------- */
491
492 if (file)
493 {
494 RB(skipheader) (s, slen, file) ;
495 }
496
497 /* ---------------------------------------------------------------------- */
498 /* read the column pointers and check them */
499 /* ---------------------------------------------------------------------- */
500
501 if (!RB(iread) (file, ncol+1, 1, Ap, s, slen))
502 {
503 return (RBIO_CP_IOERROR) ; /* I/O error reading column pointers */
504 }
505 if (Ap [0] != 0 || Ap [ncol] != nnz)
506 {
507 return (RBIO_CP_INVALID) ; /* column pointers invalid */
508 }
509 for (j = 1 ; j <= ncol ; j++)
510 {
511 if (Ap [j] < Ap [j-1])
512 {
513 return (RBIO_CP_INVALID) ;
514 }
515 }
516
517 /* ---------------------------------------------------------------------- */
518 /* read the row indices and check them */
519 /* ---------------------------------------------------------------------- */
520
521 if (!RB(iread) (file, nnz, 1, Ai, s, slen))
522 {
523 return (RBIO_ROW_IOERROR) ; /* I/O error reading row indices */
524 }
525
526 for (i = 0 ; i < nrow ; i++)
527 {
528 w [i] = -1 ;
529 }
530
531 for (j = 0 ; j < ncol ; j++)
532 {
533 ilast = -1 ;
534 for (p = Ap [j] ; p < Ap [j+1] ; p++)
535 {
536 i = Ai [p] ;
537 if (i < 0 || i >= nrow)
538 {
539 return (RBIO_ROW_INVALID) ; /* row index out of range */
540 }
541 if (w [i] == j)
542 {
543 return (RBIO_DUPLICATE) ; /* duplicate entry in matrix */
544 }
545 w [i]= j ;
546 if (i < ilast)
547 {
548 return (RBIO_JUMBLED) ; /* row indices unsorted */
549 }
550 ilast = i ;
551 }
552 }
553
554 /* ---------------------------------------------------------------------- */
555 /* read the values */
556 /* ---------------------------------------------------------------------- */
557
558 if (!RB(xread) (file, nnz, mkind, Ax, Az, s, slen))
559 {
560 return (RBIO_VALUE_IOERROR) ; /* I/O error reading values */
561 }
562
563 /* ---------------------------------------------------------------------- */
564 /* construct or check the upper triangular part for symmetric matrices */
565 /* ---------------------------------------------------------------------- */
566
567 if (skind > 0)
568 {
569
570 /* ------------------------------------------------------------------ */
571 /* check the matrix and compute new column counts */
572 /* ------------------------------------------------------------------ */
573
574 for (j = 0 ; j <= ncol ; j++)
575 {
576 w [j] = 0 ;
577 }
578
579 for (j = 0 ; j < ncol ; j++)
580 {
581 for (p = Ap [j] ; p < Ap [j+1] ; p++)
582 {
583 i = Ai [p] ;
584 if (i == j)
585 {
586 /* diagonal entry, only appears as A(j,j) */
587 w [j]++ ;
588 }
589 else if (i > j)
590 {
591 /* entry in lower triangular part, A(i,j) will be */
592 /* duplicated as A(j,i), so count it in column i and j */
593 w [i]++ ;
594 w [j]++ ;
595 }
596 else
597 {
598 /* error: entry in upper triangular part */
599 return (RBIO_EXTRANEOUS) ;
600 }
601 }
602 }
603
604 /* ------------------------------------------------------------------ */
605 /* construct the upper trianglar part, if requested */
606 /* ------------------------------------------------------------------ */
607
608 if (build_upper)
609 {
610
611 /* -------------------------------------------------------------- */
612 /* compute the new column pointers */
613 /* -------------------------------------------------------------- */
614
615 cp [0] = 0 ;
616 for (j = 1 ; j <= ncol ; j++)
617 {
618 cp [j] = cp [j-1] + w [j-1] ;
619 }
620
621 /* -------------------------------------------------------------- */
622 /* shift the matrix by adding gaps to the top of each column */
623 /* -------------------------------------------------------------- */
624
625 for (j = ncol-1 ; j >= 0 ; j--)
626 {
627 /* number of entries in lower tri. part (including diagonal) */
628 llen = Ap [j+1] - Ap [j] ;
629
630 /* number of entries in entire column */
631 alen = cp [j+1] - cp [j] ;
632
633 /* move the column from Ai [Ap [j] ... Ap [j+1]-1] down to
634 Ai [cp [j+1]-llen ... cp[j+1]-1], leaving a gap
635 at Ai [Ap [j] ... cp [j+1]-llen] */
636
637 for (k = 1 ; k <= llen ; k++)
638 {
639 psrc = Ap [j+1] - k ;
640 pdst = cp [j+1] - k ;
641 Ai [pdst] = Ai [psrc] ;
642 RB(get_entry) (mkind, Ax, Az, psrc, &xr, &xz) ;
643 RB(put_entry) (mkind, Ax, Az, pdst, xr, xz) ;
644 }
645 }
646
647 /* -------------------------------------------------------------- */
648 /* populate the upper triangular part */
649 /* -------------------------------------------------------------- */
650
651 /* create temporary column pointers to point to the gaps */
652 for (j = 0 ; j < ncol ; j++)
653 {
654 w [j] = cp [j] ;
655 }
656
657 for (j = 0 ; j < ncol ; j++)
658 {
659 /* scan the entries in the lower tri. part, in
660 Ai [cp[j+1]-llen ... cp[j+1]-1] */
661 llen = Ap [j+1] - Ap [j] ;
662
663 for (k = 1 ; k <= llen ; k++)
664 {
665 /* get the A(i,j) entry in strictly lower triangular part */
666 psrc = cp [j+1] - k ;
667 i = Ai [psrc] ;
668 if (i != j)
669 {
670 /* get the numerical value of the A(i,j) entry */
671 RB(get_entry) (mkind, Ax, Az, psrc, &xr, &xz) ;
672
673 /* negate for Hermitian and skew-symmetric */
674 if (skind == 2)
675 {
676 xz = -xz ; /* complex Hermitian */
677 }
678 else if (skind == 3)
679 {
680 xr = -xr ; /* skew symmetric */
681 xz = -xz ;
682 }
683
684 /* add A(j,i) as the next entry in column i */
685 pdst = w [i]++ ;
686 Ai [pdst] = j ;
687 RB(put_entry) (mkind, Ax, Az, pdst, xr, xz) ;
688 }
689 }
690 }
691
692 /* -------------------------------------------------------------- */
693 /* finalize the column pointers */
694 /* -------------------------------------------------------------- */
695
696 for (j = 0 ; j <= ncol ; j++)
697 {
698 Ap [j] = cp [j] ;
699 }
700 }
701 }
702
703 /* ---------------------------------------------------------------------- */
704 /* matrix is valid */
705 /* ---------------------------------------------------------------------- */
706
707 return (RBIO_OK) ;
708 }
709
710
711 /* -------------------------------------------------------------------------- */
712 /* RBzcount: count the number of explicit zeros */
713 /* -------------------------------------------------------------------------- */
714
RB(zcount)715 PRIVATE Int RB(zcount) /* return number of explicit zeros in A */
716 (
717 Int nnz, /* number of entries in A */
718 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
719 double *Ax, /* NULL, size nnz, or 2*nnz */
720 double *Az /* NULL or size nnz */
721 )
722 {
723 double xr, xz ;
724 Int p, znz = 0 ;
725 for (p = 0 ; p < nnz ; p++)
726 {
727 RB(get_entry) (mkind, Ax, Az, p, &xr, &xz) ;
728 if (xr == 0 && xz == 0)
729 {
730 znz++ ;
731 }
732 }
733 return (znz) ;
734 }
735
736
737 /* -------------------------------------------------------------------------- */
738 /* RBextract: extract explicit zeros */
739 /* -------------------------------------------------------------------------- */
740
741 /* The matrix is A is split into A (the nonzeros) and Z (the pattern of
742 explicit zero entries in A). On input, A may have explicit zeros. On
743 output they are removed from A. If Zp is NULL on input, these explicit
744 zeros are removed from A and simply discarded. */
745
RB(extract)746 PRIVATE Int RB(extract) /* return # of explicit zero entries */
747 (
748 /* input */
749 Int ncol,
750 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
751
752 /* input/output */
753 Int *Ap, /* size ncol+1, column pointers A A */
754 Int *Ai, /* size nnz=Ap[ncol], row indices of A */
755 double *Ax, /* NULL, size nnz, or 2*nnz */
756 double *Az, /* NULL or size nnz */
757
758 /* output */
759 Int *Zp, /* size ncol+1, column pointers for Z */
760 Int *Zi /* size znz = Zp [ncol] = # of zeros in A on input */
761 )
762 {
763 double xr, xz ;
764 Int p, i, j, pa, pz ;
765
766 pa = 0 ;
767 pz = 0 ;
768
769 for (j = 0 ; j < ncol ; j++)
770 {
771 /* save the new start of column j of A */
772 p = Ap [j] ;
773 Ap [j] = pa ;
774 if (Zp) Zp [j] = pz ;
775
776 /* split column j of A */
777 for ( ; p < Ap [j+1] ; p++)
778 {
779 i = Ai [p] ;
780 RB(get_entry) (mkind, Ax, Az, p, &xr, &xz) ;
781 if (xr == 0 && xz == 0)
782 {
783 /* copy into Z, if Z exists */
784 if (Zp) Zi [pz++] = i ;
785 }
786 else
787 {
788 /* copy into A */
789 Ai [pa] = i ;
790 RB(put_entry) (mkind, Ax, Az, pa++, xr, xz) ;
791 }
792 }
793 }
794 Ap [ncol] = pa ;
795 if (Zp) Zp [ncol] = pz ;
796 return (pz) ;
797 }
798
799
800 /* -------------------------------------------------------------------------- */
801 /* RBreadline: read a line from the file */
802 /* -------------------------------------------------------------------------- */
803
RB(readline)804 PRIVATE Int RB(readline) /* return length of string, or -1 on error or EOF */
805 (
806 char *s, /* buffer to read into */
807 Int slen, /* s [0..slen] is of length slen+1 */
808 FILE *file /* file to read from */
809 )
810 {
811 Int len, ok ;
812 if (!file) file = stdin ; /* file defaults to stdin */
813 ok = (fgets (s, slen, file) != NULL) ; /* read line in s [0..slen-1] */
814 s [slen] = '\0' ; /* ensure line is terminated */
815 len = ok ? ((Int) strlen (s)) : (-1) ; /* get the length */
816 return ((len < slen) ? len : (-1)) ; /* return len, or -1 on error */
817 }
818
819
820 /* -------------------------------------------------------------------------- */
821 /* RBsubstring: extract a substring */
822 /* -------------------------------------------------------------------------- */
823
RB(substring)824 PRIVATE void RB(substring)
825 (
826 /* input */
827 char *s, /* input string */
828 Int len, /* length of string s [0..len] */
829 Int start, /* extract t [0:len-1] from s [start:start+len-1] */
830 Int tlen, /* length of substring t, excluding null terminator */
831
832 /* output */
833 char *t /* size tlen+1 */
834 )
835 {
836 Int i, end, k ;
837 end = MIN (start+tlen, len) ; /* do not go past the end of s */
838 k = 0 ;
839 for (i = start ; i < end ; i++)
840 {
841 t [k++] = s [i] ;
842 }
843 t [k] = '\0' ; /* ensure t is null-terminated */
844 }
845
846
847 /* -------------------------------------------------------------------------- */
848 /* RBxtoken: get the next token from a string, returning it as a double */
849 /* -------------------------------------------------------------------------- */
850
851 /* On input, the token to return from s [0..len] starts at s[k]. On output, k
852 is updated so that the s[k] is the start of the token after this one. */
853
RB(xtoken)854 PRIVATE Int RB(xtoken) /* TRUE if token found, FALSE othewise */
855 (
856 /* input/output */
857 char *s, /* parse the next token in s [k..len] and update k */
858 Int len, /* length of s (input only) */
859 Int *k, /* start parsing at s [k] */
860 /* output */
861 double *x /* value of the token, or 0 if not found */
862 )
863 {
864 Int start ;
865
866 *x = 0 ;
867
868 /* consume leading spaces, if any */
869 while ((*k) < len && s [*k] == ' ')
870 {
871 (*k)++ ;
872 }
873
874 /* the token starts here, if present */
875 start = (*k) ;
876
877 if (s [start] == '\0')
878 {
879 /* the end of s has been reached; there is no token */
880 return (FALSE) ;
881 }
882
883 /* find where the token ends */
884 while ((*k) < len && s [*k] != ' ')
885 {
886 (*k)++ ;
887 }
888
889 /* terminate the current token. Note this might be s [len] */
890 if (s [*k] != '\0')
891 {
892 s [(*k)++] = '\0' ;
893 }
894
895 /* parse the current token and return its value */
896 return (sscanf (s+start, "%lg", x) == 1) ;
897 }
898
899
900 /* -------------------------------------------------------------------------- */
901 /* RBitoken: get the next token from a string, returning it as an Int */
902 /* -------------------------------------------------------------------------- */
903
904 /* Same as RBxtoken, except for the data type of the 4th parameter. The
905 return value i is checked for integer overflow of i+1. */
906
RB(itoken)907 PRIVATE Int RB(itoken)
908 (
909 /* input/output */
910 char *s, /* parse the next token in s [k..len] and update k */
911 Int len, /* length of s (input only) */
912 Int *k, /* start parsing at s [k] */
913 /* output */
914 Int *i /* value of the token, or 0 if not found */
915 )
916 {
917 double x ;
918 Int ok = RB(xtoken) (s, len, k, &x) ;
919 *i = (Int) x ; /* convert to integer */
920 return (ok && ((double) ((*i)+1) == (x+1))) ;
921 }
922
923 /* -------------------------------------------------------------------------- */
924 /* RBprune_space: remove trailing space from a string */
925 /* -------------------------------------------------------------------------- */
926
RB(prune_space)927 PRIVATE void RB(prune_space)
928 (
929 /* input/output */
930 char *s
931 )
932 {
933 Int k ;
934 for (k = strlen (s) - 1 ; k >= 0 ; k--)
935 {
936 if (isspace (s [k]))
937 {
938 s [k] = '\0' ;
939 }
940 else
941 {
942 return ;
943 }
944 }
945 }
946
947
948 /* -------------------------------------------------------------------------- */
949 /* RBheader: read Rutherford/Boeing header lines */
950 /* -------------------------------------------------------------------------- */
951
952 /*
953 The Rutherford/Boeing matrix type is a 3-character string:
954
955 (1) R: real, C: complex, P: pattern only, I: integer
956 mkind: 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged
957 Note that this function does not return mkind == 4, since the
958 split/merged Ax and Az are considered by this function.
959 It must be selected by the caller.
960
961 (2) S: symmetric, U: unsymmetric, H: Hermitian, Z: skew symmetric,
962 R: rectangular
963 skind: R: -1, U: 0, S: 1, H: 2, Z: 3
964
965 (3) A: assembled, E: element form
966 nelnz = 0 for A, number of elements for E
967
968 pattern matrices are given numerical values of 1 (except PZA).
969 PZA matrices have +1 in the lower triangular part and -1 in
970 the upper triangular part.
971 */
972
RB(header)973 PRIVATE Int RB(header) /* 0: success, < 0: error, > 0: warning */
974 (
975 /* input */
976 FILE *file, /* must be already open for reading */
977
978 /* output */
979 char title [73], /* title, from first line of header */
980 char key [9], /* 8-character key, from first header line */
981 char mtype [4], /* RUA, RSA, PUA, PSA, RRA, etc */
982 Int *nrow, /* A is nrow-by-ncol */
983 Int *ncol,
984 Int *nnz, /* number of entries in A (tril(A) if symmetric) */
985 Int *nelnz, /* number of finite-elements */
986 char ptrfmt [21], /* Fortran format for column pointers */
987 char indfmt [21], /* Fortran format for row indices */
988 char valfmt [21], /* Fortran format for numerical values */
989 Int *mkind, /* R__: 0, P__: 1, C__: 2, I__: 3 */
990 Int *skind, /* _R_: -1, _U_: 0, _S_: 1, _H_: 2, _Z_: 3 */
991 Int *fem, /* __A: false, __E: true */
992 char *s, /* first line of data after the header */
993 Int slen /* s is of length slen+1 */
994 )
995 {
996 Int len, totcrd, ptrcrd, indcrd, valcrd, k, ok ;
997
998 /* ---------------------------------------------------------------------- */
999 /* get the 1st line: title, key */
1000 /* ---------------------------------------------------------------------- */
1001
1002 len = RB(readline) (s, slen, file) ;
1003 if (len < 0) return (RBIO_HEADER_IOERROR) ;
1004
1005 RB(substring) (s, len, 0, 72, title) ;
1006 RB(substring) (s, len, 72, 8, key) ;
1007
1008 if (title [71] == '|')
1009 {
1010 /* remove the marker placed by RBwrite, if present */
1011 title [71] = '\0' ;
1012 }
1013
1014 /* remove trailing spaces from the key and the title */
1015 RB(prune_space) (title) ;
1016 RB(prune_space) (key) ;
1017
1018 /* ---------------------------------------------------------------------- */
1019 /* get the 2nd line: totcrd, ptrcrd, indcrd, valcrd */
1020 /* ---------------------------------------------------------------------- */
1021
1022 len = RB(readline) (s, slen, file) ;
1023 if (len < 0) return (RBIO_HEADER_IOERROR) ;
1024
1025 k = 0 ;
1026 ok = TRUE ;
1027 ok = ok && RB(itoken) (s, len, &k, &totcrd) ;
1028 ok = ok && RB(itoken) (s, len, &k, &ptrcrd) ;
1029 ok = ok && RB(itoken) (s, len, &k, &indcrd) ;
1030 ok = ok && RB(itoken) (s, len, &k, &valcrd) ;
1031
1032 /* ---------------------------------------------------------------------- */
1033 /* get the 3rd line: mtype (RUA, etc), nrow, ncol, nnz, nelnz */
1034 /* ---------------------------------------------------------------------- */
1035
1036 len = RB(readline)(s, slen, file) ;
1037 if (len < 0) return (RBIO_HEADER_IOERROR) ;
1038
1039 RB(substring) (s, len, 0, 3, mtype) ;
1040
1041 k = 3 ;
1042 ok = ok && RB(itoken) (s, len, &k, nrow) ;
1043 ok = ok && RB(itoken) (s, len, &k, ncol) ;
1044 ok = ok && RB(itoken) (s, len, &k, nnz) ;
1045 RB(itoken) (s, len, &k, nelnz) ;
1046
1047 if (!ok || *nrow <= 0 || *ncol <= 0 || *nnz <= 0)
1048 {
1049 return (RBIO_DIM_INVALID) ;
1050 }
1051
1052 /* ---------------------------------------------------------------------- */
1053 /* get the 4th line: ptrfmt, indfmt, valfmt */
1054 /* ---------------------------------------------------------------------- */
1055
1056 len = RB(readline) (s, slen, file) ;
1057 if (len < 0) return (RBIO_HEADER_IOERROR) ;
1058
1059 RB(substring) (s, len, 0, 16, ptrfmt) ;
1060 RB(substring) (s, len, 16, 16, indfmt) ;
1061 RB(substring) (s, len, 32, 20, valfmt) ;
1062
1063 /* ---------------------------------------------------------------------- */
1064 /* skip the Harwell/Boeing header line 5, if present */
1065 /* ---------------------------------------------------------------------- */
1066
1067 len = RB(readline) (s, slen, file) ;
1068 if (len <= 0) return (RBIO_HEADER_IOERROR) ;
1069
1070 if ((s [0] == 'F') || (s [0] == 'f') || (s [0] == 'M') || (s [0] == 'm'))
1071 {
1072 len = RB(readline) (s, slen, file) ;
1073 if (len <= 0) return (RBIO_HEADER_IOERROR) ;
1074 }
1075
1076 /* ---------------------------------------------------------------------- */
1077 /* determine if real, pattern, integer, or complex */
1078 /* ---------------------------------------------------------------------- */
1079
1080 if (mtype [0] == 'R' || mtype [0] == 'r')
1081 {
1082 *mkind = 0 ; /* real */
1083 }
1084 else if (mtype [0] == 'P' || mtype [0] == 'p')
1085 {
1086 *mkind = 1 ; /* pattern */
1087 }
1088 else if (mtype [0] == 'C' || mtype [0] == 'c')
1089 {
1090 *mkind = 2 ; /* complex (assume split-complex) */
1091 }
1092 else if (mtype [0] == 'I' || mtype [0] == 'i')
1093 {
1094 *mkind = 3 ; /* integer */
1095 }
1096 else
1097 {
1098 return (RBIO_TYPE_INVALID) ; /* error: invalid matrix mtype */
1099 }
1100
1101 /* ---------------------------------------------------------------------- */
1102 /* determine the symmetry property */
1103 /* ---------------------------------------------------------------------- */
1104
1105 if (mtype [1] == 'R' || mtype [1] == 'r')
1106 {
1107 *skind = -1 ; /* rectangular: RRA, PRA, IRA, and CRA matrices */
1108 }
1109 else if (mtype [1] == 'U' || mtype [1] == 'u')
1110 {
1111 *skind = 0 ; /* unsymmetric: RUA, PUA, IUA, and CUA matrices */
1112 }
1113 else if (mtype [1] == 'S' || mtype [1] == 's')
1114 {
1115 *skind = 1 ; /* symmetric: RSA, PSA, ISA, and CSA matrices */
1116 }
1117 else if (mtype [1] == 'H' || mtype [1] == 'h')
1118 {
1119 *skind = 2 ; /* Hermitian: CHA (PHA, IHA, and RHA are OK too) */
1120 }
1121 else if (mtype [1] == 'Z' || mtype [1] == 'z')
1122 {
1123 *skind = 3 ; /* skew symmetric: RZA, PZA, IZA, and CZA */
1124 }
1125 else
1126 {
1127 return (RBIO_TYPE_INVALID) ; /* error: invalid matrix mtype */
1128 }
1129
1130 /* ---------------------------------------------------------------------- */
1131 /* assembled vs elemental matrices (**A vs **E) */
1132 /* ---------------------------------------------------------------------- */
1133
1134 if (mtype [2] == 'A' || mtype [2] == 'a')
1135 {
1136 /* assembled - ignore nelnz */
1137 *fem = FALSE ;
1138 *nelnz = 0 ;
1139 }
1140 else if (mtype [2] == 'E' || mtype [2] == 'e')
1141 {
1142 /* finite-element */
1143 *fem = TRUE ;
1144 }
1145 else
1146 {
1147 return (RBIO_TYPE_INVALID) ; /* error: invalid matrix mtype */
1148 }
1149
1150 /* ---------------------------------------------------------------------- */
1151 /* assembled matrices must be square if skind is not R */
1152 /* ---------------------------------------------------------------------- */
1153
1154 if (!(*fem) && (*skind) != -1 && (*nrow) != (*ncol))
1155 {
1156 return (RBIO_DIM_INVALID) ; /* error: invalid matrix dimensions */
1157 }
1158
1159 /* ---------------------------------------------------------------------- */
1160 /* matrix header is valid */
1161 /* ---------------------------------------------------------------------- */
1162
1163 return (0) ;
1164 }
1165
1166
1167 /* -------------------------------------------------------------------------- */
1168 /* RBread1i: read a single integer value from the file */
1169 /* -------------------------------------------------------------------------- */
1170
RB(read1i)1171 PRIVATE Int RB(read1i) /* TRUE if OK, false otherwise */
1172 (
1173 FILE *file, /* file to read from (must be already open) */
1174 char *s, /* buffer to use */
1175 Int *len, /* strlen(s) */
1176 Int slen, /* size of s */
1177 Int *k, /* read position in s */
1178 Int *x /* value read from the file */
1179 )
1180 {
1181 /* read the token from the current line */
1182 Int ok = RB(itoken) (s, *len, k, x) ;
1183 if (!ok)
1184 {
1185 /* input line is exhausted; get the next one */
1186 *len = RB(readline) (s, slen, file) ;
1187 if (*len < 0) return (FALSE) ;
1188 /* read first token from the new line */
1189 *k = 0 ;
1190 ok = RB(itoken) (s, *len, k, x) ;
1191 }
1192 return (ok) ;
1193 }
1194
1195
1196 /* -------------------------------------------------------------------------- */
1197 /* RBread1x: read a single real value from the file */
1198 /* -------------------------------------------------------------------------- */
1199
RB(read1x)1200 PRIVATE Int RB(read1x) /* TRUE if OK, false otherwise */
1201 (
1202 FILE *file, /* file to read from (must be already open) */
1203 char *s, /* buffer to use */
1204 Int *len, /* strlen(s) */
1205 Int slen, /* size of s */
1206 Int *k, /* read position in s */
1207 double *x /* value read from the file */
1208 )
1209 {
1210 /* read the token from the current line */
1211 Int ok = RB(xtoken) (s, *len, k, x) ;
1212 if (!ok)
1213 {
1214 /* input line is exhausted; get the next one */
1215 *len = RB(readline) (s, slen, file) ;
1216 if (*len < 0) return (FALSE) ;
1217 /* read first token from the new line */
1218 *k = 0 ;
1219 ok = RB(xtoken) (s, *len, k, x) ;
1220 }
1221 return (ok) ;
1222 }
1223
1224
1225 /* -------------------------------------------------------------------------- */
1226 /* RBiread: read n integers from the file */
1227 /* -------------------------------------------------------------------------- */
1228
RB(iread)1229 PRIVATE Int RB(iread) /* TRUE if OK, FALSE otherwise */
1230 (
1231 /* input */
1232 FILE *file, /* file to read from (must be already open) */
1233 Int n, /* number of integers to read */
1234 Int offset, /* offset to subtract from each integer */
1235
1236 /* output */
1237 Int *A, /* read integers from file into A [0..n-1] */
1238
1239 /* input/workspace */
1240 char *s, /* first line of input may be present in s */
1241 Int slen /* s is of size slen+1 */
1242 )
1243 {
1244 Int k, len, i, x, ok = TRUE ;
1245
1246 /* the first token will be read from the current line, s, on input */
1247 len = strlen (s) ;
1248 k = 0 ;
1249
1250 for (i = 0 ; ok && i < n ; i++)
1251 {
1252 /* read the next integer */
1253 ok = RB(read1i) (file, s, &len, slen, &k, &x) ;
1254 A [i] = x - offset ;
1255 }
1256 s [0] = '\0' ; /* ignore remainder of current input string */
1257 return (ok) ;
1258 }
1259
1260
1261 /* -------------------------------------------------------------------------- */
1262 /* RBxread: read n real or complex values from the file */
1263 /* -------------------------------------------------------------------------- */
1264
RB(xread)1265 PRIVATE Int RB(xread) /* TRUE if OK, FALSE otherwise */
1266 (
1267 /* input */
1268 FILE *file, /* file to read from (must be already open) */
1269 Int n, /* number of values to read (2*n if complex) */
1270 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
1271
1272 /* output */
1273 double *Ax, /* read reals from file into Ax [0..n-1] */
1274 double *Az, /* read reals from file into Az [0..n-1], may be NULL */
1275
1276 /* input/workspace */
1277 char *s, /* first line of input may be present in s */
1278 Int slen /* s is of size slen+1 */
1279 )
1280 {
1281 double xr, xz ;
1282 Int k, len, i, ok = TRUE ;
1283
1284 /* the first token will be read from the current line, s, on input */
1285 len = strlen (s) ;
1286 k = 0 ;
1287 xr = 1 ;
1288 xz = 0 ;
1289
1290 for (i = 0 ; ok && i < n ; i++)
1291 {
1292 /* read the real and imaginary part of the kth entry */
1293 if (ok && mkind != 1)
1294 {
1295 /* read the real part, if present */
1296 ok = RB(read1x) (file, s, &len, slen, &k, &xr) ;
1297 }
1298 if (ok && (mkind == 2 || mkind == 4))
1299 {
1300 /* read the imaginary part, if present */
1301 ok = RB(read1x) (file, s, &len, slen, &k, &xz) ;
1302 }
1303 /* store the ith entry in Ax and/or Az */
1304 RB(put_entry) (mkind, Ax, Az, i, xr, xz) ;
1305 }
1306 s [0] = '\0' ; /* ignore remainder of current input string */
1307 return (ok) ;
1308 }
1309
1310
1311 /* -------------------------------------------------------------------------- */
1312 /* RBread: read a Rutherford/Boeing matrix from a file */
1313 /* -------------------------------------------------------------------------- */
1314
1315 /*
1316 Space is allocated as needed for the output matrix A, the zero pattern Z,
1317 and required workspace. To free the result of this function, do on the
1318 following:
1319
1320 if (Ap) free (Ap) ;
1321
1322 or
1323
1324 SuiteSparse_free (Ap) ;
1325
1326 etc, for Ap, Ai, Ax, Az, Zp, and Zi.
1327
1328 Format of the numerical values:
1329 If Ax is NULL on input, only the pattern of the matrix is returned. If
1330 Ax is not NULL but p_Az is NULL, then the real and imaginary parts (if
1331 complex) are both returned in Ax. If both Ax and p_Az are non-NULL,
1332 imaginary parts (if present) are returned in Az.
1333
1334 For real and integer matrices (R__, I__), Az is always returned as NULL.
1335
1336 For pattern matrices (P__), Ax and Az are always returned as NULL.
1337
1338 For complex matrices (C__), either Ax is used if p_Az is NULL
1339 (merged-complex) or both Ax and Az are used (split-complex).
1340
1341 Only assembled matrices (__A) are handled. Finite-element matrices
1342 (__E) are not handled.
1343 */
1344
RB(read)1345 PUBLIC Int RB(read) /* 0: OK, < 0: error, > 0: warning */
1346 (
1347 /* input */
1348 char *filename, /* filename to read from */
1349 Int build_upper, /* if true, construct upper part for sym. matrices */
1350 Int zero_handling, /* 0: do nothing, 1: prune zeros, 2: extract zeros */
1351
1352 /* output */
1353 char title [73],
1354 char key [9],
1355 char mtype [4], /* RUA, RSA, PUA, PSA, RRA, etc */
1356 Int *nrow, /* A is nrow-by-ncol */
1357 Int *ncol,
1358 Int *mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
1359 Int *skind, /* R: -1, U: 0, S: 1, H: 2, Z: 3 */
1360 Int *asize, /* Ai array has size asize*sizeof(double) */
1361 Int *znz, /* number of explicit zeros removed from A */
1362
1363 /* output: these are malloc'ed below and must be freed by the caller */
1364 Int **p_Ap, /* column pointers of A */
1365 Int **p_Ai, /* row indices of A */
1366 double **p_Ax, /* real values (ignored if NULL) of A */
1367 double **p_Az, /* imaginary values (ignored if NULL) of A */
1368 Int **p_Zp, /* column pointers of Z */
1369 Int **p_Zi /* row indices of Z */
1370 )
1371 {
1372 Int nnz, nelnz, status, fem ;
1373 Int *w, *cp, *Ap, *Ai, *Zp, *Zi ;
1374 double *Ax, *Az ;
1375 FILE *file = NULL ; /* read from stdin if NULL */
1376 int ok ;
1377 char s [SLEN+1], ptrfmt [21], indfmt [21], valfmt [21] ;
1378
1379 /* ---------------------------------------------------------------------- */
1380 /* check inputs */
1381 /* ---------------------------------------------------------------------- */
1382
1383 if (p_Ap) *p_Ap = NULL ;
1384 if (p_Ai) *p_Ai = NULL ;
1385 if (p_Ax) *p_Ax = NULL ;
1386 if (p_Az) *p_Az = NULL ;
1387 if (p_Zp) *p_Zp = NULL ;
1388 if (p_Zi) *p_Zi = NULL ;
1389
1390 if (!title || !key || !mtype || !p_Ap || !p_Ai || !nrow || !ncol || !mkind
1391 || !skind || (zero_handling == 2 && (!p_Zp || !p_Zi)) || !znz || !asize)
1392 {
1393 /* one or more required arguments are missing (NULL) */
1394 return (RBIO_ARG_ERROR) ;
1395 }
1396
1397 /* ---------------------------------------------------------------------- */
1398 /* read the header */
1399 /* ---------------------------------------------------------------------- */
1400
1401 if (filename)
1402 {
1403 file = fopen (filename, "r") ;
1404 if (file == NULL)
1405 {
1406 return (RBIO_FILE_IOERROR) ; /* cannot open file */
1407 }
1408 }
1409
1410 status = RB(header) (file, title, key, mtype, nrow, ncol, &nnz, &nelnz,
1411 ptrfmt, indfmt, valfmt, mkind, skind, &fem, s, SLEN) ;
1412
1413 /* close the file, so if there is an error the file is not left open */
1414 if (filename) fclose (file) ;
1415
1416 if (status != 0)
1417 {
1418 return (status) ;
1419 }
1420
1421 if (fem)
1422 {
1423 return (RBIO_UNSUPPORTED) ;
1424 }
1425
1426 /* ---------------------------------------------------------------------- */
1427 /* adjust mkind based on the presence of Ax and Az input arguments */
1428 /* ---------------------------------------------------------------------- */
1429
1430 if (!p_Ax)
1431 {
1432 /* the numerical values cannot be returned, so this is read as a
1433 pattern-only matrix */
1434 *mkind = 1 ;
1435 }
1436
1437 if (*mkind == 2 && !p_Az)
1438 {
1439 /* The matrix is complex. If p_Az is NULL, then the imaginary part
1440 will be returned in Ax, and the matrix is thus 4:merged-complex */
1441 *mkind = 4 ;
1442 }
1443
1444 /* ---------------------------------------------------------------------- */
1445 /* allocate space for A */
1446 /* ---------------------------------------------------------------------- */
1447
1448 *asize = ((build_upper) ? 2 : 1) * MAX (nnz, 1) ;
1449 Ap = (Int *) SuiteSparse_malloc ((*ncol) + 1, sizeof (Int)) ;
1450 Ai = (Int *) SuiteSparse_malloc (*asize, sizeof (Int)) ;
1451 ok = (Ap != NULL && Ai != NULL) ;
1452 Ax = NULL ;
1453 Az = NULL ;
1454 Zp = NULL ;
1455 Zi = NULL ;
1456
1457 if (*mkind == 0 || *mkind == 3 || *mkind == 1)
1458 {
1459 /* return A as real, integer, or pattern */
1460 if (p_Ax)
1461 {
1462 Ax = (double *) SuiteSparse_malloc (*asize, sizeof (double)) ;
1463 ok = ok && (Ax != NULL) ;
1464 }
1465 }
1466 else if (*mkind == 2)
1467 {
1468 /* return A as split-complex */
1469 Ax = (double *) SuiteSparse_malloc (*asize, sizeof (double)) ;
1470 Az = (double *) SuiteSparse_malloc (*asize, sizeof (double)) ;
1471 ok = ok && (Ax != NULL && Az != NULL) ;
1472 }
1473 else /* if (*mkind == 4) */
1474 {
1475 /* return A as merged-complex */
1476 Ax = (double *) SuiteSparse_malloc (*asize, 2*sizeof (double)) ;
1477 ok = ok && (Ax != NULL) ;
1478 }
1479
1480 /* ---------------------------------------------------------------------- */
1481 /* allocate workspace */
1482 /* ---------------------------------------------------------------------- */
1483
1484 cp = (Int *) SuiteSparse_malloc ((*ncol) + 1, sizeof (Int)) ;
1485 w = (Int *) SuiteSparse_malloc (MAX (*nrow, *ncol) + 1, sizeof (Int)) ;
1486 ok = ok && (cp != NULL && w != NULL) ;
1487
1488 if (!ok)
1489 {
1490 FREE_ALL ;
1491 return (RBIO_OUT_OF_MEMORY) ;
1492 }
1493
1494 /* ---------------------------------------------------------------------- */
1495 /* read the matrix */
1496 /* ---------------------------------------------------------------------- */
1497
1498 if (filename)
1499 {
1500 file = fopen (filename, "r") ;
1501 if (file == NULL)
1502 {
1503 FREE_ALL ;
1504 return (RBIO_FILE_IOERROR) ; /* cannot reopen file */
1505 }
1506 }
1507
1508 status = RB(read2) (file, *nrow, *ncol, nnz, *mkind, *skind, build_upper,
1509 Ap, Ai, Ax, Az, w, cp, s, SLEN) ;
1510
1511 if (filename) fclose (file) ;
1512
1513 if (status != 0)
1514 {
1515 FREE_ALL ;
1516 return (status) ;
1517 }
1518
1519 /* ---------------------------------------------------------------------- */
1520 /* free workspace */
1521 /* ---------------------------------------------------------------------- */
1522
1523 FREE_WORK ;
1524
1525 /* ---------------------------------------------------------------------- */
1526 /* prune or extract exact zeros */
1527 /* ---------------------------------------------------------------------- */
1528
1529 if (zero_handling == 2)
1530 {
1531 /* allocate the Z matrix */
1532 *znz = RB(zcount) (Ap [*ncol], *mkind, Ax, Az) ;
1533 Zp = (Int *) SuiteSparse_malloc ((*ncol) + 1, sizeof (Int)) ;
1534 Zi = (Int *) SuiteSparse_malloc (*znz, sizeof (Int)) ;
1535 if (Zp == NULL || Zi == NULL)
1536 {
1537 FREE_ALL ;
1538 return (RBIO_OUT_OF_MEMORY) ;
1539 }
1540 /* remove zeros from A and store their pattern in Z */
1541 RB(extract) (*ncol, *mkind, Ap, Ai, Ax, Az, Zp, Zi) ;
1542 }
1543 else if (zero_handling == 1)
1544 {
1545 /* remove zeros from A and discard them */
1546 *znz = RB(extract) (*ncol, *mkind, Ap, Ai, Ax, Az, NULL, NULL) ;
1547 }
1548 else
1549 {
1550 *znz = 0 ;
1551 }
1552
1553 /* ---------------------------------------------------------------------- */
1554 /* return results */
1555 /* ---------------------------------------------------------------------- */
1556
1557 if (p_Ap) *p_Ap = Ap ;
1558 if (p_Ai) *p_Ai = Ai ;
1559 if (p_Ax) *p_Ax = Ax ;
1560 if (p_Az) *p_Az = Az ;
1561 if (p_Zp) *p_Zp = Zp ;
1562 if (p_Zi) *p_Zi = Zi ;
1563 return (RBIO_OK) ;
1564 }
1565
1566
1567 /* -------------------------------------------------------------------------- */
1568 /* RBreadraw: read the raw contents of a Rutherford/Boeing file */
1569 /* -------------------------------------------------------------------------- */
1570
RB(readraw)1571 PUBLIC Int RB(readraw) /* 0: OK, < 0: error, > 0: warning */
1572 (
1573 /* input */
1574 char *filename, /* filename to read from */
1575
1576 /* output */
1577 char title [73],
1578 char key [9],
1579 char mtype [4], /* RUA, RSA, PUA, PSA, RRA, etc */
1580 Int *nrow, /* A is nrow-by-ncol */
1581 Int *ncol,
1582 Int *nnz, /* size of Ai */
1583 Int *nelnz,
1584 Int *mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
1585 Int *skind, /* R: -1, U: 0, S: 1, H: 2, Z: 3 */
1586 Int *fem, /* 0:__A, 1:__E */
1587 Int *xsize, /* size of Ax */
1588
1589 /* output: these are malloc'ed below and must be freed by the caller */
1590 Int **p_Ap, /* size ncol+1, column pointers of A */
1591 Int **p_Ai, /* size nnz, row indices of A */
1592 double **p_Ax /* size xsize, numerical values of A */
1593 )
1594 {
1595 FILE *file = NULL ; /* read from stdin if NULL */
1596 Int *Ap, *Ai ;
1597 double *Ax ;
1598 Int status ;
1599 int ok ;
1600 char s [SLEN+1], ptrfmt [21], indfmt [21], valfmt [21] ;
1601
1602 /* ---------------------------------------------------------------------- */
1603 /* check inputs */
1604 /* ---------------------------------------------------------------------- */
1605
1606 if (p_Ap) *p_Ap = NULL ;
1607 if (p_Ai) *p_Ai = NULL ;
1608 if (p_Ax) *p_Ax = NULL ;
1609
1610 if (!title || !key || !mtype || !nrow || !ncol || !nnz || !nelnz || !mkind
1611 || !skind || !fem || !xsize || !p_Ap || !p_Ai || !p_Ax)
1612 {
1613 /* one or more required arguments are missing (NULL) */
1614 return (RBIO_ARG_ERROR) ;
1615 }
1616
1617 /* ---------------------------------------------------------------------- */
1618 /* read the header */
1619 /* ---------------------------------------------------------------------- */
1620
1621 if (filename)
1622 {
1623 file = fopen (filename, "r") ;
1624 if (file == NULL)
1625 {
1626 return (RBIO_FILE_IOERROR) ; /* cannot open file */
1627 }
1628 }
1629
1630 status = RB(header) (file, title, key, mtype, nrow, ncol, nnz, nelnz,
1631 ptrfmt, indfmt, valfmt, mkind, skind, fem, s, SLEN) ;
1632
1633 /* close the file, so if there is an error the file is not left open */
1634 if (filename) fclose (file) ;
1635
1636 if (status != 0)
1637 {
1638 return (status) ;
1639 }
1640
1641 /* ---------------------------------------------------------------------- */
1642 /* allocate space for Ap, Ai, and Ax */
1643 /* ---------------------------------------------------------------------- */
1644
1645 Ap = (Int *) SuiteSparse_malloc ((*ncol) + 1, sizeof (Int)) ;
1646 Ai = (Int *) SuiteSparse_malloc (*nnz, sizeof (Int)) ;
1647 ok = (Ap != NULL && Ai != NULL) ;
1648
1649 if (*mkind == 1)
1650 {
1651 /* A is pattern-only */
1652 *xsize = 0 ;
1653 Ax = NULL ;
1654 }
1655 else
1656 {
1657 /* A has numerical values */
1658 *xsize = ((*fem) ? (*nelnz) : (*nnz)) * (((*mkind) == 2) ? 2 : 1) ;
1659 Ax = (double *) SuiteSparse_malloc (*xsize, sizeof (double)) ;
1660 ok = ok && (Ax != NULL) ;
1661 }
1662
1663 if (!ok)
1664 {
1665 FREE_RAW ;
1666 return (RBIO_OUT_OF_MEMORY) ;
1667 }
1668
1669 /* ---------------------------------------------------------------------- */
1670 /* read the matrix */
1671 /* ---------------------------------------------------------------------- */
1672
1673 if (filename)
1674 {
1675 file = fopen (filename, "r") ;
1676 if (file == NULL)
1677 {
1678 FREE_RAW ;
1679 return (RBIO_FILE_IOERROR) ; /* cannot reopen file */
1680 }
1681 RB(skipheader) (s, SLEN, file) ; /* skip past the header */
1682 }
1683
1684 if (!RB(iread) (file, (*ncol)+1, 1, Ap, s, SLEN))
1685 {
1686 FREE_RAW ;
1687 if (filename) fclose (file) ;
1688 return (RBIO_CP_IOERROR) ; /* I/O error reading column pointers */
1689 }
1690
1691 if (!RB(iread) (file, *nnz, 1, Ai, s, SLEN))
1692 {
1693 FREE_RAW ;
1694 if (filename) fclose (file) ;
1695 return (RBIO_ROW_IOERROR) ; /* I/O error reading row indices */
1696 }
1697
1698 if (*mkind != 1)
1699 {
1700 if (!RB(xread) (file, *xsize, 0, Ax, NULL, s, SLEN))
1701 {
1702 FREE_RAW ;
1703 if (filename) fclose (file) ;
1704 return (RBIO_VALUE_IOERROR) ; /* I/O error reading values */
1705 }
1706 }
1707
1708 /* ---------------------------------------------------------------------- */
1709 /* return results */
1710 /* ---------------------------------------------------------------------- */
1711
1712 if (p_Ap) *p_Ap = Ap ;
1713 if (p_Ai) *p_Ai = Ai ;
1714 if (p_Ax) *p_Ax = Ax ;
1715 if (filename) fclose (file) ;
1716 return (RBIO_OK) ;
1717 }
1718
1719
1720 /* -------------------------------------------------------------------------- */
1721 /* RBfix_mkind_in: adjust mkind_in, based on presence of Ax and Az */
1722 /* -------------------------------------------------------------------------- */
1723
RB(fix_mkind_in)1724 PRIVATE Int RB(fix_mkind_in) /* return revised mkind */
1725 (
1726 Int mkind_in, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
1727 double *Ax,
1728 double *Az
1729 )
1730 {
1731 if (!Ax)
1732 {
1733 /* matrix must be 1:pattern */
1734 mkind_in = 1 ;
1735 }
1736 if (mkind_in == 2 && !Az)
1737 {
1738 /* matrix must be 4:merged-complex */
1739 mkind_in = 4 ;
1740 }
1741 return (mkind_in) ;
1742 }
1743
1744
1745 /* -------------------------------------------------------------------------- */
1746 /* RBwrite */
1747 /* -------------------------------------------------------------------------- */
1748
RB(write)1749 PUBLIC Int RB(write) /* 0:OK, < 0: error, > 0: warning */
1750 (
1751 /* input */
1752 char *filename, /* filename to write to (stdout if NULL) */
1753 char *title, /* title (72 char max), may be NULL */
1754 char *key, /* key (8 char max), may be NULL */
1755 Int nrow, /* A is nrow-by-ncol */
1756 Int ncol,
1757 Int *Ap, /* size ncol+1, column pointers */
1758 Int *Ai, /* size anz=Ap[ncol], row indices (sorted) */
1759 double *Ax, /* size anz or 2*anz, numerical values (binary if NULL) */
1760 double *Az, /* size anz, imaginary part (real if NULL) */
1761 Int *Zp, /* size ncol+1, column pointers for Z (or NULL) */
1762 Int *Zi, /* size znz=Zp[ncol], row indices for Z (or NULL) */
1763 Int mkind_in, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
1764
1765 /* output */
1766 char mtype [4] /* matrix type (RUA, RSA, etc), may be NULL */
1767 )
1768 {
1769 double xmin, xmax, zmin, zmax ;
1770 Int *w, *cp ;
1771 FILE *file = NULL ; /* write to stdout if NULL */
1772 Int mkind, skind, zmkind, zskind, nnz2, vals, valn, indn, ptrn, valcrd,
1773 indcrd, ptrcrd, totcrd, anz, is_int, fmt, nbuf, j, njumbled, nzeros,
1774 znz, asize, status, nelnz = 0 ;
1775 int ok ;
1776 char zmtype [4], indfmt [21], indcfm [21], valfmt [21], valcfm [21],
1777 ptrfmt [21], ptrcfm [21], mtype2 [4] ;
1778
1779 /* ---------------------------------------------------------------------- */
1780 /* check inputs */
1781 /* ---------------------------------------------------------------------- */
1782
1783 if (!title || !key || !Ap || !Ax)
1784 {
1785 /* one or more required arguments are missing (NULL) */
1786 return (RBIO_ARG_ERROR) ;
1787 }
1788
1789 mkind_in = RB(fix_mkind_in) (mkind_in, Ax, Az) ;
1790
1791 anz = Ap ? (Ap [MAX (ncol,0)]) : 1 ;
1792 status = RB(ok) (nrow, ncol, anz, Ap, Ai, Ax, Az, NULL, mkind_in,
1793 &njumbled, &nzeros) ;
1794 if (status != RBIO_OK)
1795 {
1796 /* A matrix is corrupted */
1797 return (status) ;
1798 }
1799
1800 if (Zp != NULL)
1801 {
1802 znz = Zp [MAX (ncol,0)] ;
1803 status = RB(ok) (nrow, ncol, znz, Zp, Zi, NULL, NULL, NULL, 3,
1804 &njumbled, &nzeros) ;
1805 if (status != RBIO_OK)
1806 {
1807 /* Z matrix is corrupted */
1808 return (status) ;
1809 }
1810 }
1811
1812 if (mtype == NULL)
1813 {
1814 /* use internal array for mtype; do not return to caller */
1815 mtype = mtype2 ;
1816 }
1817
1818 /* ---------------------------------------------------------------------- */
1819 /* clear the format strings */
1820 /* ---------------------------------------------------------------------- */
1821
1822 RB(fill) (valfmt, 20, ' ') ;
1823 RB(fill) (valcfm, 20, ' ') ;
1824 RB(fill) (indfmt, 20, ' ') ;
1825 RB(fill) (indcfm, 20, ' ') ;
1826 RB(fill) (ptrfmt, 20, ' ') ;
1827 RB(fill) (ptrcfm, 20, ' ') ;
1828 indn = 0 ;
1829 valn = 0 ;
1830 ptrn = 0 ;
1831
1832 /* ---------------------------------------------------------------------- */
1833 /* allocate workspace */
1834 /* ---------------------------------------------------------------------- */
1835
1836 w = SuiteSparse_malloc (MAX (nrow, ncol) + 1, sizeof (Int)) ;
1837 cp = SuiteSparse_malloc (ncol + 1, sizeof (Int)) ;
1838 if (cp == NULL || w == NULL)
1839 {
1840 FREE_WORK ;
1841 return (RBIO_OUT_OF_MEMORY) ;
1842 }
1843
1844 /* ---------------------------------------------------------------------- */
1845 /* determine the matrix type (RSA, RUA, etc) of A+Z */
1846 /* ---------------------------------------------------------------------- */
1847
1848 RB(kind) (nrow, ncol, Ap, Ai, Ax, Az, mkind_in, &mkind, &skind, mtype,
1849 &xmin, &xmax, cp) ;
1850
1851 /* now use mkind instead of mkind_in */
1852
1853 if (Zp != NULL && Zp [ncol] == 0)
1854 {
1855 /* Z has no entries; ignore it */
1856 Zp = NULL ;
1857 Zi = NULL ;
1858 }
1859
1860 if (Zp != NULL)
1861 {
1862 /* determine if Z is symmetric or not */
1863 RB(kind) (nrow, ncol, Zp, Zi, NULL, NULL, 3, &zmkind, &zskind, zmtype,
1864 &zmin, &zmax, cp) ;
1865 if (zskind == 0)
1866 {
1867 /* Z is square and unsymmetric; force A unsymmetric too */
1868 mtype [1] = 'u' ;
1869 skind = 0 ;
1870 }
1871 }
1872
1873 /* ---------------------------------------------------------------------- */
1874 /* determine the required precision for the numerical values */
1875 /* ---------------------------------------------------------------------- */
1876
1877 asize = ((mkind == 4) ? 2 : 1) * anz ;
1878 is_int = (mkind == 3) ;
1879 if (mkind != 1)
1880 {
1881 /* A is real, split-complex, integer or merged-complex: check Ax */
1882 fmt = RB(format) (asize, Ax, is_int, xmin, xmax, 0,
1883 valfmt, valcfm, &valn) ;
1884 }
1885 if (mkind == 2)
1886 {
1887 /* A is split-complex: check Az */
1888 fmt = RB(format) (anz, Az, FALSE, 0, 0, fmt, valfmt, valcfm, &valn) ;
1889 }
1890
1891 /* ---------------------------------------------------------------------- */
1892 /* determine the number of entries in the matrix A+Z */
1893 /* ---------------------------------------------------------------------- */
1894
1895 /* task 1 does not write to the file */
1896 ok = RB(writeTask) (NULL, 1, nrow, ncol, mkind, skind, Ap, Ai, Ax, Az, Zp,
1897 Zi, indcfm, indn, valcfm, valn, &nnz2, w, cp) ;
1898 if (nnz2 <= 0)
1899 {
1900 FREE_WORK ;
1901 return (RBIO_DIM_INVALID) ; /* matrix has no entries to print */
1902 }
1903
1904 /* determine pointer format. ncol+1 integers, in range 1 to nnz2+1 */
1905 RB(iformat) (1, nnz2+1, ptrfmt, ptrcfm, &ptrn) ;
1906 ptrcrd = RB(cards) (ncol+1, ptrn) ;
1907
1908 /* determine row index format. nnz2 integers, in range 1 to nrow */
1909 RB(iformat) (1, nrow, indfmt, indcfm, &indn) ;
1910 indcrd = RB(cards) (nnz2, indn) ;
1911
1912 /* determine how many lines for the numerical values */
1913 if (mkind == 0 || mkind == 3)
1914 {
1915 /* real or integer */
1916 vals = 1 ;
1917 }
1918 else if (mkind == 1)
1919 {
1920 /* pattern */
1921 vals = 0 ;
1922 }
1923 else
1924 {
1925 /* complex (merged or split) */
1926 vals = 2 ;
1927 }
1928 valcrd = RB(cards) (vals*nnz2, valn) ;
1929
1930 /* ---------------------------------------------------------------------- */
1931 /* determine total number of cards */
1932 /* ---------------------------------------------------------------------- */
1933
1934 totcrd = ptrcrd + indcrd + valcrd ;
1935
1936 /* ---------------------------------------------------------------------- */
1937 /* open the file */
1938 /* ---------------------------------------------------------------------- */
1939
1940 if (filename)
1941 {
1942 file = fopen (filename, "w") ;
1943 if (file == NULL)
1944 {
1945 FREE_WORK ;
1946 return (RBIO_FILE_IOERROR) ; /* cannot open file */
1947 }
1948 }
1949
1950 /* ---------------------------------------------------------------------- */
1951 /* write the header */
1952 /* ---------------------------------------------------------------------- */
1953
1954 ok = fprintf (file, "%-71.71s|%-8.8s\n",
1955 title ? title : "", key ? key : "") > 0;
1956 ok = ok && fprintf (file, "%14" IDD "%14" IDD "%14" IDD "%14" IDD "\n",
1957 totcrd, ptrcrd, indcrd, valcrd) > 0 ;
1958 ok = ok && fprintf (file,
1959 "%3s %14" IDD "%14" IDD "%14" IDD "%14" IDD "\n",
1960 mtype, nrow, ncol, nnz2, nelnz) > 0 ;
1961 ok = ok && fprintf (file, "%.16s%.16s%.20s\n", ptrfmt, indfmt, valfmt) > 0 ;
1962 if (!ok)
1963 {
1964 /* file I/O error */
1965 FREE_WORK ;
1966 if (filename) fclose (file) ;
1967 return (RBIO_HEADER_IOERROR) ;
1968 }
1969
1970 /* ---------------------------------------------------------------------- */
1971 /* write the column pointers (convert to 1-based) */
1972 /* ---------------------------------------------------------------------- */
1973
1974 nbuf = 0 ;
1975 for (j = 0 ; ok && j <= ncol ; j++)
1976 {
1977 ok = RB(iprint) (file, ptrcfm, 1+ cp [j], ptrn, &nbuf) ;
1978 }
1979 ok = ok && fprintf (file, "\n") > 0 ;
1980 if (!ok)
1981 {
1982 /* file I/O error */
1983 FREE_WORK ;
1984 if (filename) fclose (file) ;
1985 return (RBIO_CP_IOERROR) ;
1986 }
1987
1988 /* ---------------------------------------------------------------------- */
1989 /* write the row indices (convert to 1-based) */
1990 /* ---------------------------------------------------------------------- */
1991
1992 ok = RB(writeTask) (file, 2, nrow, ncol, mkind, skind, Ap, Ai, Ax, Az,
1993 Zp, Zi, indcfm, indn, valcfm, valn, &nnz2, w, cp) ;
1994 if (!ok)
1995 {
1996 /* file I/O error */
1997 FREE_WORK ;
1998 if (filename) fclose (file) ;
1999 return (RBIO_ROW_IOERROR) ;
2000 }
2001
2002 /* ---------------------------------------------------------------------- */
2003 /* write the numerical values */
2004 /* ---------------------------------------------------------------------- */
2005
2006 if (mkind != 1)
2007 {
2008 ok = RB(writeTask) (file, 3, nrow, ncol, mkind, skind, Ap, Ai,
2009 Ax, Az, Zp, Zi, indcfm, indn, valcfm, valn, &nnz2, w, cp) ;
2010 }
2011
2012 /* ---------------------------------------------------------------------- */
2013 /* free workspace and close the file */
2014 /* ---------------------------------------------------------------------- */
2015
2016 FREE_WORK ;
2017 if (filename) fclose (file) ;
2018 return (ok ? RBIO_OK : RBIO_VALUE_IOERROR) ;
2019 }
2020
2021
2022 /* -------------------------------------------------------------------------- */
2023 /* RBkind: determine the type of a sparse matrix */
2024 /* -------------------------------------------------------------------------- */
2025
RB(kind)2026 PUBLIC Int RB(kind) /* 0: OK, < 0: error, > 0: warning */
2027 (
2028 /* input */
2029 Int nrow, /* A is nrow-by-ncol */
2030 Int ncol,
2031 Int *Ap, /* Ap [0...ncol]: column pointers */
2032 Int *Ai, /* Ai [0...nnz-1]: row indices */
2033 double *Ax, /* Ax [0...nnz-1]: real values. Az holds imaginary part */
2034 double *Az, /* if real, Az is NULL. if complex, Az is non-NULL */
2035 Int mkind_in, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
2036
2037 /* output */
2038 Int *mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
2039 Int *skind, /* r: -1 (rectangular), u: 0 (unsymmetric), s: 1 symmetric,
2040 h: 2 (Hermitian), z: 3 (skew symmetric) */
2041 char mtype [4], /* rua, psa, rra, cha, etc */
2042 double *xmin, /* smallest value */
2043 double *xmax, /* largest value */
2044
2045 /* workspace: allocated internally if NULL */
2046 Int *cp /* workspace of size ncol+1, undefined on input and output*/
2047 )
2048 {
2049 Int nnz, is_h, is_z, is_s, k, p, i, j, pt, get_workspace ;
2050 Int *w = NULL ;
2051 double aij_real, aij_imag, aji_real, aji_imag ;
2052
2053 /* ---------------------------------------------------------------------- */
2054 /* check inputs */
2055 /* ---------------------------------------------------------------------- */
2056
2057 if (!Ap || !Ai || !mkind || !skind || !mtype || !xmin || !xmax
2058 || ncol < 0 || nrow < 0)
2059 {
2060 /* one or more required arguments are missing (NULL) or invalid */
2061 return (RBIO_ARG_ERROR) ;
2062 }
2063
2064 /* ---------------------------------------------------------------------- */
2065 /* allocate workspace, if needed */
2066 /* ---------------------------------------------------------------------- */
2067
2068 get_workspace = (cp == NULL) ;
2069 if (get_workspace)
2070 {
2071 cp = (Int *) SuiteSparse_malloc (ncol + 1, sizeof (Int)) ;
2072 }
2073 if (cp == NULL)
2074 {
2075 return (RBIO_OUT_OF_MEMORY) ;
2076 }
2077
2078 /* ---------------------------------------------------------------------- */
2079 /* determine numeric type (I*A, R*A, P*A, C*A) */
2080 /* ---------------------------------------------------------------------- */
2081
2082 mkind_in = RB(fix_mkind_in) (mkind_in, Ax, Az) ;
2083
2084 mtype [3] = '\0' ;
2085 nnz = Ap [ncol] ;
2086 *xmin = 0 ;
2087 *xmax = 0 ;
2088
2089 if (mkind_in == 2 || mkind_in == 4)
2090 {
2091 /* A is a complex matrix (C*A), split or merged */
2092 mtype [0] = 'c' ;
2093 *mkind = mkind_in ;
2094 }
2095 else if (mkind_in == 1)
2096 {
2097 /* A is a pattern-only matrix (P*A) */
2098 mtype [0] = 'p' ;
2099 *mkind = 1 ;
2100 }
2101 else
2102 {
2103 /* The input matrix A is said to be 0:real or 3:integer, */
2104 /* Ax is not NULL */
2105
2106 /* select P** format if all entries are equal to 1 */
2107 /* select I** format if all entries are integer and */
2108 /* between -99,999,999 and +999,999,999 */
2109
2110 Int is_p = TRUE ;
2111 Int is_int = TRUE ;
2112 *xmin = Ax [0] ;
2113 *xmax = Ax [0] ;
2114
2115 for (p = 0 ; (is_p || is_int) && p < nnz ; p++)
2116 {
2117 if (is_p)
2118 {
2119 if (Ax [p] != 1)
2120 {
2121 is_p = FALSE ;
2122 }
2123 }
2124 if (is_int)
2125 {
2126 double x = Ax [p] ;
2127 k = (Int) x ;
2128 *xmin = MIN (x, *xmin) ;
2129 *xmax = MAX (x, *xmax) ;
2130 if (((double) k) != x)
2131 {
2132 /* the entry is not an integer */
2133 is_int = FALSE ;
2134 }
2135 if (x < -99999999 || x >= 999999999)
2136 {
2137 /* use real format for really big integers */
2138 is_int = FALSE ;
2139 }
2140 }
2141 }
2142
2143 if (is_p)
2144 {
2145 /* pattern-only matrix (P*A) */
2146 mtype [0] = 'p' ;
2147 *mkind = 1 ;
2148 }
2149 else if (is_int)
2150 {
2151 /* integer matrix (I*A) */
2152 mtype [0] = 'i' ;
2153 *mkind = 3 ;
2154 }
2155 else
2156 {
2157 /* real matrix (R*A) */
2158 mtype [0] = 'r' ;
2159 *mkind = 0 ;
2160 }
2161 }
2162
2163 /* only assembled matrices are handled */
2164 mtype [2] = 'a' ;
2165
2166 /* ---------------------------------------------------------------------- */
2167 /* determine symmetry (*RA, *UA, *SA, *HA, *ZA) */
2168 /* ---------------------------------------------------------------------- */
2169
2170 /* Note that A must have sorted columns for this method to work. */
2171
2172 if (nrow != ncol)
2173 {
2174 /* rectangular matrix (*RA), no need to check values or pattern */
2175 mtype [1] = 'r' ;
2176 *skind = -1 ;
2177 if (get_workspace) FREE_WORK ;
2178 return (RBIO_OK) ;
2179 }
2180
2181 /* if complex, the matrix is Hermitian until proven otherwise */
2182 is_h = (*mkind == 2 || *mkind == 4) ;
2183
2184 /* the matrix is symmetric until proven otherwise */
2185 is_s = TRUE ;
2186
2187 /* a non-pattern matrix is skew symmetric until proven otherwise */
2188 is_z = (*mkind != 1) ;
2189
2190 /* if this method returns early, the matrix is unsymmetric */
2191 mtype [1] = 'u' ;
2192 *skind = 0 ;
2193
2194 /* initialize the munch pointers (cp) */
2195 for (j = 0 ; j <= ncol ; j++)
2196 {
2197 cp [j] = Ap [j] ;
2198 }
2199
2200 for (j = 0 ; j < ncol ; j++)
2201 {
2202
2203 /* consider all entries not yet munched in column j */
2204 for (p = cp [j] ; p < Ap [j+1] ; p++)
2205 {
2206 /* get the row index of A(i,j) */
2207 i = Ai [p] ;
2208
2209 if (i < j)
2210 {
2211 /* entry A(i,j) is unmatched, matrix is unsymmetric */
2212 if (get_workspace) FREE_WORK ;
2213 return (RBIO_OK) ;
2214 }
2215
2216 /* get the A(j,i) entry, if it exists, and munch it */
2217 pt = cp [i]++ ;
2218
2219 if (pt >= Ap [i+1] || Ai [pt] != j)
2220 {
2221 /* entry A(j,i) doesn't exist, matrix is unsymmetric */
2222 if (get_workspace) FREE_WORK ;
2223 return (RBIO_OK) ;
2224 }
2225
2226 /* A(j,i) exists; check its value with A(i,j) */
2227 RB(get_entry) (*mkind, Ax, Az, p, &aij_real, &aij_imag) ;
2228 RB(get_entry) (*mkind, Ax, Az, pt, &aji_real, &aji_imag) ;
2229
2230 if (aij_real != aji_real || aij_imag != aji_imag)
2231 {
2232 is_s = FALSE ; /* the matrix cannot be *SA */
2233 }
2234 if (aij_real != -aji_real || aij_imag != -aji_imag)
2235 {
2236 is_z = FALSE ; /* the matrix cannot be *ZA */
2237 }
2238 if (aij_real != aji_real || aij_imag != -aji_imag)
2239 {
2240 is_h = FALSE ; /* the matrix cannot be *HA */
2241 }
2242
2243 if (! (is_s || is_z || is_h))
2244 {
2245 /* matrix is unsymmetric; terminate the test */
2246 if (get_workspace) FREE_WORK ;
2247 return (RBIO_OK) ;
2248 }
2249 }
2250 }
2251
2252 /* ---------------------------------------------------------------------- */
2253 /* return the symmetry */
2254 /* ---------------------------------------------------------------------- */
2255
2256 if (is_h)
2257 {
2258 /* Hermitian matrix (*HA) */
2259 mtype [1] = 'h' ;
2260 *skind = 2 ;
2261 }
2262 else if (is_s)
2263 {
2264 /* symmetric matrix (*SA) */
2265 mtype [1] = 's' ;
2266 *skind = 1 ;
2267 }
2268 else if (is_z)
2269 {
2270 /* skew symmetric matrix (*ZA) */
2271 mtype [1] = 'z' ;
2272 *skind = 3 ;
2273 }
2274
2275 /* ---------------------------------------------------------------------- */
2276 /* free workspace, if allocated */
2277 /* ---------------------------------------------------------------------- */
2278
2279 if (get_workspace) FREE_WORK ;
2280 return (RBIO_OK) ;
2281 }
2282
2283
2284 /* -------------------------------------------------------------------------- */
2285 /* RBwrite data formats */
2286 /* -------------------------------------------------------------------------- */
2287
2288 #define NFORMAT 19
2289
2290 /* Fortran formats */
2291 static char *F_format [NFORMAT] = {
2292 "(8E9.1) ",
2293 "(8E10.2) ",
2294 "(7E11.3) ",
2295 "(6E12.4) ",
2296 "(6E13.5) ",
2297 "(5E14.6) ",
2298 "(5E15.7) ",
2299 "(5E16.8) ",
2300 "(4E17.9) ",
2301 "(4E18.10) ",
2302 "(4E19.11) ",
2303 "(4E20.12) ",
2304 "(3E21.13) ",
2305 "(3E22.14) ",
2306 "(3E23.15) ",
2307 "(3E24.16) ",
2308 "(3E25.17) ",
2309 "(3E26.18) ",
2310 "(2E30.18E3) " } ;
2311
2312 /* corresponding C formats to use */
2313 static char *C_format [NFORMAT] = {
2314 "%9.1E",
2315 "%10.2E",
2316 "%11.3E",
2317 "%12.4E",
2318 "%13.5E",
2319 "%14.6E",
2320 "%15.7E",
2321 "%16.8E",
2322 "%17.9E",
2323 "%18.10E",
2324 "%19.11E",
2325 "%20.12E",
2326 "%21.13E",
2327 "%22.14E",
2328 "%23.15E",
2329 "%24.16E",
2330 "%25.17E",
2331 "%26.18E",
2332 "%30.18E" } ;
2333
2334 /* Number of entries per line for each of the formats */
2335 static Int entries_per_line [NFORMAT] = {
2336 8,
2337 8,
2338 7,
2339 6,
2340 6,
2341 5,
2342 5,
2343 5,
2344 4,
2345 4,
2346 4,
2347 4,
2348 3,
2349 3,
2350 3,
2351 3,
2352 3,
2353 3,
2354 2} ;
2355
2356
2357 /* -------------------------------------------------------------------------- */
2358 /* RBformat: determine the format required for an array of values */
2359 /* -------------------------------------------------------------------------- */
2360
2361 /*
2362 This function ensures that a sufficiently wide format is used that
2363 accurately represent the data. It also ensures that when printed,
2364 the numerical values all have at least one blank space between them.
2365 */
2366
RB(format)2367 PRIVATE Int RB(format) /* return format to use (index in F_, C_format) */
2368 (
2369 /* input */
2370 Int nnz, /* number of nonzeros */
2371 double *x, /* of size nnz */
2372 Int is_int, /* true if integer format is to be used */
2373 double xmin, /* minimum value of x */
2374 double xmax, /* maximum value of x */
2375 Int fmt, /* initial format to use (index into F_format, ...) */
2376
2377 /* output */
2378 char valfmt [21], /* Fortran format to use */
2379 char valcfm [21], /* C format to use */
2380 Int *valn /* number of entries per line */
2381 )
2382 {
2383 Int i ;
2384 double a, b ;
2385 char s [1024] ;
2386
2387 if (is_int)
2388 {
2389
2390 /* ------------------------------------------------------------------ */
2391 /* use an integer format */
2392 /* ------------------------------------------------------------------ */
2393
2394 RB(iformat) (xmin, xmax, valfmt, valcfm, valn) ;
2395 return (-1) ;
2396
2397 }
2398 else
2399 {
2400
2401 /* ------------------------------------------------------------------ */
2402 /* find the required precision for a real or complex matrix */
2403 /* ------------------------------------------------------------------ */
2404
2405 fmt = 0 ;
2406 for (i = 0 ; i < nnz ; i++)
2407 {
2408
2409 /* determine if the matrix has huge values, tiny values, or NaN's */
2410 a = ABS (x [i]) ;
2411 if (a != 0)
2412 {
2413 if (ISNAN (a) || a < 1e-90 || a > 1e90)
2414 {
2415 fmt = NFORMAT-1 ;
2416 break ;
2417 }
2418 }
2419
2420 a = x ? x [i] : 1 ;
2421 for ( ; fmt < NFORMAT-1 ; fmt++)
2422 {
2423 /* write the value to a string, read back in, and check,
2424 * using the kth format */
2425 sprintf (s, C_format [fmt], a) ;
2426 b = 0 ;
2427 sscanf (s, "%lg", &b) ;
2428 if (s [0] == ' ' && a == b)
2429 {
2430 /* success, use this format (or wider) for all numbers */
2431 break ;
2432 }
2433 }
2434 }
2435
2436 strncpy (valfmt, F_format [fmt], 21) ;
2437 strncpy (valcfm, C_format [fmt], 21) ;
2438 *valn = entries_per_line [fmt] ;
2439 return (fmt) ;
2440 }
2441 }
2442
2443
2444 /* -------------------------------------------------------------------------- */
2445 /* RBwriteTask: write portions of the matrix to the file */
2446 /* -------------------------------------------------------------------------- */
2447
2448 /*
2449 task 0: just count the total number of entries in the matrix. No file I/O.
2450 task 1: do task 0, and also construct w and cp. No file I/O.
2451 task 2: write the row indices
2452 task 3: write the numerical values
2453 */
2454
RB(writeTask)2455 PRIVATE Int RB(writeTask) /* returns TRUE if successful, FALSE on failure */
2456 (
2457 /* input */
2458 FILE *file, /* file to print to (already open) */
2459 Int task, /* 0 to 3 (see above) */
2460 Int nrow, /* A is nrow-by-ncol */
2461 Int ncol,
2462 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
2463 Int skind, /* -1:rect, 0:unsym, 1:sym, 2:hermitian, 3:skew */
2464 Int *Ap, /* size ncol+1, column pointers */
2465 Int *Ai, /* size anz=Ap[ncol], row indices */
2466 double *Ax, /* size anz, real values */
2467 double *Az, /* size anz, imaginary part (may be NULL) */
2468 Int *Zp, /* size ncol+1, column pointers for Z (may be NULL) */
2469 Int *Zi, /* size Zp[ncol], row indices for Z */
2470 char *indcfm, /* C format for indices */
2471 Int indn, /* # of indices per line */
2472 char *valcfm, /* C format for values */
2473 Int valn, /* # of values per line */
2474
2475 /* output */
2476 Int *nnz, /* number of entries that will be printed to the file */
2477
2478 /* workspace */
2479 Int *w, /* size MAX(nrow,ncol)+1 */
2480 Int *cp /* size ncol+1 */
2481 )
2482 {
2483 double xr, xz ;
2484 Int j, pa, pz, paend, pzend, ia, iz, i, nbuf, ok ;
2485
2486 /* ---------------------------------------------------------------------- */
2487 /* clear the nonzero counts */
2488 /* ---------------------------------------------------------------------- */
2489
2490 *nnz = 0 ;
2491 for (j = 0 ; j < ncol ; j++)
2492 {
2493 w [j] = 0 ;
2494 }
2495
2496 /* ---------------------------------------------------------------------- */
2497 /* print, or count, each column */
2498 /* ---------------------------------------------------------------------- */
2499
2500 nbuf = 0 ; /* number of characters in current line */
2501 ok = TRUE ;
2502 for (j = 0 ; ok && j < ncol ; j++)
2503 {
2504
2505 /* find the set union of A (:,j) and Z (:,j) */
2506 pa = Ap [j] ;
2507 pz = Zp ? Zp [j] : 0 ;
2508 paend = Ap [j+1] ;
2509 pzend = Zp ? Zp [j+1] : 0 ;
2510
2511 /* repeat while entries still exist in A(:,j) or Z(:,j) */
2512 while (ok)
2513 {
2514 /* get the next entry from A(:,j) */
2515 ia = (pa < paend) ? Ai [pa] : nrow ;
2516
2517 /* get the next entry from Z(:,j) */
2518 iz = (pz < pzend) ? Zi [pz] : nrow ;
2519
2520 /* exit loop if neither entry is present */
2521 if (ia >= nrow && iz >= nrow) break ;
2522
2523 if (ia < iz)
2524 {
2525 /* get A (i,j) */
2526 i = ia ;
2527 RB(get_entry) (mkind, Ax, Az, pa, &xr, &xz) ;
2528 pa++ ;
2529 }
2530 else if (iz < ia)
2531 {
2532 /* get Z (i,j) */
2533 i = iz ;
2534 xr = 0 ;
2535 xz = 0 ;
2536 pz++ ;
2537 }
2538 else
2539 {
2540 /* get A (i,j), and delete its matched Z(i,j) */
2541 i = ia ;
2542 RB(get_entry) (mkind, Ax, Az, pa, &xr, &xz) ;
2543 pa++ ;
2544 pz++ ;
2545 }
2546
2547 if (skind <= 0 || i >= j)
2548 {
2549 /* consider the (i,j) entry with value (xr,xz) */
2550 (*nnz)++ ;
2551 if (task == 1)
2552 {
2553 /* only determining nonzero counts */
2554 w [j]++ ;
2555 }
2556 else if (task == 2)
2557 {
2558 /* printing the row indices (convert to 1-based) */
2559 ok = RB(iprint) (file, indcfm, 1+ i, indn, &nbuf) ;
2560 }
2561 else if (task == 3)
2562 {
2563 /* printing the numerical values */
2564 ok = RB(xprint) (file, valcfm, xr, valn, mkind, &nbuf) ;
2565 if (ok && (mkind == 2 || mkind == 4))
2566 {
2567 ok = RB(xprint) (file, valcfm, xz, valn, mkind, &nbuf) ;
2568 }
2569 }
2570 }
2571 }
2572 }
2573
2574 /* ---------------------------------------------------------------------- */
2575 /* determine the new column pointers, or finish printing */
2576 /* ---------------------------------------------------------------------- */
2577
2578 if (ok)
2579 {
2580 if (task == 1)
2581 {
2582 cp [0] = 0 ;
2583 for (j = 0 ; j < ncol ; j++)
2584 {
2585 cp [j+1] = cp [j] + w [j] ;
2586 }
2587 }
2588 else
2589 {
2590 ok = (fprintf (file ? file : stdout, "\n") > 0) ;
2591 }
2592 }
2593 return (ok) ;
2594 }
2595
2596
2597 /* -------------------------------------------------------------------------- */
2598 /* RBiprint: print one integer value to the file */
2599 /* -------------------------------------------------------------------------- */
2600
RB(iprint)2601 PRIVATE Int RB(iprint) /* returns TRUE if OK, FALSE otherwise */
2602 (
2603 /* input */
2604 FILE *file, /* which file to write to */
2605 char *indcfm, /* C format to use */
2606 Int i, /* value to write */
2607 Int indn, /* number of entries to write per line */
2608
2609 /* input/output */
2610 Int *nbuf /* number of entries written to current line */
2611 )
2612 {
2613 Int ok = TRUE ;
2614 if (!file) file = stdout ; /* file defaults to stdout */
2615 if (*nbuf >= indn)
2616 {
2617 *nbuf = 0 ;
2618 ok = (fprintf (file, "\n") > 0) ;
2619 }
2620 ok = ok && (fprintf (file, indcfm, i) > 0) ;
2621 (*nbuf)++ ;
2622 return (ok) ;
2623 }
2624
2625
2626 /* -------------------------------------------------------------------------- */
2627 /* RBxprint: print one real value to the file */
2628 /* -------------------------------------------------------------------------- */
2629
RB(xprint)2630 PRIVATE Int RB(xprint) /* returns TRUE if OK, FALSE otherwise */
2631 (
2632 /* input */
2633 FILE *file, /* which file to write to */
2634 char *valcfm, /* C format to use */
2635 double x, /* value to write */
2636 Int valn, /* number of entries to write per line */
2637 Int mkind, /* 0:R, 1:P: 2:Csplit, 3:I, 4:Cmerged */
2638
2639 /* input/output */
2640 Int *nbuf /* number of entries written to current line */
2641 )
2642 {
2643 Int ok = TRUE ;
2644 if (!file) file = stdout ; /* file defaults to stdout */
2645 if (mkind == 3)
2646 {
2647 /* write out the value as an integer */
2648 ok = RB(iprint) (file, valcfm, (Int) x, valn, nbuf) ;
2649 }
2650 else
2651 {
2652 /* write out the value as a real */
2653 if (*nbuf >= valn)
2654 {
2655 *nbuf = 0 ;
2656 ok = (fprintf (file, "\n") > 0) ;
2657 }
2658 ok = ok && (fprintf (file, valcfm, x) > 0) ;
2659 (*nbuf)++ ;
2660 }
2661 return (ok) ;
2662 }
2663
2664
2665 /* -------------------------------------------------------------------------- */
2666 /* RBiformat: determine format for printing an integer */
2667 /* -------------------------------------------------------------------------- */
2668
RB(iformat)2669 PRIVATE void RB(iformat)
2670 (
2671 /* input */
2672 double xmin, /* smallest integer to print */
2673 double xmax, /* largest integer to print */
2674
2675 /* output */
2676 char indfmt [21], /* Fortran format to use */
2677 char indcfm [21], /* C format to use */
2678 Int *indn /* number of entries per line */
2679 )
2680 {
2681 if (xmin >= 0 && xmax <= 9.)
2682 {
2683 strncpy (indfmt, "(40I2) ", 21) ;
2684 strncpy (indcfm, "%2" IDD , 21) ;
2685 *indn = 40 ;
2686 }
2687 else if (xmin >= -9. && xmax <= 99.)
2688 {
2689 strncpy (indfmt, "(26I3) ", 21) ;
2690 strncpy (indcfm, "%3" IDD , 21) ;
2691 *indn = 26 ;
2692 }
2693 else if (xmin >= -99. && xmax <= 999.)
2694 {
2695 strncpy (indfmt, "(20I4) ", 21) ;
2696 strncpy (indcfm, "%4" IDD , 21) ;
2697 *indn = 20 ;
2698 }
2699 else if (xmin >= -999. && xmax <= 9999.)
2700 {
2701 strncpy (indfmt, "(16I5) ", 21) ;
2702 strncpy (indcfm, "%5" IDD , 21) ;
2703 *indn = 16 ;
2704 }
2705 else if (xmin >= -9999. && xmax <= 99999.)
2706 {
2707 strncpy (indfmt, "(13I6) ", 21) ;
2708 strncpy (indcfm, "%6" IDD , 21) ;
2709 *indn = 13 ;
2710 }
2711 else if (xmin >= -99999. && xmax <= 999999.)
2712 {
2713 strncpy (indfmt, "(11I7) ", 21) ;
2714 strncpy (indcfm, "%7" IDD , 21) ;
2715 *indn = 11 ;
2716 }
2717 else if (xmin >= -999999. && xmax <= 9999999.)
2718 {
2719 strncpy (indfmt, "(10I8) ", 21) ;
2720 strncpy (indcfm, "%8" IDD , 21) ;
2721 *indn = 10 ;
2722 }
2723 else if (xmin >= -9999999. && xmax <= 99999999.)
2724 {
2725 strncpy (indfmt, "(8I9) ", 21) ;
2726 strncpy (indcfm, "%9" IDD , 21) ;
2727 *indn = 8 ;
2728 }
2729 else if (xmin >= -99999999. && xmax <= 999999999.)
2730 {
2731 strncpy (indfmt, "(8I10) ", 21) ;
2732 strncpy (indcfm, "%10" IDD , 21) ;
2733 *indn = 8 ;
2734 }
2735 else if (xmin >= -999999999. && xmax <= 9999999999.)
2736 {
2737 strncpy (indfmt, "(7I11) ", 21) ;
2738 strncpy (indcfm, "%11" IDD , 21) ;
2739 *indn = 7 ;
2740 }
2741 else if (xmin >= -9999999999. && xmax <= 99999999999.)
2742 {
2743 strncpy (indfmt, "(6I12) ", 21) ;
2744 strncpy (indcfm, "%12" IDD , 21) ;
2745 *indn = 6 ;
2746 }
2747 else if (xmin >= -99999999999. && xmax <= 999999999999.)
2748 {
2749 strncpy (indfmt, "(6I13) ", 21) ;
2750 strncpy (indcfm, "%13" IDD , 21) ;
2751 *indn = 6 ;
2752 }
2753 else if (xmin >= -999999999999. && xmax <= 9999999999999.)
2754 {
2755 strncpy (indfmt, "(5I14) ", 21) ;
2756 strncpy (indcfm, "%14" IDD , 21) ;
2757 *indn = 5 ;
2758 }
2759 else
2760 {
2761 strncpy (indfmt, "(5I15) ", 21) ;
2762 strncpy (indcfm, "%15" IDD , 21) ;
2763 *indn = 5 ;
2764 }
2765 }
2766
2767
2768 /* -------------------------------------------------------------------------- */
2769 /* RBcards: determine number of cards required */
2770 /* -------------------------------------------------------------------------- */
2771
RB(cards)2772 PRIVATE Int RB(cards)
2773 (
2774 Int nitems, /* number of items to print */
2775 Int nperline /* number of items per line */
2776 )
2777 {
2778 return ((nitems == 0) ? 0 : (((nitems-1) / nperline) + 1)) ;
2779 }
2780
2781
2782 /* -------------------------------------------------------------------------- */
2783 /* RBfill: fill a string */
2784 /* -------------------------------------------------------------------------- */
2785
RB(fill)2786 PRIVATE void RB(fill)
2787 (
2788 char *s, /* string to fill */
2789 Int len, /* length of s (including trailing '\0') */
2790 char c /* character to fill s with */
2791 )
2792 {
2793 Int i ;
2794 for (i = 0 ; i < len ; i++)
2795 {
2796 s [i] = c ;
2797 }
2798 s [len-1] = '\0' ;
2799 }
2800
2801
2802 /* -------------------------------------------------------------------------- */
2803 /* RBok: verify a sparse matrix */
2804 /* -------------------------------------------------------------------------- */
2805
RB(ok)2806 PUBLIC Int RB(ok) /* 0:OK, < 0: error, > 0: warning */
2807 (
2808 /* inputs, not modified */
2809 Int nrow, /* number of rows */
2810 Int ncol, /* number of columns */
2811 Int nzmax, /* max # of entries */
2812 Int *Ap, /* size ncol+1, column pointers */
2813 Int *Ai, /* size nz = Ap [ncol], row indices */
2814 double *Ax, /* real part, or both if merged-complex */
2815 double *Az, /* imaginary part for split-complex */
2816 char *As, /* logical matrices (useful for MATLAB caller only) */
2817 Int mkind, /* 0:real, 1:logical/pattern, 2:split-complex, 3:integer,
2818 4:merged-complex */
2819
2820 /* outputs, not defined on input */
2821 Int *p_njumbled, /* # of jumbled row indices (-1 if not computed) */
2822 Int *p_nzeros /* number of explicit zeros (-1 if not computed) */
2823 )
2824 {
2825 double xr, xz ;
2826 Int i, j, p, pend, njumbled, nzeros, ilast ;
2827
2828 /* ---------------------------------------------------------------------- */
2829 /* in case of early return */
2830 /* ---------------------------------------------------------------------- */
2831
2832 if (p_njumbled) *p_njumbled = -1 ;
2833 if (p_nzeros ) *p_nzeros = -1 ;
2834
2835 if (mkind < 0 || mkind > 4)
2836 {
2837 return (RBIO_MKIND_INVALID) ;
2838 }
2839
2840 /* ---------------------------------------------------------------------- */
2841 /* check the dimensions */
2842 /* ---------------------------------------------------------------------- */
2843
2844 if (nrow < 0 || ncol < 0 || nzmax < 0)
2845 {
2846 return (RBIO_DIM_INVALID) ;
2847 }
2848
2849 /* ---------------------------------------------------------------------- */
2850 /* check the column pointers */
2851 /* ---------------------------------------------------------------------- */
2852
2853 if (Ap == NULL || Ap [0] != 0)
2854 {
2855 /* column pointers invalid */
2856 return (RBIO_CP_INVALID) ;
2857 }
2858 for (j = 0 ; j < ncol ; j++)
2859 {
2860 p = Ap [j] ;
2861 pend = Ap [j+1] ;
2862 if (pend < p || pend > nzmax)
2863 {
2864 /* column pointers not monotonically non-decreasing */
2865 return (RBIO_CP_INVALID) ;
2866 }
2867 }
2868
2869 /* ---------------------------------------------------------------------- */
2870 /* check the row indices and numerical values */
2871 /* ---------------------------------------------------------------------- */
2872
2873 if (Ai == NULL)
2874 {
2875 /* row indices not present */
2876 return (RBIO_ROW_INVALID) ;
2877 }
2878
2879 njumbled = 0 ;
2880 nzeros = 0 ;
2881
2882 for (j = 0 ; j < ncol ; j++)
2883 {
2884 ilast = -1 ;
2885 for (p = Ap [j] ; p < Ap [j+1] ; p++)
2886 {
2887 i = Ai [p] ;
2888 if (i < 0 || i >= nrow)
2889 {
2890 /* row indices out of range */
2891 return (RBIO_ROW_INVALID) ;
2892 }
2893 if (i <= ilast)
2894 {
2895 /* row indices unsorted, or duplicates present */
2896 njumbled++ ;
2897 }
2898 if (mkind == 1 && As)
2899 {
2900 xr = (double) (As [p]) ;
2901 xz = 0 ;
2902 }
2903 else
2904 {
2905 RB(get_entry) (mkind, Ax, Az, p, &xr, &xz) ;
2906 }
2907 if (xr == 0 && xz == 0)
2908 {
2909 /* an explicit zero is present */
2910 nzeros++ ;
2911 }
2912 ilast = i ;
2913 }
2914 }
2915
2916 /* ---------------------------------------------------------------------- */
2917 /* return results */
2918 /* ---------------------------------------------------------------------- */
2919
2920 if (p_njumbled) *p_njumbled = njumbled ;
2921 if (p_nzeros ) *p_nzeros = nzeros ;
2922 return ((njumbled > 0) ? RBIO_JUMBLED : RBIO_OK) ;
2923 }
2924