1 /*
2 * Name: lpio.c
3 * Author: Pietro Belotti
4 * Purpose: read compact .bz file into sparse matrix data structures
5 *
6 * This code is published under the Eclipse Public License (EPL).
7 * See http://www.eclipse.org/legal/epl-v10.html
8 *
9 */
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include <math.h>
15 #include <bzlib.h>
16
17 #include "sparse.h"
18 #include "lpio.h"
19 #include "rtr.h"
20 #include "misc.h"
21
22 #ifdef RTR_MPI
23 #include <mpi.h>
24
25 enum {RTR_INIT_SIZE, RTR_INIT_INT, RTR_INIT_DOUBLE};
26 #endif
27
28 /*
29 * scan sequence of bzip2'ed files (each contains a subset of rows),
30 * partition into chunks and send them
31 */
32
read_problem(sparseLP * lp,char ** filename)33 int read_problem (sparseLP *lp, /* LP data */
34 char **filename /* file name */
35 ) {
36
37 char first = 1;
38
39 FILE *f;
40
41 int status;
42
43 BZFILE *bzf;
44
45 int rownnz, curnnz,
46 totrow, currow = 0,
47 nnz, true_nnz,
48 i, j, best,
49 pos_i = 0, pos_d = 0,
50 size_d, size_i,
51 size_rl, size_ru,
52 target_nnz,
53 n_left_nnz;
54
55 int bcast [4];
56 int pos;
57
58 static int *bufi;
59 static double *bufd;
60
61 double lhsmax, lhsmin;
62
63 double *rlb;
64 double *rub;
65 double *rhs;
66
67 double coe;
68
69 double *lb, *ub;
70
71 if (lp == NULL) {
72
73 /*
74 This function is called at the end with lp == NULL. The sole
75 purpose of this call is to provide a mechanism to free bufi and
76 bufd. That's why those variables are declared static, too.
77 */
78
79 if (bufi) free (bufi);
80 if (bufd) free (bufd);
81
82 return 0;
83 }
84
85 /*
86 * Read chunks of matrix from .bz2 file
87 */
88
89 if (lp -> my_id == 0) {
90
91 lp->ntaut = lp->niis = 0;
92
93 /*
94 * init bz2'ed archive
95 */
96
97 if (strcmp (*filename, "-")) f = fopen (*filename, "r");
98 else f = stdin;
99
100 if ((f == NULL) || ferror (f)) {
101 printf ("unable to open file %s\n", *filename);
102 return -1;
103 }
104
105 printf ("%6.2f Opening bz2 file %s\n", CoinCpuTime(), *filename); fflush (stdout);
106
107 bzf = BZ2_bzReadOpen (&status, f, 1, 0, NULL, 0);
108
109 /*
110 * read scalar parameters (file header)
111 */
112
113 bzgetint (bzf, &(lp -> c0));
114 bzgetint (bzf, &(lp -> r0));
115 bzgetint (bzf, &nnz);
116 bzgetint (bzf, &best);
117
118 n_left_nnz = true_nnz = nnz;
119
120 /*
121 * Send this initial information to slaves so that they prepare
122 * the tables of the right size
123 */
124
125 bcast [0] = nnz;
126 bcast [1] = lp -> c0;
127 bcast [2] = lp -> r0;
128 bcast [3] = lp -> ncpus;
129
130 #ifdef RTR_MPI
131 MPI_Bcast (bcast, 4, MPI_INT, 0, MPI_COMM_WORLD); // 0 is the master node
132 #endif
133
134 /*
135 * create double and int buffers for the data to be sent to slaves
136 */
137
138 /* target_nnz = nnz / lp -> ncpus + 1; */
139
140 size_i = size_d = 0;
141
142 bufd = NULL;
143 bufi = NULL;
144
145 /*
146 * store variable upper/lower bounds (not alternate, but in
147 * separate areas of the array)
148 */
149
150 for (i=lp->c0; i>0; i--, pos_d++) {
151
152 reallocate_double (pos_d + lp -> c0, &size_d, &bufd);
153
154 bzgetdbl (bzf, bufd + pos_d);
155 bzgetdbl (bzf, bufd + pos_d + lp -> c0);
156 }
157
158 lb = bufd;
159 ub = bufd + lp->c0;
160
161 /*
162 * allocate initial space for rhs (can be expanded with realloc)
163 */
164
165 size_rl = size_ru = 0;
166
167 rlb = rub = NULL;
168
169 /*
170 * read row parameters and send chunk to each slave
171 */
172
173 /* printf ("%6.2f: spreading matrix chunks...\n", CoinCpuTime()); */
174
175 target_nnz = nnz / lp -> ncpus + 1;
176
177 for (totrow = i = 0; i < lp -> ncpus; i++) {
178
179 if (!(i % CHUNKS_PER_LINE)) {
180 printf ("%6.2f chunks %d..%d\n", CoinCpuTime (), i, mymin (lp->ncpus-1, i + CHUNKS_PER_LINE - 1));
181 /* fflush (stdout); */
182 }
183
184 pos_d = 2 * lp->c0;
185 pos_i = 0;
186
187 /*
188 * read each row and enqueue coefficients and indices
189 */
190
191 for (curnnz = currow = 0;
192
193 ((i == lp->ncpus-1) ||
194 (curnnz < target_nnz)) &&
195 (totrow < lp -> r0);
196
197 currow++, totrow++) {
198
199 /*
200 * read this row's nnz
201 */
202
203 bzgetint (bzf, &rownnz);
204
205 /* printf ("%d: %d\n", currow, rownnz); fflush (stdout); */
206
207 /* printf ("line %d has %d nnzs (to be drawn from %d) %d < %d/%d\n",currow, rownnz, n_left_nnz); */
208
209 n_left_nnz -= rownnz;
210
211 curnnz += rownnz;
212
213 reallocate_int (pos_i, &size_i, &bufi);
214
215 bufi [pos_i++] = rownnz;
216
217 reallocate_double (currow, &size_rl, &rlb);
218 reallocate_double (currow, &size_ru, &rub);
219
220 /*
221 * read b and c in b <= Ax <= c
222 */
223
224 bzgetdbl (bzf, rlb + currow);
225 bzgetdbl (bzf, rub + currow);
226
227 if ((rlb [currow] > -1e20) && (rub [currow] < 1e20) && first) {
228
229 first = 0;
230 printf ("warning: range constraints (first found at %d: [%f,%f]) not yet implemented\n",
231 currow, rlb [currow], rub [currow]);
232 }
233
234 /*
235 * read coefficients and indices
236 */
237
238 lhsmin = lhsmax = 0;
239
240 #ifdef RTR_USE_PRAGMAS
241 #pragma execution_frequency (very_high)
242 #endif
243
244 for (j=rownnz; j>0; --j, pos_i++, pos_d++) {
245
246 bzgetint (bzf, &pos);
247 if ((bzgetdbl (bzf, &coe) == 1) && (*++filename)) {
248
249 BZ2_bzReadClose (&status, bzf);
250
251 fclose (f);
252
253 if (strcmp (*filename, "-")) f = fopen (*filename, "r");
254 else f = stdin;
255
256 if ((f == NULL) || ferror (f)) {
257 printf ("unable to open file %s\n", *filename);
258 return -1;
259 }
260
261 printf ("%6.2f Opening bz2 file %s\n", CoinCpuTime(), *filename); fflush (stdout);
262
263 bzf = BZ2_bzReadOpen( &status, f, 1, 0, NULL, 0);
264
265 /* printf ("acabou!\n"); */
266 }
267
268 /* #pragma disjoint (*bufd, *bufi, coe, pos, pos_d, pos_i, lhsmax, lhsmin) */
269
270 {
271 reallocate_double (pos_d, &size_d, &bufd);
272 reallocate_int (pos_i, &size_i, &bufi);
273
274 lb = bufd;
275 ub = bufd + lp->c0;
276
277 bufd [pos_d] = coe;
278 bufi [pos_i] = pos;
279
280 if (!(lp->noprep)) {
281
282 if (coe > 1e-20) {
283 lhsmax += coe * ub [pos];
284 lhsmin += coe * lb [pos];
285 }
286 else if (coe < -1e-20) {
287 lhsmax += coe * lb [pos];
288 lhsmin += coe * ub [pos];
289 }
290 }
291 }
292 }
293
294 if (!(lp->noprep)) {
295
296 if (((rlb [currow] > -1e20) && (lhsmax < rlb [currow])) ||
297 ((rub [currow] < 1e20) && (lhsmin > rub [currow]))) { /* iis of size 1 */
298
299 /* printf ("iis %d: %.3e < %.3e or %.3e > %.3e\n", currow, lhsmax, rlb [currow], lhsmin, rub [currow]); */
300
301 pos_i -= (rownnz + 1);
302 pos_d -= rownnz;
303 curnnz -= rownnz;
304 --currow;
305
306 ++ (lp -> niis);
307
308 true_nnz -= rownnz;
309 target_nnz = (true_nnz) / lp -> ncpus;
310 }
311 else if (((rub [currow] > 1e20) || (lhsmax <= rub [currow])) &&
312 ((rlb [currow] < -1e20) || (lhsmin >= rlb [currow]))) { /* tautology */
313
314 pos_i -= (rownnz + 1);
315 pos_d -= rownnz;
316 curnnz -= rownnz;
317 --currow;
318
319 ++ (lp -> ntaut);
320 true_nnz -= rownnz;
321 target_nnz = (true_nnz) / lp -> ncpus;
322 }
323 }
324 }
325
326 /*
327 * Add constraints rhs's to the end of bufd
328 */
329
330 /* #pragma disjoint (*bufd, *rlb, *rub) */
331
332 {
333 reallocate_double (pos_d + 2 * currow, &size_d, &bufd);
334
335 for (j=currow; j>0; j--) bufd [pos_d++] = *rlb++;
336 for (j=currow; j>0; j--) bufd [pos_d++] = *rub++;
337
338 rlb -= currow;
339 rub -= currow;
340 }
341
342 lb = bufd;
343 ub = bufd + lp->c0;
344
345 bcast [0] = pos_d;
346 bcast [1] = pos_i;
347 bcast [2] = currow;
348
349 #ifdef RTR_MPI
350 if (i < lp -> ncpus - 1) {
351 MPI_Send (bcast, 3, MPI_INT, i+1, RTR_INIT_SIZE, MPI_COMM_WORLD);
352 MPI_Send (bufd, pos_d, MPI_DOUBLE, i+1, RTR_INIT_DOUBLE, MPI_COMM_WORLD);
353 MPI_Send (bufi, pos_i, MPI_INT, i+1, RTR_INIT_INT, MPI_COMM_WORLD);
354 }
355 #endif
356
357 bufd = (double *) realloc (bufd, (size_d = pos_d) * sizeof (double));
358 bufi = (int *) realloc (bufi, (size_i = pos_i) * sizeof (int));
359 rlb = (double *) realloc (rlb, (size_rl = currow) * sizeof (double));
360 rub = (double *) realloc (rub, (size_ru = currow) * sizeof (double));
361
362 lb = bufd;
363 ub = bufd + lp->c0;
364
365 /*
366 if (currow < 1e4) printf ("[%d,", currow);
367 else if (currow < 1e7) printf ("[%dk,", currow/(int)1e3);
368 else if (currow < 1e10) printf ("[%dm,", currow/(int)1e6);
369 else printf ("[%dg,", currow/(int)1e9);
370
371 if (pos_i <= 1e4) printf ("%d] ", pos_i - currow);
372 else if (pos_i <= 1e7) printf ("%dk] ", (pos_i - currow)/(int)1e3);
373 else if (pos_i <= 1e7) printf ("%dm] ", (pos_i - currow)/(int)1e6);
374 else printf ("%dg] ", (pos_i - currow)/(int)1e9); fflush (stdout);
375
376 if (!((i+1) % CHUNKS_PER_LINE)) printf ("\n");
377 */
378
379 if (i < lp->ncpus - 1)
380 target_nnz = n_left_nnz / (lp->ncpus - i - 1);
381 }
382
383 printf ("%6.2f done. %d rows, %d columns, %d nonzero", CoinCpuTime (), lp->r0, lp->c0, nnz);
384
385 if (!(lp -> noprep) && (lp->niis || lp->ntaut))
386 printf ("\n %d iis, %d tautologies, %d nonzero eliminated", lp->niis, lp->ntaut, nnz - true_nnz);
387
388 printf ("\n");
389
390 free (rlb);
391 free (rub);
392
393 BZ2_bzReadClose (&status, bzf);
394
395 fclose (f);
396
397 } else {
398
399 #ifdef RTR_MPI
400 MPI_Status status;
401
402 MPI_Bcast (bcast, 4, MPI_INT, 0, MPI_COMM_WORLD);
403
404 nnz = bcast [0];
405 lp -> c0 = bcast [1];
406 lp -> r0 = bcast [2];
407 lp -> ncpus = bcast [3];
408
409 //
410 // Stay idle until my own data comes...
411 //
412
413 MPI_Recv (bcast, 3, MPI_INT, 0, RTR_INIT_SIZE, MPI_COMM_WORLD, &status);
414
415 pos_d = size_d = bcast [0];
416 pos_i = size_i = bcast [1];
417 currow = bcast [2];
418
419 bufi = (int *) malloc (size_i * sizeof (int));
420 bufd = (double *) malloc (size_d * sizeof (double));
421
422 MPI_Recv (bufd, size_d, MPI_DOUBLE, 0, RTR_INIT_DOUBLE, MPI_COMM_WORLD, &status);
423 MPI_Recv (bufi, size_i, MPI_INT, 0, RTR_INIT_INT, MPI_COMM_WORLD, &status);
424 #endif
425 }
426
427 /*
428 * create sparse representation of local submatrix A_j
429 */
430
431 lp -> rk = currow;
432 lp -> nnzk = pos_d - 2 * (lp -> rk + lp -> c0);
433
434 lp -> lb = bufd;
435 lp -> ub = lp -> lb + lp -> c0;
436
437 lp -> rlb = lp -> ub + lp -> c0 + lp -> nnzk;
438 lp -> rub = lp -> rlb + lp -> rk;
439
440 rhs = lp -> rhs = (double *) malloc (currow * sizeof (double));
441
442 rlb = lp -> rlb;
443 rub = lp -> rub;
444
445 for (j=currow; j>0; j--, rub++, rlb++, rhs++)
446
447 if (*rlb < -1e20) *rhs = *rub;
448 else *rhs = *rlb;
449
450 rhs -= currow;
451 rlb -= currow;
452 rub -= currow;
453
454 lp -> ic = (double **) malloc (currow * sizeof (double *));
455 lp -> ip = (int **) malloc (currow * sizeof (int *));
456 lp -> il = (int *) malloc (currow * sizeof (int));
457
458 lp -> chosen = (char *) malloc (currow * sizeof (char));
459
460 pos_i = 0;
461 pos_d = 2 * lp -> c0;
462
463 {
464 double **ic = lp -> ic;
465 int **ip = lp -> ip;
466 int *il = lp -> il;
467 char *ch = lp -> chosen;
468
469 for (j=currow; j>0; j--) {
470
471 *ch++ = 0;
472 *ic++ = bufd + pos_d;
473 *ip++ = bufi + pos_i + 1;
474 *il++ = bufi [pos_i];
475
476 pos_d += bufi [pos_i];
477 pos_i += bufi [pos_i] + 1;
478 }
479
480 ic -= currow;
481 il -= currow;
482
483 /*
484 * Normalization: the constraint must be of the form ax >= b with
485 * ||a||=1
486 */
487
488 for (i=lp->rk; i>0; i--, ic++, il++, rlb++, rub++, rhs++) {
489
490 register double norm = get_norm (*ic, *il);
491
492 if (norm > EPSILON) {
493
494 register double *icl = *ic;
495
496 if (*rub < 1e29) {
497
498 norm = - norm;
499 *rhs = *rub / norm;
500 }
501 else *rhs = *rlb / norm;
502
503 for (j = *il; j--;) *icl++ /= norm;
504 }
505 }
506 }
507
508 /* if (lp->my_id==0) printf ("%6.2f create transposed matrix\n", CoinCpuTime ()); */
509
510 create_transpose (lp);
511
512 return status;
513 }
514