1 /**
2 @file
3
4 @ingroup cudd
5
6 @brief Function to read a matrix in Harwell format.
7
8 @author Fabio Somenzi
9
10 @copyright@parblock
11 Copyright (c) 1995-2015, Regents of the University of Colorado
12
13 All rights reserved.
14
15 Redistribution and use in source and binary forms, with or without
16 modification, are permitted provided that the following conditions
17 are met:
18
19 Redistributions of source code must retain the above copyright
20 notice, this list of conditions and the following disclaimer.
21
22 Redistributions in binary form must reproduce the above copyright
23 notice, this list of conditions and the following disclaimer in the
24 documentation and/or other materials provided with the distribution.
25
26 Neither the name of the University of Colorado nor the names of its
27 contributors may be used to endorse or promote products derived from
28 this software without specific prior written permission.
29
30 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
31 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
32 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
33 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
34 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
35 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
36 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
37 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
38 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
39 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
40 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
41 POSSIBILITY OF SUCH DAMAGE.
42 @endparblock
43
44 */
45
46 #include "util.h"
47 #include "cuddInt.h"
48
49 /*---------------------------------------------------------------------------*/
50 /* Constant declarations */
51 /*---------------------------------------------------------------------------*/
52
53
54 /*---------------------------------------------------------------------------*/
55 /* Stucture declarations */
56 /*---------------------------------------------------------------------------*/
57
58
59 /*---------------------------------------------------------------------------*/
60 /* Type declarations */
61 /*---------------------------------------------------------------------------*/
62
63
64 /*---------------------------------------------------------------------------*/
65 /* Variable declarations */
66 /*---------------------------------------------------------------------------*/
67
68
69 /*---------------------------------------------------------------------------*/
70 /* Macro declarations */
71 /*---------------------------------------------------------------------------*/
72
73 /** \cond */
74
75 /*---------------------------------------------------------------------------*/
76 /* Static function prototypes */
77 /*---------------------------------------------------------------------------*/
78
79 /** \endcond */
80
81
82 /*---------------------------------------------------------------------------*/
83 /* Definition of exported functions */
84 /*---------------------------------------------------------------------------*/
85
86 /**
87 @brief Reads in a matrix in the format of the Harwell-Boeing
88 benchmark suite.
89
90 @details The variables are ordered as follows:
91 <blockquote>
92 x\[0\] y\[0\] x\[1\] y\[1\] ...
93 </blockquote>
94 0 is the most significant bit. On input, nx and ny hold the numbers
95 of row and column variables already in existence.
96
97 @return 1 on success; 0 otherwise.
98
99 @sideeffect On output, nx and ny hold the numbers of row and column
100 variables actually used by the matrix. m and n are set to the
101 numbers of rows and columns of the matrix. Their values on input
102 are immaterial. The %ADD for the sparse matrix is returned in E, and
103 its reference count is > 0.
104
105 @see Cudd_addRead Cudd_bddRead
106
107 */
108 int
Cudd_addHarwell(FILE * fp,DdManager * dd,DdNode ** E,DdNode *** x,DdNode *** y,DdNode *** xn,DdNode *** yn_,int * nx,int * ny,int * m,int * n,int bx,int sx,int by,int sy,int pr)109 Cudd_addHarwell(
110 FILE * fp /**< pointer to the input file */,
111 DdManager * dd /**< %DD manager */,
112 DdNode ** E /**< characteristic function of the graph */,
113 DdNode *** x /**< array of row variables */,
114 DdNode *** y /**< array of column variables */,
115 DdNode *** xn /**< array of complemented row variables */,
116 DdNode *** yn_ /**< array of complemented column variables */,
117 int * nx /**< number or row variables */,
118 int * ny /**< number or column variables */,
119 int * m /**< number of rows */,
120 int * n /**< number of columns */,
121 int bx /**< first index of row variables */,
122 int sx /**< step of row variables */,
123 int by /**< first index of column variables */,
124 int sy /**< step of column variables */,
125 int pr /**< verbosity level */)
126 {
127 DdNode *one, *zero;
128 DdNode *w;
129 DdNode *cubex, *cubey, *minterm1;
130 int u, v, err, i, j, nv;
131 double val;
132 /* local copies of x, y, xn, yn_ */
133 DdNode **lx = NULL, **ly = NULL, **lxn = NULL, **lyn = NULL;
134 int lnx, lny; /* local copies of nx and ny */
135 char title[73], key[9], mxtype[4], rhstyp[4];
136 int totcrd, ptrcrd, indcrd, valcrd, rhscrd,
137 nrow, ncol, nnzero, neltvl,
138 nrhs, nrhsix;
139 int *colptr, *rowind;
140 #if 0
141 int nguess, nexact;
142 int *rhsptr, *rhsind;
143 #endif
144
145 if (*nx < 0 || *ny < 0) return(0);
146
147 one = DD_ONE(dd);
148 zero = DD_ZERO(dd);
149
150 /* Read the header */
151 err = fscanf(fp, "%72c %8c", title, key);
152 if (err == EOF) {
153 return(0);
154 } else if (err != 2) {
155 return(0);
156 }
157 title[72] = (char) 0;
158 key[8] = (char) 0;
159
160 err = fscanf(fp, "%d %d %d %d %d", &totcrd, &ptrcrd, &indcrd,
161 &valcrd, &rhscrd);
162 if (err == EOF) {
163 return(0);
164 } else if (err != 5) {
165 return(0);
166 }
167
168 err = fscanf(fp, "%3s %d %d %d %d", mxtype, &nrow, &ncol,
169 &nnzero, &neltvl);
170 if (err == EOF) {
171 return(0);
172 } else if (err != 5) {
173 return(0);
174 }
175
176 /* Skip FORTRAN formats */
177 if (rhscrd == 0) {
178 err = fscanf(fp, "%*s %*s %*s \n");
179 } else {
180 err = fscanf(fp, "%*s %*s %*s %*s \n");
181 }
182 if (err == EOF) {
183 return(0);
184 } else if (err != 0) {
185 return(0);
186 }
187
188 /* Print out some stuff if requested to be verbose */
189 if (pr>0) {
190 (void) fprintf(dd->out,"%s: type %s, %d rows, %d columns, %d entries\n", key,
191 mxtype, nrow, ncol, nnzero);
192 if (pr>1) (void) fprintf(dd->out,"%s\n", title);
193 }
194
195 /* Check matrix type */
196 if (mxtype[0] != 'R' || mxtype[1] != 'U' || mxtype[2] != 'A') {
197 (void) fprintf(dd->err,"%s: Illegal matrix type: %s\n",
198 key, mxtype);
199 return(0);
200 }
201 if (neltvl != 0) return(0);
202
203 /* Read optional 5-th line */
204 if (rhscrd != 0) {
205 err = fscanf(fp, "%3c %d %d", rhstyp, &nrhs, &nrhsix);
206 if (err == EOF) {
207 return(0);
208 } else if (err != 3) {
209 return(0);
210 }
211 rhstyp[3] = (char) 0;
212 if (rhstyp[0] != 'F') {
213 (void) fprintf(dd->err,
214 "%s: Sparse right-hand side not yet supported\n", key);
215 return(0);
216 }
217 if (pr>0) (void) fprintf(dd->out,"%d right-hand side(s)\n", nrhs);
218 } else {
219 nrhs = 0;
220 }
221
222 /* Compute the number of variables */
223
224 /* row and column numbers start from 0 */
225 u = nrow - 1;
226 for (i=0; u > 0; i++) {
227 u >>= 1;
228 }
229 lnx = i;
230 if (nrhs == 0) {
231 v = ncol - 1;
232 } else {
233 v = 2* (ddMax(ncol, nrhs) - 1);
234 }
235 for (i=0; v > 0; i++) {
236 v >>= 1;
237 }
238 lny = i;
239
240 /* Allocate or reallocate arrays for variables as needed */
241 if (*nx == 0) {
242 if (lnx > 0) {
243 *x = lx = ALLOC(DdNode *,lnx);
244 if (lx == NULL) {
245 dd->errorCode = CUDD_MEMORY_OUT;
246 return(0);
247 }
248 *xn = lxn = ALLOC(DdNode *,lnx);
249 if (lxn == NULL) {
250 dd->errorCode = CUDD_MEMORY_OUT;
251 return(0);
252 }
253 } else {
254 *x = *xn = NULL;
255 }
256 } else if (lnx > *nx) {
257 *x = lx = REALLOC(DdNode *, *x, lnx);
258 if (lx == NULL) {
259 dd->errorCode = CUDD_MEMORY_OUT;
260 return(0);
261 }
262 *xn = lxn = REALLOC(DdNode *, *xn, lnx);
263 if (lxn == NULL) {
264 dd->errorCode = CUDD_MEMORY_OUT;
265 return(0);
266 }
267 } else {
268 lx = *x;
269 lxn = *xn;
270 }
271 if (*ny == 0) {
272 if (lny >0) {
273 *y = ly = ALLOC(DdNode *,lny);
274 if (ly == NULL) {
275 dd->errorCode = CUDD_MEMORY_OUT;
276 return(0);
277 }
278 *yn_ = lyn = ALLOC(DdNode *,lny);
279 if (lyn == NULL) {
280 dd->errorCode = CUDD_MEMORY_OUT;
281 return(0);
282 }
283 } else {
284 *y = *yn_ = NULL;
285 }
286 } else if (lny > *ny) {
287 *y = ly = REALLOC(DdNode *, *y, lny);
288 if (ly == NULL) {
289 dd->errorCode = CUDD_MEMORY_OUT;
290 return(0);
291 }
292 *yn_ = lyn = REALLOC(DdNode *, *yn_, lny);
293 if (lyn == NULL) {
294 dd->errorCode = CUDD_MEMORY_OUT;
295 return(0);
296 }
297 } else {
298 ly = *y;
299 lyn = *yn_;
300 }
301
302 /* Create new variables as needed */
303 for (i= *nx,nv=bx+(*nx)*sx; i < lnx; i++,nv+=sx) {
304 do {
305 dd->reordered = 0;
306 lx[i] = cuddUniqueInter(dd, nv, one, zero);
307 } while (dd->reordered == 1);
308 if (lx[i] == NULL) {
309 if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) {
310 dd->timeoutHandler(dd, dd->tohArg);
311 }
312 return(0);
313 }
314 cuddRef(lx[i]);
315 do {
316 dd->reordered = 0;
317 lxn[i] = cuddUniqueInter(dd, nv, zero, one);
318 } while (dd->reordered == 1);
319 if (lxn[i] == NULL) {
320 if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) {
321 dd->timeoutHandler(dd, dd->tohArg);
322 }
323 return(0);
324 }
325 cuddRef(lxn[i]);
326 }
327 for (i= *ny,nv=by+(*ny)*sy; i < lny; i++,nv+=sy) {
328 do {
329 dd->reordered = 0;
330 ly[i] = cuddUniqueInter(dd, nv, one, zero);
331 } while (dd->reordered == 1);
332 if (ly[i] == NULL) {
333 if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) {
334 dd->timeoutHandler(dd, dd->tohArg);
335 }
336 return(0);
337 }
338 cuddRef(ly[i]);
339 do {
340 dd->reordered = 0;
341 lyn[i] = cuddUniqueInter(dd, nv, zero, one);
342 } while (dd->reordered == 1);
343 if (lyn[i] == NULL) {
344 if (dd->errorCode == CUDD_TIMEOUT_EXPIRED && dd->timeoutHandler) {
345 dd->timeoutHandler(dd, dd->tohArg);
346 }
347 return(0);
348 }
349 cuddRef(lyn[i]);
350 }
351
352 /* Update matrix parameters */
353 *nx = lnx;
354 *ny = lny;
355 *m = nrow;
356 if (nrhs == 0) {
357 *n = ncol;
358 } else {
359 *n = (1 << (lny - 1)) + nrhs;
360 }
361
362 /* Read structure data */
363 colptr = ALLOC(int, ncol+1);
364 if (colptr == NULL) {
365 dd->errorCode = CUDD_MEMORY_OUT;
366 return(0);
367 }
368 rowind = ALLOC(int, nnzero);
369 if (rowind == NULL) {
370 dd->errorCode = CUDD_MEMORY_OUT;
371 return(0);
372 }
373
374 for (i=0; i<ncol+1; i++) {
375 err = fscanf(fp, " %d ", &u);
376 if (err == EOF){
377 FREE(colptr);
378 FREE(rowind);
379 return(0);
380 } else if (err != 1) {
381 FREE(colptr);
382 FREE(rowind);
383 return(0);
384 }
385 colptr[i] = u - 1;
386 }
387 if (colptr[0] != 0) {
388 (void) fprintf(dd->err,"%s: Unexpected colptr[0] (%d)\n",
389 key,colptr[0]);
390 FREE(colptr);
391 FREE(rowind);
392 return(0);
393 }
394 for (i=0; i<nnzero; i++) {
395 err = fscanf(fp, " %d ", &u);
396 if (err == EOF){
397 FREE(colptr);
398 FREE(rowind);
399 return(0);
400 } else if (err != 1) {
401 FREE(colptr);
402 FREE(rowind);
403 return(0);
404 }
405 rowind[i] = u - 1;
406 }
407
408 *E = zero; cuddRef(*E);
409
410 for (j=0; j<ncol; j++) {
411 v = j;
412 cubey = one; cuddRef(cubey);
413 for (nv = lny - 1; nv>=0; nv--) {
414 if (v & 1) {
415 w = Cudd_addApply(dd, Cudd_addTimes, cubey, ly[nv]);
416 } else {
417 w = Cudd_addApply(dd, Cudd_addTimes, cubey, lyn[nv]);
418 }
419 if (w == NULL) {
420 Cudd_RecursiveDeref(dd, cubey);
421 FREE(colptr);
422 FREE(rowind);
423 return(0);
424 }
425 cuddRef(w);
426 Cudd_RecursiveDeref(dd, cubey);
427 cubey = w;
428 v >>= 1;
429 }
430 for (i=colptr[j]; i<colptr[j+1]; i++) {
431 u = rowind[i];
432 err = fscanf(fp, " %lf ", &val);
433 if (err == EOF || err != 1){
434 Cudd_RecursiveDeref(dd, cubey);
435 FREE(colptr);
436 FREE(rowind);
437 return(0);
438 }
439 /* Create new Constant node if necessary */
440 cubex = cuddUniqueConst(dd, (CUDD_VALUE_TYPE) val);
441 if (cubex == NULL) {
442 Cudd_RecursiveDeref(dd, cubey);
443 FREE(colptr);
444 FREE(rowind);
445 return(0);
446 }
447 cuddRef(cubex);
448
449 for (nv = lnx - 1; nv>=0; nv--) {
450 if (u & 1) {
451 w = Cudd_addApply(dd, Cudd_addTimes, cubex, lx[nv]);
452 } else {
453 w = Cudd_addApply(dd, Cudd_addTimes, cubex, lxn[nv]);
454 }
455 if (w == NULL) {
456 Cudd_RecursiveDeref(dd, cubey);
457 Cudd_RecursiveDeref(dd, cubex);
458 FREE(colptr);
459 FREE(rowind);
460 return(0);
461 }
462 cuddRef(w);
463 Cudd_RecursiveDeref(dd, cubex);
464 cubex = w;
465 u >>= 1;
466 }
467 minterm1 = Cudd_addApply(dd, Cudd_addTimes, cubey, cubex);
468 if (minterm1 == NULL) {
469 Cudd_RecursiveDeref(dd, cubey);
470 Cudd_RecursiveDeref(dd, cubex);
471 FREE(colptr);
472 FREE(rowind);
473 return(0);
474 }
475 cuddRef(minterm1);
476 Cudd_RecursiveDeref(dd, cubex);
477 w = Cudd_addApply(dd, Cudd_addPlus, *E, minterm1);
478 if (w == NULL) {
479 Cudd_RecursiveDeref(dd, cubey);
480 FREE(colptr);
481 FREE(rowind);
482 return(0);
483 }
484 cuddRef(w);
485 Cudd_RecursiveDeref(dd, minterm1);
486 Cudd_RecursiveDeref(dd, *E);
487 *E = w;
488 }
489 Cudd_RecursiveDeref(dd, cubey);
490 }
491 FREE(colptr);
492 FREE(rowind);
493
494 /* Read right-hand sides */
495 for (j=0; j<nrhs; j++) {
496 v = j + (1<< (lny-1));
497 cubey = one; cuddRef(cubey);
498 for (nv = lny - 1; nv>=0; nv--) {
499 if (v & 1) {
500 w = Cudd_addApply(dd, Cudd_addTimes, cubey, ly[nv]);
501 } else {
502 w = Cudd_addApply(dd, Cudd_addTimes, cubey, lyn[nv]);
503 }
504 if (w == NULL) {
505 Cudd_RecursiveDeref(dd, cubey);
506 return(0);
507 }
508 cuddRef(w);
509 Cudd_RecursiveDeref(dd, cubey);
510 cubey = w;
511 v >>= 1;
512 }
513 for (i=0; i<nrow; i++) {
514 u = i;
515 err = fscanf(fp, " %lf ", &val);
516 if (err == EOF || err != 1){
517 Cudd_RecursiveDeref(dd, cubey);
518 return(0);
519 }
520 /* Create new Constant node if necessary */
521 if (val == (double) 0.0) continue;
522 cubex = cuddUniqueConst(dd, (CUDD_VALUE_TYPE) val);
523 if (cubex == NULL) {
524 Cudd_RecursiveDeref(dd, cubey);
525 return(0);
526 }
527 cuddRef(cubex);
528
529 for (nv = lnx - 1; nv>=0; nv--) {
530 if (u & 1) {
531 w = Cudd_addApply(dd, Cudd_addTimes, cubex, lx[nv]);
532 } else {
533 w = Cudd_addApply(dd, Cudd_addTimes, cubex, lxn[nv]);
534 }
535 if (w == NULL) {
536 Cudd_RecursiveDeref(dd, cubey);
537 Cudd_RecursiveDeref(dd, cubex);
538 return(0);
539 }
540 cuddRef(w);
541 Cudd_RecursiveDeref(dd, cubex);
542 cubex = w;
543 u >>= 1;
544 }
545 minterm1 = Cudd_addApply(dd, Cudd_addTimes, cubey, cubex);
546 if (minterm1 == NULL) {
547 Cudd_RecursiveDeref(dd, cubey);
548 Cudd_RecursiveDeref(dd, cubex);
549 return(0);
550 }
551 cuddRef(minterm1);
552 Cudd_RecursiveDeref(dd, cubex);
553 w = Cudd_addApply(dd, Cudd_addPlus, *E, minterm1);
554 if (w == NULL) {
555 Cudd_RecursiveDeref(dd, cubey);
556 return(0);
557 }
558 cuddRef(w);
559 Cudd_RecursiveDeref(dd, minterm1);
560 Cudd_RecursiveDeref(dd, *E);
561 *E = w;
562 }
563 Cudd_RecursiveDeref(dd, cubey);
564 }
565
566 return(1);
567
568 } /* end of Cudd_addHarwell */
569
570
571 /*---------------------------------------------------------------------------*/
572 /* Definition of internal functions */
573 /*---------------------------------------------------------------------------*/
574
575 /*---------------------------------------------------------------------------*/
576 /* Definition of static functions */
577 /*---------------------------------------------------------------------------*/
578
579