1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "gretl_mpi.h"
22 #include "gretl_typemap.h"
23 #include <mpi.h>
24
25 #ifdef WIN32
26 # include <windows.h>
27 #else
28 # include <dlfcn.h>
29 #endif
30
31 /* Support for MPI in libgretl. On systems other than MS Windows
32 We get the MPI symbols that we need from the address space of
33 the calling program to avoid introducing a hard dependency on
34 libmpi; on Windows we call LoadLibrary() on msmpi.dll to the
35 same effect.
36
37 To use functions in this translation unit from elsewhere in
38 libgretl one must first call gretl_MPI_init(), and then
39 guard subsequent calls with "if (gretl_mpi_initialized())".
40 Otherwise you'll get a segfault.
41
42 For future reference, the MPI variants define the following
43 specific symbols in mpi.h:
44
45 Open MPI: OMPI_MAJOR_VERSION
46 MPICH: MPICH_VERSION
47 MS-MPI: MSMPI_VER
48 */
49
50 #ifdef OMPI_MAJOR_VERSION
51 /* external constants that need to be loaded from the
52 Open MPI library */
53 static MPI_Comm mpi_comm_world;
54 static MPI_Datatype mpi_double;
55 static MPI_Datatype mpi_int;
56 static MPI_Datatype mpi_unsigned;
57 static MPI_Datatype mpi_byte;
58 static MPI_Op mpi_sum;
59 static MPI_Op mpi_prod;
60 static MPI_Op mpi_max;
61 static MPI_Op mpi_min;
62 #else
63 /* It seems that MPICH and MS-MPI just define these symbols
64 as integer values in the header */
65 # define mpi_comm_world MPI_COMM_WORLD
66 # define mpi_double MPI_DOUBLE
67 # define mpi_int MPI_INT
68 # define mpi_unsigned MPI_UNSIGNED
69 # define mpi_byte MPI_BYTE
70 # define mpi_sum MPI_SUM
71 # define mpi_prod MPI_PROD
72 # define mpi_max MPI_MAX
73 # define mpi_min MPI_MIN
74 #endif
75
76 enum {
77 TAG_MATRIX_INFO = 1,
78 TAG_MATRIX_VAL,
79 TAG_SCALAR_VAL,
80 TAG_INT_VAL,
81 TAG_ARRAY_LEN,
82 TAG_SERIES_LEN,
83 TAG_SERIES_VAL,
84 TAG_LIST_LEN,
85 TAG_LIST_VAL,
86 TAG_STR_LEN,
87 TAG_STR_VAL,
88 TAG_BMEMB_INFO,
89 TAG_ARRAY_INFO,
90 TAG_BUNDLE_SIZE
91 };
92
93 #define MI_LEN 5 /* matrix info length */
94
95 #define MPI_DEBUG 0
96
97 static void *MPIhandle; /* handle to the MPI library */
98 static int gretl_MPI_err; /* initialization error record */
99 static int gretl_MPI_initted; /* are we initialized or not? */
100
101 /* renamed, pointerized versions of the MPI functions we need */
102 static int (*mpi_comm_rank) (MPI_Comm, int *);
103 static int (*mpi_comm_size) (MPI_Comm, int *);
104 static int (*mpi_error_class) (int, int *);
105 static int (*mpi_error_string) (int, char *, int *);
106 static int (*mpi_reduce) (void *, void *, int, MPI_Datatype, MPI_Op,
107 int, MPI_Comm);
108 static int (*mpi_allreduce) (void *, void *, int, MPI_Datatype, MPI_Op,
109 MPI_Comm);
110 static int (*mpi_bcast) (void *, int, MPI_Datatype, int, MPI_Comm);
111 static int (*mpi_send) (void *, int, MPI_Datatype, int, int, MPI_Comm);
112 static int (*mpi_recv) (void *, int, MPI_Datatype, int, int, MPI_Comm,
113 MPI_Status *);
114 static int (*mpi_barrier) (MPI_Comm);
115 static int (*mpi_probe) (int, int, MPI_Comm, MPI_Status *);
116 static double (*mpi_wtime) (void);
117 static int (*mpi_initialized) (int *);
118
119 static int gretl_matrix_bcast (gretl_matrix **pm, int id, int root);
120 static int gretl_array_bcast (gretl_array **pa, int id, int root);
121 static int gretl_bundle_send (gretl_bundle *b, int dest);
122 static gretl_bundle *gretl_bundle_receive (int source, int *err);
123
124 static void *mpi_receive_element (int source, GretlType etype,
125 int *err);
126
mpiget(void * handle,const char * name,int * err)127 static void *mpiget (void *handle, const char *name, int *err)
128 {
129 #ifdef WIN32
130 void *p = GetProcAddress(handle, name);
131 #else
132 void *p = dlsym(handle, name);
133 #endif
134
135 if (p == NULL) {
136 printf("mpiget: couldn't find '%s'\n", name);
137 *err += 1;
138 }
139
140 #if MPI_DEBUG
141 else {
142 printf("mpiget: '%s' -> %p\n", name, p);
143 }
144 #endif
145
146 return p;
147 }
148
gretl_MPI_init(void)149 int gretl_MPI_init (void)
150 {
151 int err = 0;
152
153 if (gretl_MPI_initted) {
154 return gretl_MPI_err;
155 }
156
157 #if MPI_DEBUG
158 printf("Loading MPI symbols...\n");
159 #endif
160
161 #ifdef WIN32
162 MPIhandle = LoadLibrary("msmpi.dll");
163 #else
164 MPIhandle = dlopen(NULL, RTLD_NOW);
165 #endif
166
167 if (MPIhandle == NULL) {
168 err = E_EXTERNAL;
169 goto bailout;
170 }
171
172 mpi_comm_rank = mpiget(MPIhandle, "MPI_Comm_rank", &err);
173 mpi_comm_size = mpiget(MPIhandle, "MPI_Comm_size", &err);
174 mpi_error_class = mpiget(MPIhandle, "MPI_Error_class", &err);
175 mpi_error_string = mpiget(MPIhandle, "MPI_Error_string", &err);
176 mpi_reduce = mpiget(MPIhandle, "MPI_Reduce", &err);
177 mpi_allreduce = mpiget(MPIhandle, "MPI_Allreduce", &err);
178 mpi_bcast = mpiget(MPIhandle, "MPI_Bcast", &err);
179 mpi_send = mpiget(MPIhandle, "MPI_Send", &err);
180 mpi_recv = mpiget(MPIhandle, "MPI_Recv", &err);
181 mpi_probe = mpiget(MPIhandle, "MPI_Probe", &err);
182 mpi_barrier = mpiget(MPIhandle, "MPI_Barrier", &err);
183 mpi_wtime = mpiget(MPIhandle, "MPI_Wtime", &err);
184 mpi_initialized = mpiget(MPIhandle, "MPI_Initialized", &err);
185
186 #ifdef OMPI_MAJOR_VERSION
187 if (!err) {
188 mpi_comm_world = (MPI_Comm) mpiget(MPIhandle, "ompi_mpi_comm_world", &err);
189 mpi_double = (MPI_Datatype) mpiget(MPIhandle, "ompi_mpi_double", &err);
190 mpi_int = (MPI_Datatype) mpiget(MPIhandle, "ompi_mpi_int", &err);
191 mpi_unsigned = (MPI_Datatype) mpiget(MPIhandle, "ompi_mpi_unsigned", &err);
192 mpi_byte = (MPI_Datatype) mpiget(MPIhandle, "ompi_mpi_byte", &err);
193 mpi_sum = (MPI_Op) mpiget(MPIhandle, "ompi_mpi_op_sum", &err);
194 mpi_prod = (MPI_Op) mpiget(MPIhandle, "ompi_mpi_op_prod", &err);
195 mpi_max = (MPI_Op) mpiget(MPIhandle, "ompi_mpi_op_max", &err);
196 mpi_min = (MPI_Op) mpiget(MPIhandle, "ompi_mpi_op_min", &err);
197 }
198 #endif
199
200 if (err) {
201 close_plugin(MPIhandle);
202 MPIhandle = NULL;
203 err = E_EXTERNAL;
204 }
205
206 bailout:
207
208 #if MPI_DEBUG
209 printf("load_MPI_symbols: returning %d\n", err);
210 #endif
211
212 gretl_MPI_err = err;
213 gretl_MPI_initted = 1;
214
215 return err;
216 }
217
gretl_mpi_error(int * err)218 static void gretl_mpi_error (int *err)
219 {
220 char msg[BUFSIZ];
221 int len, id = 0;
222
223 mpi_comm_rank(mpi_comm_world, &id);
224 mpi_error_string(*err, msg, &len);
225 gretl_errmsg_sprintf("%3d: %s", id, msg);
226
227 *err = E_EXTERNAL;
228 }
229
dim_error(int * dims,int n)230 static int dim_error (int *dims, int n)
231 {
232 int i, d0 = 0;
233 int err = 0;
234
235 /* allow for the possibility that some nodes
236 carry null matrices: d0 will be the common
237 non-zero dimension */
238
239 for (i=0; i<n; i++) {
240 if (dims[i] > 0) {
241 d0 = dims[i];
242 break;
243 }
244 }
245
246 /* check that all non-null matrices are conformable */
247
248 for (i=0; i<n; i++) {
249 if (dims[i] != d0 && dims[i] != 0) {
250 err = E_NONCONF;
251 break;
252 }
253 }
254
255 return err;
256 }
257
set_dim_max(int * dims,int * dtest,int n)258 static void set_dim_max (int *dims, int *dtest, int n)
259 {
260 int i;
261
262 /* find the max row dimension of matrices with
263 cols > 0, or vice versa */
264
265 for (i=0; i<n; i++) {
266 if (dtest[i] > 0) {
267 if (dims[i] > dims[n]) {
268 dims[n] = dims[i];
269 }
270 }
271 }
272 }
273
matrix_dims_check(int * rows,int * cols,int n,Gretl_MPI_Op op)274 static int matrix_dims_check (int *rows, int *cols, int n,
275 Gretl_MPI_Op op)
276 {
277 int match_rows =
278 (op == GRETL_MPI_SUM || op == GRETL_MPI_PROD ||
279 op == GRETL_MPI_HCAT);
280 int match_cols =
281 (op == GRETL_MPI_SUM || op == GRETL_MPI_PROD ||
282 op == GRETL_MPI_VCAT);
283 int err = 0;
284
285 if (match_rows) {
286 err = dim_error(rows, n);
287 }
288
289 if (!err && match_cols) {
290 err = dim_error(cols, n);
291 }
292
293 if (!err) {
294 set_dim_max(rows, cols, n);
295 set_dim_max(cols, rows, n);
296 /* for now we'll consider a null matrix reduction
297 as an error (though _maybe_ we want to permit
298 this?)
299 */
300 if (rows[n] == 0 || cols[n] == 0) {
301 err = E_DATA;
302 }
303 }
304
305 return err;
306 }
307
fill_matrix_info(int * rc,const gretl_matrix * m)308 static void fill_matrix_info (int *rc, const gretl_matrix *m)
309 {
310 rc[0] = m->rows;
311 rc[1] = m->cols;
312 rc[2] = m->is_complex;
313 rc[3] = gretl_matrix_get_t1(m);
314 rc[4] = gretl_matrix_get_t2(m);
315 }
316
matrix_reduce_alloc(int * rows,int * cols,int n,Gretl_MPI_Op op,gretl_matrix ** pm,double ** px)317 static int matrix_reduce_alloc (int *rows, int *cols,
318 int n, Gretl_MPI_Op op,
319 gretl_matrix **pm,
320 double **px)
321 {
322 int rtotal = 0, ctotal = 0;
323 int i, err = 0;
324
325 /* Note: in the arrays @rows and @cols the elements 0 to
326 n-1 are the values for each node, and the elements n
327 are the maxima.
328
329 The tasks here are (a) to allocate a matrix of the right
330 size to hold the result of the "reduce" operation and (b)
331 to allocate workspace big enough to hold the values of
332 the largest individual matrix.
333 */
334
335 if (op == GRETL_MPI_SUM || op == GRETL_MPI_PROD) {
336 rtotal = rows[n];
337 ctotal = cols[n];
338 } else if (op == GRETL_MPI_HCAT) {
339 /* there's a known common number of rows, but
340 we need to determine the total number of
341 columns in the result
342 */
343 rtotal = rows[n];
344 for (i=0; i<n; i++) {
345 if (rows[i] > 0) {
346 ctotal += cols[i];
347 }
348 }
349 } else if (op == GRETL_MPI_VCAT) {
350 /* there's a known common number of columns, but
351 we need to determine the total number of rows
352 in the result
353 */
354 ctotal = cols[n];
355 for (i=0; i<n; i++) {
356 if (cols[i] > 0) {
357 rtotal += rows[i];
358 }
359 }
360 }
361
362 if (rtotal == 0 || ctotal == 0) {
363 err = E_DATA;
364 } else if (op == GRETL_MPI_PROD) {
365 *pm = gretl_unit_matrix_new(rtotal, ctotal);
366 } else {
367 *pm = gretl_zero_matrix_new(rtotal, ctotal);
368 }
369
370 if (*pm == NULL) {
371 err = E_ALLOC;
372 } else {
373 size_t xsize = rows[n] * cols[n];
374
375 *px = malloc(xsize * sizeof **px);
376 if (*px == NULL) {
377 gretl_matrix_free(*pm);
378 *pm = NULL;
379 err = E_ALLOC;
380 }
381 }
382
383 return err;
384 }
385
matrix_reduce_step(gretl_matrix * mtarg,double * restrict src,int n,Gretl_MPI_Op op,int * offset)386 static int matrix_reduce_step (gretl_matrix *mtarg,
387 double * restrict src,
388 int n, Gretl_MPI_Op op,
389 int *offset)
390 {
391 double * restrict targ = mtarg->val;
392 int i;
393
394 if (op == GRETL_MPI_SUM) {
395 for (i=0; i<n; i++) {
396 targ[i] += src[i];
397 }
398 } else if (op == GRETL_MPI_PROD) {
399 for (i=0; i<n; i++) {
400 targ[i] *= src[i];
401 }
402 } else if (op == GRETL_MPI_HCAT) {
403 int k = *offset;
404
405 for (i=0; i<n; i++) {
406 targ[k++] = src[i];
407 }
408 *offset = k;
409 } else if (op == GRETL_MPI_VCAT) {
410 int rmin = *offset;
411 int nrows = n / mtarg->cols;
412 int rmax = rmin + nrows;
413 int j, k = 0;
414
415 for (j=0; j<mtarg->cols; j++) {
416 for (i=rmin; i<rmax; i++) {
417 gretl_matrix_set(mtarg, i, j, src[k++]);
418 }
419 }
420 *offset += nrows;
421 }
422
423 return 0;
424 }
425
invalid_rank_error(int r)426 static int invalid_rank_error (int r)
427 {
428 gretl_errmsg_sprintf(_("Invalid MPI rank %d"), r);
429 return E_DATA;
430 }
431
gretl_comm_check(int root,int * idp,int * npp)432 static int gretl_comm_check (int root, int *idp, int *npp)
433 {
434 int np;
435
436 mpi_comm_size(mpi_comm_world, &np);
437 if (root < 0 || root >= np) {
438 return invalid_rank_error(root);
439 }
440
441 if (npp != NULL) {
442 *npp = np;
443 }
444 if (idp != NULL) {
445 mpi_comm_rank(mpi_comm_world, idp);
446 }
447
448 return 0;
449 }
450
gretl_matrix_mpi_reduce(gretl_matrix * sm,gretl_matrix ** pm,Gretl_MPI_Op op,int root,gretlopt opt)451 int gretl_matrix_mpi_reduce (gretl_matrix *sm,
452 gretl_matrix **pm,
453 Gretl_MPI_Op op,
454 int root,
455 gretlopt opt)
456 {
457 gretl_matrix *rm = NULL;
458 double *val = NULL;
459 int *rows = NULL;
460 int *cols = NULL;
461 int rc[MI_LEN] = {0};
462 int id, np;
463 int err = 0;
464
465 if (op != GRETL_MPI_SUM &&
466 op != GRETL_MPI_PROD &&
467 op != GRETL_MPI_HCAT &&
468 op != GRETL_MPI_VCAT) {
469 return E_DATA;
470 }
471
472 err = gretl_comm_check(root, &id, &np);
473 if (err) {
474 return err;
475 }
476
477 if (id != root) {
478 /* send matrix dimensions to root */
479 fill_matrix_info(rc, sm);
480 err = mpi_send(rc, MI_LEN, mpi_int, root, TAG_MATRIX_INFO,
481 mpi_comm_world);
482 } else {
483 /* root: gather dimensions from other processes,
484 check for conformability and allocate storage
485 */
486 int i;
487
488 rows = malloc((np+1) * sizeof *rows);
489 cols = malloc((np+1) * sizeof *cols);
490 if (rows == NULL || cols == NULL) {
491 err = E_ALLOC;
492 }
493
494 /* initialize record of row/col maxima */
495 rows[np] = 0;
496 cols[np] = 0;
497
498 for (i=0; i<np && !err; i++) {
499 if (i == root) {
500 if (sm != NULL) {
501 rows[i] = sm->rows;
502 cols[i] = sm->cols;
503 } else {
504 rows[i] = cols[i] = 0;
505 }
506 } else {
507 err = mpi_recv(rc, MI_LEN, mpi_int, i, TAG_MATRIX_INFO,
508 mpi_comm_world, MPI_STATUS_IGNORE);
509 if (!err) {
510 rows[i] = rc[0];
511 cols[i] = rc[1];
512 }
513 }
514 }
515
516 if (!err) {
517 err = matrix_dims_check(rows, cols, np, op);
518 }
519
520 if (!err) {
521 err = matrix_reduce_alloc(rows, cols, np, op,
522 &rm, &val);
523 }
524 }
525
526 if (!err) {
527 if (id != root) {
528 /* send data to root */
529 int sendsize = rc[0] * rc[1];
530
531 if (sendsize > 0) {
532 err = mpi_send(sm->val, sendsize, mpi_double, root,
533 TAG_MATRIX_VAL, mpi_comm_world);
534 }
535 } else {
536 /* root gathers and processes data */
537 int i, recvsize, offset = 0;
538
539 for (i=0; i<np && !err; i++) {
540 recvsize = rows[i] * cols[i];
541 if (recvsize == 0) {
542 continue;
543 }
544 if (i == root) {
545 if (sm != NULL) {
546 err = matrix_reduce_step(rm, sm->val, recvsize, op,
547 &offset);
548 }
549 } else {
550 err = mpi_recv(val, recvsize, mpi_double, i,
551 TAG_MATRIX_VAL, mpi_comm_world,
552 MPI_STATUS_IGNORE);
553 if (!err) {
554 err = matrix_reduce_step(rm, val, recvsize, op,
555 &offset);
556 }
557 }
558 }
559 }
560 }
561
562 if (id == root) {
563 /* clean up */
564 free(val);
565 free(rows);
566 free(cols);
567 /* handle return value */
568 if (!err) {
569 *pm = rm;
570 } else {
571 gretl_matrix_free(rm);
572 }
573 }
574
575 if (!err && (opt & OPT_A)) {
576 err = gretl_matrix_bcast(pm, id, root);
577 }
578
579 return err;
580 }
581
gretl_array_mpi_reduce(gretl_array * sa,gretl_array ** pa,Gretl_MPI_Op op,int root)582 int gretl_array_mpi_reduce (gretl_array *sa,
583 gretl_array **pa,
584 Gretl_MPI_Op op,
585 int root)
586 {
587 gretl_array *a = NULL;
588 GretlType atype, etype;
589 int *nmvec = NULL;
590 int id, np, nm;
591 int i, j;
592 int err = 0;
593
594 if (op != GRETL_MPI_ACAT) {
595 return E_DATA;
596 }
597
598 err = gretl_comm_check(root, &id, &np);
599 if (err) {
600 return err;
601 }
602
603 atype = gretl_array_get_type(sa);
604 etype = gretl_array_get_content_type(sa);
605 nm = gretl_array_get_length(sa);
606
607 if (id != root) {
608 /* send size of our array to root */
609 err = mpi_send(&nm, 1, mpi_int, root, TAG_ARRAY_LEN,
610 mpi_comm_world);
611 } else {
612 /* root: gather and record array sizes from other processes;
613 allocate composite array
614 */
615 int ntotal = 0;
616
617 nmvec = malloc(np * sizeof *nmvec);
618 if (nmvec == NULL) {
619 err = E_ALLOC;
620 }
621 for (i=0; i<np && !err; i++) {
622 if (i == root) {
623 nmvec[i] = nm;
624 } else {
625 err = mpi_recv(&nmvec[i], 1, mpi_int, i, TAG_ARRAY_LEN,
626 mpi_comm_world, MPI_STATUS_IGNORE);
627 }
628 ntotal += nmvec[i];
629 }
630 if (!err) {
631 a = gretl_array_new(atype, ntotal, &err);
632 }
633 }
634
635 if (!err) {
636 void *data;
637
638 if (id != root) {
639 /* send our elements to root */
640 for (j=0; j<nm; j++) {
641 data = gretl_array_get_data(sa, j);
642 err = gretl_mpi_send(data, etype, root);
643 }
644 } else {
645 /* root: gather elements from other processes
646 and pack into big array
647 */
648 int k = 0;
649
650 for (i=0; i<np && !err; i++) {
651 for (j=0; j<nmvec[i] && !err; j++) {
652 if (i == root) {
653 data = gretl_array_get_element(sa, j, NULL, &err);
654 err = gretl_array_set_element(a, k++, data, etype, 1);
655 } else {
656 data = mpi_receive_element(i, etype, &err);
657 if (!err) {
658 err = gretl_array_set_data(a, k++, data);
659 }
660 }
661 }
662 }
663 }
664 }
665
666 if (id == root) {
667 free(nmvec);
668 /* handle return value */
669 if (!err) {
670 *pa = a;
671 } else {
672 gretl_array_destroy(a);
673 }
674 }
675
676 return err;
677 }
678
gretl_scalar_mpi_reduce(double x,double * xp,Gretl_MPI_Op op,int root,gretlopt opt)679 int gretl_scalar_mpi_reduce (double x,
680 double *xp,
681 Gretl_MPI_Op op,
682 int root,
683 gretlopt opt)
684 {
685 MPI_Op mpi_op;
686 int np, ret;
687
688 mpi_comm_size(mpi_comm_world, &np);
689 if (root < 0 || root >= np) {
690 return invalid_rank_error(root);
691 }
692
693 if (op == GRETL_MPI_SUM) {
694 mpi_op = mpi_sum;
695 } else if (op == GRETL_MPI_PROD) {
696 mpi_op = mpi_prod;
697 } else if (op == GRETL_MPI_MAX) {
698 mpi_op = mpi_max;
699 } else if (op == GRETL_MPI_MIN) {
700 mpi_op = mpi_min;
701 } else {
702 return E_DATA;
703 }
704
705 if (opt & OPT_A) {
706 ret = mpi_allreduce(&x, xp, 1, mpi_double, mpi_op,
707 mpi_comm_world);
708 } else {
709 ret = mpi_reduce(&x, xp, 1, mpi_double, mpi_op,
710 root, mpi_comm_world);
711 }
712
713 return ret;
714 }
715
maybe_date_matrix(gretl_matrix * m,int * rc)716 static void maybe_date_matrix (gretl_matrix *m, int *rc)
717 {
718 if (rc[3] >= 0 && rc[4] >= rc[3]) {
719 gretl_matrix_set_t1(m, rc[3]);
720 gretl_matrix_set_t2(m, rc[4]);
721 }
722 }
723
gretl_matrix_bcast(gretl_matrix ** pm,int id,int root)724 static int gretl_matrix_bcast (gretl_matrix **pm, int id, int root)
725 {
726 gretl_matrix *m = NULL;
727 int rc[MI_LEN];
728 int err = 0;
729
730 if (id == root) {
731 m = *pm;
732 fill_matrix_info(rc, m);
733 }
734
735 /* broadcast the matrix dimensions first */
736 err = mpi_bcast(rc, MI_LEN, mpi_int, root, mpi_comm_world);
737
738 if (!err && id != root) {
739 int r = rc[0];
740 int c = rc[1];
741 int cmplx = rc[2];
742
743 /* everyone but root needs to allocate space */
744 if (cmplx) {
745 *pm = m = gretl_cmatrix_new(r, c);
746 } else {
747 *pm = m = gretl_matrix_alloc(r, c);
748 }
749 if (m == NULL) {
750 return E_ALLOC;
751 } else {
752 maybe_date_matrix(m, rc);
753 }
754 }
755
756 if (!err) {
757 /* broadcast the matrix content */
758 int n = rc[0] * rc[1];
759
760 if (rc[2]) {
761 n *= 2;
762 }
763
764 /* FIXME we can get a hang here with 100% CPU if
765 a worker bombs out on bcast(); in that case
766 it seems that root's call never returns --
767 or maybe not before some looong time-out.
768 */
769 err = mpi_bcast(m->val, n, mpi_double, root,
770 mpi_comm_world);
771 }
772
773 if (err) {
774 gretl_mpi_error(&err);
775 }
776
777 return err;
778 }
779
gretl_series_bcast(double ** px,int len,int id,int root)780 static int gretl_series_bcast (double **px, int len,
781 int id, int root)
782 {
783 double *x = NULL;
784 int err = 0;
785
786 if (id == root) {
787 x = *px;
788 if (x == NULL) {
789 return E_DATA;
790 }
791 } else {
792 /* everyone else needs to allocate space */
793 x = malloc(len * sizeof *x);
794 if (x == NULL) {
795 return E_ALLOC;
796 }
797 }
798
799 if (!err) {
800 /* broadcast the series */
801 err = mpi_bcast(x, len, mpi_double, root, mpi_comm_world);
802 }
803
804 if (err) {
805 gretl_mpi_error(&err);
806 } else if (id != root) {
807 *px = x;
808 }
809
810 return err;
811 }
812
gretl_list_bcast(int ** plist,int id,int root)813 static int gretl_list_bcast (int **plist, int id, int root)
814 {
815 int *list = NULL;
816 int len = 0;
817 int err = 0;
818
819 if (id == root) {
820 list = *plist;
821 if (list == NULL) {
822 return E_DATA;
823 } else {
824 len = list[0];
825 }
826 }
827
828 /* broadcast the list length first */
829 err = mpi_bcast(&len, 1, mpi_int, root, mpi_comm_world);
830
831 if (!err && id != root) {
832 /* everyone but root needs to allocate space */
833 list = gretl_list_new(len);
834 if (list == NULL) {
835 return E_ALLOC;
836 }
837 }
838
839 if (!err) {
840 /* broadcast the list */
841 err = mpi_bcast(list, len + 1, mpi_int, root, mpi_comm_world);
842 }
843
844 if (err) {
845 gretl_mpi_error(&err);
846 } else if (id != root) {
847 *plist = list;
848 }
849
850 return err;
851 }
852
gretl_string_bcast(char ** pbuf,int id,int root)853 static int gretl_string_bcast (char **pbuf, int id, int root)
854 {
855 char *buf = NULL;
856 int bytes = 0;
857 int err = 0;
858
859 if (id == root) {
860 buf = *pbuf;
861 if (buf == NULL) {
862 return E_DATA;
863 } else {
864 bytes = strlen(buf) + 1;
865 }
866 }
867
868 /* broadcast the buffer size first */
869 err = mpi_bcast(&bytes, 1, mpi_int, root, mpi_comm_world);
870
871 if (!err && id != root) {
872 /* everyone but root needs to allocate space */
873 buf = calloc(bytes, 1);
874 if (buf == NULL) {
875 return E_ALLOC;
876 }
877 }
878
879 if (!err) {
880 /* broadcast the buffer */
881 err = mpi_bcast(buf, bytes, mpi_byte, root, mpi_comm_world);
882 }
883
884 if (err) {
885 gretl_mpi_error(&err);
886 } else if (id != root) {
887 *pbuf = buf;
888 }
889
890 return err;
891 }
892
gretl_scalar_bcast(double * px,int root)893 static int gretl_scalar_bcast (double *px, int root)
894 {
895 int err;
896
897 err = mpi_bcast(px, 1, mpi_double, root, mpi_comm_world);
898
899 if (err) {
900 gretl_mpi_error(&err);
901 }
902
903 return err;
904 }
905
gretl_int_bcast(int * pi,int root)906 static int gretl_int_bcast (int *pi, int root)
907 {
908 int err;
909
910 err = mpi_bcast(pi, 1, mpi_int, root, mpi_comm_world);
911
912 if (err) {
913 gretl_mpi_error(&err);
914 }
915
916 return err;
917 }
918
gretl_unsigned_bcast(unsigned int * pu,int root)919 static int gretl_unsigned_bcast (unsigned int *pu, int root)
920 {
921 int err;
922
923 err = mpi_bcast(pu, 1, mpi_unsigned, root, mpi_comm_world);
924
925 if (err) {
926 gretl_mpi_error(&err);
927 }
928
929 return err;
930 }
931
932 /* Compose a message indicating the type (and size, if
933 applicable) of a bundle member to be passed via MPI.
934 */
935
compose_msgbuf(char * buf,GretlType type,const char * key,int * size,void * data)936 static int compose_msgbuf (char *buf, GretlType type,
937 const char *key, int *size,
938 void *data)
939 {
940 int n = 0;
941
942 if (type == GRETL_TYPE_DOUBLE) {
943 n = sizeof(double);
944 } else if (type == GRETL_TYPE_INT) {
945 n = sizeof(int);
946 } else if (type == GRETL_TYPE_UNSIGNED) {
947 n = sizeof(unsigned int);
948 } else if (type == GRETL_TYPE_SERIES) {
949 n = *size;
950 } else if (type == GRETL_TYPE_STRING) {
951 char *s = data;
952 n = strlen(s) + 1;
953 } else if (type == GRETL_TYPE_LIST) {
954 int *list = data;
955 n = (list[0] + 1) * sizeof(int);
956 } else if (type == GRETL_TYPE_MATRIX ||
957 type == GRETL_TYPE_BUNDLE ||
958 type == GRETL_TYPE_ARRAY) {
959 n = 0;
960 } else {
961 return E_DATA;
962 }
963
964 sprintf(buf, "%d %s %d", type, key, n);
965 *size = n;
966
967 return 0;
968 }
969
970 /* Parse the information sent in @buf by root, regarding a
971 bundle member.
972 */
973
parse_msgbuf(const char * buf,char * key,GretlType * type,int * size)974 static int parse_msgbuf (const char *buf, char *key,
975 GretlType *type, int *size)
976 {
977 int t;
978
979 if (sscanf(buf, "%d %s %d", &t, key, size) != 3) {
980 return E_DATA;
981 } else {
982 *type = t;
983 return 0;
984 }
985 }
986
gretl_bundle_bcast(gretl_bundle ** pb,int id,int root)987 static int gretl_bundle_bcast (gretl_bundle **pb,
988 int id, int root)
989 {
990 gretl_bundle *b = NULL;
991 gretl_matrix *m = NULL;
992 gretl_array *a = NULL;
993 gretl_bundle *bsub = NULL;
994 gretl_array *keys = NULL;
995 double *x = NULL;
996 GretlType type;
997 char key[32];
998 void *data;
999 char msgbuf[64];
1000 int msglen;
1001 int size;
1002 int i, nk;
1003 int err = 0;
1004
1005 if (id == root) {
1006 b = *pb;
1007 nk = gretl_bundle_get_n_keys(b);
1008 keys = gretl_bundle_get_keys(b, &err);
1009 if (err) {
1010 return err;
1011 }
1012 }
1013
1014 /* broadcast the number of keys first */
1015 err = mpi_bcast(&nk, 1, mpi_int, root, mpi_comm_world);
1016
1017 if (!err && id != root) {
1018 /* everyone but root needs to start a bundle */
1019 b = gretl_bundle_new();
1020 if (b == NULL) {
1021 return E_ALLOC;
1022 }
1023 }
1024
1025 msglen = sizeof msgbuf;
1026
1027 for (i=0; i<nk && !err; i++) {
1028 /* loop across bundle keys */
1029 memset(msgbuf, 0, msglen);
1030 size = 0;
1031 data = NULL;
1032 m = NULL;
1033 a = NULL;
1034 bsub = NULL;
1035 x = NULL;
1036
1037 if (id == root) {
1038 const char *rkey = gretl_array_get_data(keys, i);
1039
1040 data = gretl_bundle_get_data(b, rkey, &type, &size, &err);
1041 if (!err) {
1042 compose_msgbuf(msgbuf, type, rkey, &size, data);
1043 if (type == GRETL_TYPE_MATRIX) {
1044 m = (gretl_matrix *) data;
1045 } else if (type == GRETL_TYPE_ARRAY) {
1046 a = (gretl_array *) data;
1047 } else if (type == GRETL_TYPE_BUNDLE) {
1048 bsub = (gretl_bundle *) data;
1049 } else if (type == GRETL_TYPE_SERIES) {
1050 x = (double *) data;
1051 }
1052 }
1053 }
1054 if (!err) {
1055 /* broadcast info on the bundle member */
1056 err = mpi_bcast(msgbuf, msglen, mpi_byte, root, mpi_comm_world);
1057 }
1058 if (!err) {
1059 err = parse_msgbuf(msgbuf, key, &type, &size);
1060 }
1061 if (err) {
1062 break;
1063 }
1064 if (type == GRETL_TYPE_MATRIX) {
1065 err = gretl_matrix_bcast(&m, id, root);
1066 if (!err && id != root) {
1067 err = gretl_bundle_donate_data(b, key, m, type, 0);
1068 }
1069 } else if (type == GRETL_TYPE_ARRAY) {
1070 err = gretl_array_bcast(&a, id, root);
1071 if (!err && id != root) {
1072 err = gretl_bundle_donate_data(b, key, a, type, 0);
1073 }
1074 } else if (type == GRETL_TYPE_BUNDLE) {
1075 err = gretl_bundle_bcast(&bsub, id, root);
1076 if (!err && id != root) {
1077 err = gretl_bundle_donate_data(b, key, bsub, type, 0);
1078 }
1079 } else if (type == GRETL_TYPE_SERIES) {
1080 err = gretl_series_bcast(&x, size, id, root);
1081 if (!err && id != root) {
1082 err = gretl_bundle_donate_data(b, key, x, type, size);
1083 }
1084 } else {
1085 /* scalar, string or list */
1086 if (id != root) {
1087 data = calloc(size, 1);
1088 if (data == NULL) {
1089 err = E_ALLOC;
1090 }
1091 }
1092 err = mpi_bcast(data, size, mpi_byte, root, mpi_comm_world);
1093 if (!err && id != root && data != NULL) {
1094 if (gretl_is_scalar_type(type)) {
1095 err = gretl_bundle_set_data(b, key, data, type, 0);
1096 free(data);
1097 } else {
1098 err = gretl_bundle_donate_data(b, key, data, type, 0);
1099 }
1100 }
1101 }
1102 }
1103
1104 gretl_array_destroy(keys);
1105
1106 mpi_barrier(mpi_comm_world);
1107
1108 if (err) {
1109 gretl_mpi_error(&err);
1110 } else if (id != root) {
1111 *pb = b;
1112 }
1113
1114 return err;
1115 }
1116
gretl_array_bcast(gretl_array ** pa,int id,int root)1117 static int gretl_array_bcast (gretl_array **pa, int id, int root)
1118 {
1119 gretl_array *a = NULL;
1120 GretlType type = 0;
1121 int nelem = 0;
1122 int st[2];
1123 int i, err = 0;
1124
1125 if (id == root) {
1126 a = *pa;
1127 st[0] = nelem = gretl_array_get_length(a);
1128 st[1] = type = gretl_array_get_type(a);
1129 }
1130
1131 /* broadcast the array size and type */
1132 err = mpi_bcast(st, 2, mpi_int, root, mpi_comm_world);
1133
1134 if (!err && id != root) {
1135 /* everyone but root needs to allocate an array */
1136 nelem = st[0];
1137 type = st[1];
1138 *pa = a = gretl_array_new(type, nelem, &err);
1139 }
1140
1141 for (i=0; i<nelem && !err; i++) {
1142 void *ptr, *data = NULL;
1143
1144 if (id == root) {
1145 data = gretl_array_get_data(a, i);
1146 }
1147 ptr = &data;
1148 if (type == GRETL_TYPE_MATRICES) {
1149 err = gretl_matrix_bcast(ptr, id, root);
1150 if (id != root) {
1151 data = *(gretl_matrix **) ptr;
1152 }
1153 } else if (type == GRETL_TYPE_STRINGS) {
1154 err = gretl_string_bcast(ptr, id, root);
1155 if (id != root) {
1156 data = *(char **) ptr;
1157 }
1158 } else if (type == GRETL_TYPE_BUNDLES) {
1159 err = gretl_bundle_bcast(ptr, id, root);
1160 if (id != root) {
1161 data = *(gretl_bundle **) ptr;
1162 }
1163 } else if (type == GRETL_TYPE_LISTS) {
1164 err = gretl_list_bcast(ptr, id, root);
1165 if (id != root) {
1166 data = *(int **) ptr;
1167 }
1168 } else if (type == GRETL_TYPE_ARRAYS) {
1169 err = gretl_array_bcast(ptr, id, root);
1170 if (id != root) {
1171 data = *(gretl_array **) ptr;
1172 }
1173 } else {
1174 /* ?? */
1175 err = E_TYPES;
1176 }
1177 if (!err && id != root) {
1178 gretl_array_set_data(a, i, data);
1179 }
1180 }
1181
1182 mpi_barrier(mpi_comm_world);
1183
1184 if (err) {
1185 gretl_mpi_error(&err);
1186 }
1187
1188 return err;
1189 }
1190
gretl_mpi_barrier(void)1191 int gretl_mpi_barrier (void)
1192 {
1193 return mpi_barrier(mpi_comm_world) != MPI_SUCCESS;
1194 }
1195
1196 /**
1197 * gretl_mpi_bcast:
1198 * @p: the location of the object to be broadcast.
1199 * @type: the type of the object to which @p points.
1200 * @root: the rank of the root process.
1201 *
1202 * Broadcasts the value referenced by @p to MPI_COMM_WORLD.
1203 * @type must be GRETL_TYPE_MATRIX, in which case
1204 * @p should be a (**gretl_matrix) pointer; GRETL_TYPE_BUNDLE,
1205 * in which case @p should be a (**gretl_bundle) pointer;
1206 * or GRETL_TYPE_DOUBLE, in which case @p should be a (*double)
1207 * pointer.
1208 *
1209 * Returns: 0 on successful completion, non-zero code otherwise.
1210 **/
1211
gretl_mpi_bcast(void * p,GretlType type,int root)1212 int gretl_mpi_bcast (void *p, GretlType type, int root)
1213 {
1214 int id, err;
1215
1216 err = gretl_comm_check(root, &id, NULL);
1217 if (err) {
1218 return err;
1219 }
1220
1221 if (type == GRETL_TYPE_DOUBLE) {
1222 return gretl_scalar_bcast((double *) p, root);
1223 } else if (type == GRETL_TYPE_INT) {
1224 return gretl_int_bcast((int *) p, root);
1225 } else if (type == GRETL_TYPE_UNSIGNED) {
1226 return gretl_unsigned_bcast((unsigned int *) p, root);
1227 } else if (type == GRETL_TYPE_MATRIX) {
1228 return gretl_matrix_bcast((gretl_matrix **) p, id, root);
1229 } else if (type == GRETL_TYPE_BUNDLE) {
1230 return gretl_bundle_bcast((gretl_bundle **) p, id, root);
1231 } else if (type == GRETL_TYPE_ARRAY) {
1232 return gretl_array_bcast((gretl_array **) p, id, root);
1233 } else if (type == GRETL_TYPE_STRING) {
1234 return gretl_string_bcast((char **) p, id, root);
1235 } else if (type == GRETL_TYPE_LIST) {
1236 return gretl_list_bcast((int **) p, id, root);
1237 } else {
1238 return E_DATA;
1239 }
1240 }
1241
1242 /* this "send" function is public, for convenience of callers
1243 such as the svm plugin
1244 */
1245
gretl_matrix_mpi_send(const gretl_matrix * m,int dest)1246 int gretl_matrix_mpi_send (const gretl_matrix *m, int dest)
1247 {
1248 int rc[MI_LEN];
1249 int err;
1250
1251 fill_matrix_info(rc, m);
1252
1253 err = mpi_send(rc, MI_LEN, mpi_int, dest, TAG_MATRIX_INFO,
1254 mpi_comm_world);
1255
1256 if (!err) {
1257 int n = m->rows * m->cols;
1258
1259 if (m->is_complex) {
1260 n *= 2;
1261 }
1262 err = mpi_send(m->val, n, mpi_double, dest, TAG_MATRIX_VAL,
1263 mpi_comm_world);
1264 }
1265
1266 if (err) {
1267 gretl_mpi_error(&err);
1268 }
1269
1270 return err;
1271 }
1272
gretl_string_send(char * s,int dest)1273 static int gretl_string_send (char *s, int dest)
1274 {
1275 int n = strlen(s) + 1;
1276 int err;
1277
1278 err = mpi_send(&n, 1, mpi_int, dest, TAG_STR_LEN,
1279 mpi_comm_world);
1280
1281 if (!err) {
1282 err = mpi_send(s, n, mpi_byte, dest, TAG_STR_VAL,
1283 mpi_comm_world);
1284 }
1285
1286 if (err) {
1287 gretl_mpi_error(&err);
1288 }
1289
1290 return err;
1291 }
1292
gretl_list_send(int * list,int dest)1293 static int gretl_list_send (int *list, int dest)
1294 {
1295 int n = list[0];
1296 int err;
1297
1298 err = mpi_send(&n, 1, mpi_int, dest, TAG_LIST_LEN,
1299 mpi_comm_world);
1300
1301 if (!err) {
1302 err = mpi_send(list, n+1, mpi_int, dest, TAG_LIST_VAL,
1303 mpi_comm_world);
1304 }
1305
1306 if (err) {
1307 gretl_mpi_error(&err);
1308 }
1309
1310 return err;
1311 }
1312
gretl_scalar_send(double * px,int dest)1313 static int gretl_scalar_send (double *px, int dest)
1314 {
1315 int err;
1316
1317 err = mpi_send(px, 1, mpi_double, dest, TAG_SCALAR_VAL,
1318 mpi_comm_world);
1319
1320 if (err) {
1321 gretl_mpi_error(&err);
1322 }
1323
1324 return err;
1325 }
1326
gretl_int_send(int * pi,int dest)1327 static int gretl_int_send (int *pi, int dest)
1328 {
1329 int err;
1330
1331 err = mpi_send(pi, 1, mpi_int, dest, TAG_INT_VAL,
1332 mpi_comm_world);
1333
1334 if (err) {
1335 gretl_mpi_error(&err);
1336 }
1337
1338 return err;
1339 }
1340
gretl_array_send(gretl_array * a,int dest)1341 static int gretl_array_send (gretl_array *a, int dest)
1342 {
1343 GretlType type = 0;
1344 int nelem = 0;
1345 int st[2];
1346 int i, err = 0;
1347
1348 st[0] = nelem = gretl_array_get_length(a);
1349 st[1] = type = gretl_array_get_type(a);
1350
1351 /* send the array size and type */
1352 err = mpi_send(st, 2, mpi_int, dest, TAG_ARRAY_INFO,
1353 mpi_comm_world);
1354
1355 for (i=0; i<nelem && !err; i++) {
1356 void *data = gretl_array_get_data(a, i);
1357
1358 if (type == GRETL_TYPE_MATRICES) {
1359 err = gretl_matrix_mpi_send(data, dest);
1360 } else if (type == GRETL_TYPE_STRINGS) {
1361 err = gretl_string_send(data, dest);
1362 } else if (type == GRETL_TYPE_BUNDLES) {
1363 err = gretl_bundle_send(data, dest);
1364 } else if (type == GRETL_TYPE_LISTS) {
1365 err = gretl_list_send(data, dest);
1366 } else if (type == GRETL_TYPE_ARRAYS) {
1367 err = gretl_array_send(data, dest);
1368 }
1369 }
1370
1371 if (err) {
1372 gretl_mpi_error(&err);
1373 }
1374
1375 return err;
1376 }
1377
1378 /**
1379 * gretl_mpi_send:
1380 * @p: pointer to the object to be sent.
1381 * @type: the type of the object.
1382 * @dest: the MPI rank of the destination.
1383 *
1384 * Sends the value referenced by @p, of gretl type @type, to the
1385 * MPI process with rank @dest.
1386 *
1387 * Returns: 0 on successful completion, non-zero code otherwise.
1388 **/
1389
gretl_mpi_send(void * p,GretlType type,int dest)1390 int gretl_mpi_send (void *p, GretlType type, int dest)
1391 {
1392 int np;
1393
1394 mpi_comm_size(mpi_comm_world, &np);
1395 if (dest < 0 || dest >= np) {
1396 return invalid_rank_error(dest);
1397 }
1398
1399 if (type == GRETL_TYPE_DOUBLE) {
1400 return gretl_scalar_send((double *) p, dest);
1401 } else if (type == GRETL_TYPE_INT) {
1402 return gretl_int_send((int *) p, dest);
1403 } else if (type == GRETL_TYPE_MATRIX) {
1404 return gretl_matrix_mpi_send((gretl_matrix *) p, dest);
1405 } else if (type == GRETL_TYPE_BUNDLE) {
1406 return gretl_bundle_send((gretl_bundle *) p, dest);
1407 } else if (type == GRETL_TYPE_ARRAY) {
1408 return gretl_array_send((gretl_array *) p, dest);
1409 } else if (type == GRETL_TYPE_STRING) {
1410 return gretl_string_send((char *) p, dest);
1411 } else if (type == GRETL_TYPE_LIST) {
1412 return gretl_list_send((int *) p, dest);
1413 } else {
1414 return E_TYPES;
1415 }
1416 }
1417
gretl_matrix_mpi_receive(int source,int * err)1418 gretl_matrix *gretl_matrix_mpi_receive (int source,
1419 int *err)
1420 {
1421 gretl_matrix *m = NULL;
1422 int rc[MI_LEN];
1423
1424 *err = mpi_recv(rc, MI_LEN, mpi_int, source, TAG_MATRIX_INFO,
1425 mpi_comm_world, MPI_STATUS_IGNORE);
1426
1427 if (!*err) {
1428 int r = rc[0];
1429 int c = rc[1];
1430 int cmplx = rc[2];
1431 int n = r * c;
1432
1433 if (cmplx) {
1434 m = gretl_cmatrix_new(r, c);
1435 n *= 2;
1436 } else {
1437 m = gretl_matrix_alloc(r, c);
1438 }
1439
1440 if (m == NULL) {
1441 *err = E_ALLOC;
1442 return NULL;
1443 } else {
1444 *err = mpi_recv(m->val, n, mpi_double, source,
1445 TAG_MATRIX_VAL, mpi_comm_world,
1446 MPI_STATUS_IGNORE);
1447 if (*err) {
1448 maybe_date_matrix(m, rc);
1449 }
1450 }
1451 }
1452
1453 if (*err) {
1454 gretl_mpi_error(err);
1455 }
1456
1457 return m;
1458 }
1459
gretl_matrix_mpi_fill(gretl_matrix ** pm,int source)1460 int gretl_matrix_mpi_fill (gretl_matrix **pm, int source)
1461 {
1462 int rc[MI_LEN];
1463 int err;
1464
1465 if (pm == NULL) {
1466 return E_DATA;
1467 }
1468
1469 err = mpi_recv(rc, MI_LEN, mpi_int, source, TAG_MATRIX_INFO,
1470 mpi_comm_world, MPI_STATUS_IGNORE);
1471
1472 if (!err) {
1473 gretl_matrix *m = *pm;
1474 int r = rc[0];
1475 int c = rc[1];
1476 int cmplx = rc[2];
1477 int n = r * c;
1478
1479 if (m == NULL) {
1480 if (cmplx) {
1481 m = gretl_cmatrix_new(r, c);
1482 n *= 2;
1483 } else {
1484 m = gretl_matrix_alloc(r, c);
1485 }
1486 if (m == NULL) {
1487 err = E_ALLOC;
1488 } else {
1489 *pm = m;
1490 }
1491 } else if (m->rows != r || m->cols != c ||
1492 m->is_complex != cmplx) {
1493 err = E_NONCONF;
1494 }
1495
1496 if (!err) {
1497 err = mpi_recv(m->val, n, mpi_double, source,
1498 TAG_MATRIX_VAL, mpi_comm_world,
1499 MPI_STATUS_IGNORE);
1500 if (err) {
1501 maybe_date_matrix(m, rc);
1502 }
1503 }
1504 }
1505
1506 if (err) {
1507 gretl_mpi_error(&err);
1508 }
1509
1510 return err;
1511 }
1512
gretl_list_receive(int source,int * err)1513 static int *gretl_list_receive (int source, int *err)
1514 {
1515 int *list = NULL;
1516 int n;
1517
1518 *err = mpi_recv(&n, 1, mpi_int, source, TAG_LIST_LEN,
1519 mpi_comm_world, MPI_STATUS_IGNORE);
1520
1521 if (!*err) {
1522 list = gretl_list_new(n);
1523 if (list == NULL) {
1524 *err = E_ALLOC;
1525 } else {
1526 *err = mpi_recv(list, n+1, mpi_int, source,
1527 TAG_LIST_VAL, mpi_comm_world,
1528 MPI_STATUS_IGNORE);
1529 }
1530 }
1531
1532 return list;
1533 }
1534
gretl_string_receive(int source,int * err)1535 static char *gretl_string_receive (int source, int *err)
1536 {
1537 char *s = NULL;
1538 int n;
1539
1540 *err = mpi_recv(&n, 1, mpi_int, source, TAG_STR_LEN,
1541 mpi_comm_world, MPI_STATUS_IGNORE);
1542
1543 if (!*err) {
1544 s = calloc(n, 1);
1545 if (s == NULL) {
1546 *err = E_ALLOC;
1547 } else {
1548 *err = mpi_recv(s, n, mpi_byte, source,
1549 TAG_STR_VAL, mpi_comm_world,
1550 MPI_STATUS_IGNORE);
1551 }
1552 }
1553
1554 return s;
1555 }
1556
gretl_scalar_receive(int source,int * err)1557 static double gretl_scalar_receive (int source, int *err)
1558 {
1559 double x = NADBL;
1560
1561 *err = mpi_recv(&x, 1, mpi_double, source, TAG_SCALAR_VAL,
1562 mpi_comm_world, MPI_STATUS_IGNORE);
1563 if (*err) {
1564 gretl_mpi_error(err);
1565 }
1566
1567 return x;
1568 }
1569
gretl_int_receive(int source,int * err)1570 static double gretl_int_receive (int source, int *err)
1571 {
1572 int i = 0;
1573
1574 *err = mpi_recv(&i, 1, mpi_int, source, TAG_INT_VAL,
1575 mpi_comm_world, MPI_STATUS_IGNORE);
1576 if (*err) {
1577 gretl_mpi_error(err);
1578 }
1579
1580 return i;
1581 }
1582
gretl_array_receive(int source,int * err)1583 static gretl_array *gretl_array_receive (int source, int *err)
1584 {
1585 gretl_array *a = NULL;
1586 GretlType type = 0;
1587 int nelem = 0;
1588 int st[2];
1589 int i;
1590
1591 /* get the array size and type */
1592 *err = mpi_recv(st, 2, mpi_int, source, TAG_ARRAY_INFO,
1593 mpi_comm_world, MPI_STATUS_IGNORE);
1594
1595 if (!*err) {
1596 nelem = st[0];
1597 type = st[1];
1598 a = gretl_array_new(type, nelem, err);
1599 }
1600
1601 for (i=0; i<nelem && !*err; i++) {
1602 void *data = NULL;
1603
1604 if (type == GRETL_TYPE_MATRICES) {
1605 data = gretl_matrix_mpi_receive(source, err);
1606 } else if (type == GRETL_TYPE_STRINGS) {
1607 data = gretl_string_receive(source, err);
1608 } else if (type == GRETL_TYPE_BUNDLES) {
1609 data = gretl_bundle_receive(source, err);
1610 } else if (type == GRETL_TYPE_LISTS) {
1611 data = gretl_list_receive(source, err);
1612 } else if (type == GRETL_TYPE_ARRAYS) {
1613 data = gretl_array_receive(source, err);
1614 }
1615 if (data != NULL) {
1616 gretl_array_set_data(a, i, data);
1617 }
1618 }
1619
1620 return a;
1621 }
1622
type_from_status(MPI_Status * status)1623 static GretlType type_from_status (MPI_Status *status)
1624 {
1625 if (status->MPI_TAG == TAG_MATRIX_INFO) {
1626 return GRETL_TYPE_MATRIX;
1627 } else if (status->MPI_TAG == TAG_BUNDLE_SIZE) {
1628 return GRETL_TYPE_BUNDLE;
1629 } else if (status->MPI_TAG == TAG_ARRAY_INFO) {
1630 return GRETL_TYPE_ARRAY;
1631 } else if (status->MPI_TAG == TAG_SCALAR_VAL) {
1632 return GRETL_TYPE_DOUBLE;
1633 } else if (status->MPI_TAG == TAG_INT_VAL) {
1634 return GRETL_TYPE_INT;
1635 } else if (status->MPI_TAG == TAG_STR_LEN) {
1636 return GRETL_TYPE_STRING;
1637 } else if (status->MPI_TAG == TAG_LIST_LEN) {
1638 return GRETL_TYPE_LIST;
1639 } else {
1640 return GRETL_TYPE_NONE;
1641 }
1642 }
1643
gretl_mpi_receive(int source,GretlType * ptype,int * err)1644 void *gretl_mpi_receive (int source, GretlType *ptype, int *err)
1645 {
1646 static double x;
1647 static int k;
1648 void *ret = NULL;
1649 MPI_Status status;
1650 int np;
1651
1652 mpi_comm_size(mpi_comm_world, &np);
1653 if (source < 0 || source >= np) {
1654 *err = invalid_rank_error(source);
1655 return NULL;
1656 }
1657
1658 /* check for the type of thing of offer from @source */
1659 mpi_probe(source, MPI_ANY_TAG, mpi_comm_world, &status);
1660 *ptype = type_from_status(&status);
1661
1662 if (*ptype == GRETL_TYPE_DOUBLE) {
1663 x = gretl_scalar_receive(source, err);
1664 ret = &x;
1665 } else if (*ptype == GRETL_TYPE_INT) {
1666 k = gretl_int_receive(source, err);
1667 ret = &k;
1668 } else if (*ptype == GRETL_TYPE_MATRIX) {
1669 ret = gretl_matrix_mpi_receive(source, err);
1670 } else if (*ptype == GRETL_TYPE_BUNDLE) {
1671 ret = gretl_bundle_receive(source, err);
1672 } else if (*ptype == GRETL_TYPE_ARRAY) {
1673 ret = gretl_array_receive(source, err);
1674 } else if (*ptype == GRETL_TYPE_STRING) {
1675 ret = gretl_string_receive(source, err);
1676 } else if (*ptype == GRETL_TYPE_LIST) {
1677 ret = gretl_list_receive(source, err);
1678 } else {
1679 *err = E_TYPES;
1680 }
1681
1682 return ret;
1683 }
1684
mpi_receive_element(int source,GretlType etype,int * err)1685 static void *mpi_receive_element (int source, GretlType etype,
1686 int *err)
1687 {
1688 void *ret = NULL;
1689 MPI_Status status;
1690 GretlType srctype;
1691
1692 mpi_probe(source, MPI_ANY_TAG, mpi_comm_world, &status);
1693 srctype = type_from_status(&status);
1694
1695 if (srctype != etype) {
1696 *err = E_TYPES;
1697 } else if (etype == GRETL_TYPE_MATRIX) {
1698 ret = gretl_matrix_mpi_receive(source, err);
1699 } else if (etype == GRETL_TYPE_BUNDLE) {
1700 ret = gretl_bundle_receive(source, err);
1701 } else if (etype == GRETL_TYPE_ARRAY) {
1702 ret = gretl_array_receive(source, err);
1703 } else if (etype == GRETL_TYPE_STRING) {
1704 ret = gretl_string_receive(source, err);
1705 } else if (etype == GRETL_TYPE_LIST) {
1706 ret = gretl_list_receive(source, err);
1707 }
1708
1709 return ret;
1710 }
1711
1712 /* "new"-style send/receive for bundles: we do everything
1713 in memory */
1714
gretl_series_send(double * x,int n,int dest)1715 static int gretl_series_send (double *x, int n, int dest)
1716 {
1717 int err;
1718
1719 err = mpi_send(&n, 1, mpi_int, dest, TAG_SERIES_LEN,
1720 mpi_comm_world);
1721
1722 if (!err) {
1723 err = mpi_send(x, n, mpi_double, dest, TAG_SERIES_VAL,
1724 mpi_comm_world);
1725 }
1726
1727 if (err) {
1728 gretl_mpi_error(&err);
1729 }
1730
1731 return err;
1732 }
1733
gretl_series_receive(int source,int * err)1734 static double *gretl_series_receive (int source, int *err)
1735 {
1736 double *x = NULL;
1737 int n;
1738
1739 *err = mpi_recv(&n, 1, mpi_int, source, TAG_SERIES_LEN,
1740 mpi_comm_world, MPI_STATUS_IGNORE);
1741
1742 if (!*err) {
1743 x = malloc(n * sizeof *x);
1744 if (x == NULL) {
1745 *err = E_ALLOC;
1746 return NULL;
1747 } else {
1748 *err = mpi_recv(x, n, mpi_double, source,
1749 TAG_SERIES_VAL, mpi_comm_world,
1750 MPI_STATUS_IGNORE);
1751 }
1752 }
1753
1754 if (*err) {
1755 gretl_mpi_error(err);
1756 }
1757
1758 return x;
1759 }
1760
gretl_bundle_send(gretl_bundle * b,int dest)1761 static int gretl_bundle_send (gretl_bundle *b, int dest)
1762 {
1763 gretl_array *keys = NULL;
1764 GretlType type;
1765 void *data;
1766 char msgbuf[64];
1767 int msglen = sizeof msgbuf;
1768 int size;
1769 int i, nk;
1770 int err = 0;
1771
1772 nk = gretl_bundle_get_n_keys(b);
1773 if (nk > 0) {
1774 keys = gretl_bundle_get_keys(b, &err);
1775 if (err) {
1776 return err;
1777 }
1778 }
1779
1780 /* send the number of keys first */
1781 err = mpi_send(&nk, 1, mpi_int, dest, TAG_BUNDLE_SIZE,
1782 mpi_comm_world);
1783
1784 if (nk == 0) {
1785 /* the bundle is empty */
1786 return 0;
1787 }
1788
1789 for (i=0; i<nk && !err; i++) {
1790 /* loop across bundle keys */
1791 const char *key = gretl_array_get_data(keys, i);
1792
1793 data = gretl_bundle_get_data(b, key, &type, &size, &err);
1794 if (!err) {
1795 /* send info on the bundle member */
1796 compose_msgbuf(msgbuf, type, key, &size, data);
1797 err = mpi_send(msgbuf, msglen, mpi_byte, dest,
1798 TAG_BMEMB_INFO, mpi_comm_world);
1799 }
1800 if (err) {
1801 break;
1802 }
1803 if (type == GRETL_TYPE_MATRIX) {
1804 err = gretl_matrix_mpi_send(data, dest);
1805 } else if (type == GRETL_TYPE_ARRAY) {
1806 err = gretl_array_send(data, dest);
1807 } else if (type == GRETL_TYPE_BUNDLE) {
1808 err = gretl_bundle_send(data, dest);
1809 } else if (type == GRETL_TYPE_SERIES) {
1810 err = gretl_series_send(data, size, dest);
1811 } else if (type == GRETL_TYPE_STRING) {
1812 err = gretl_string_send(data, dest);
1813 } else if (type == GRETL_TYPE_LIST) {
1814 err = gretl_list_send(data, dest);
1815 } else if (type == GRETL_TYPE_DOUBLE) {
1816 err = gretl_scalar_send(data, dest);
1817 } else if (type == GRETL_TYPE_INT) {
1818 err = gretl_int_send(data, dest);
1819 }
1820 }
1821
1822 gretl_array_destroy(keys);
1823
1824 if (err) {
1825 gretl_mpi_error(&err);
1826 }
1827
1828 return err;
1829 }
1830
gretl_bundle_receive(int source,int * err)1831 gretl_bundle *gretl_bundle_receive (int source, int *err)
1832 {
1833 gretl_bundle *b = NULL;
1834 char key[32];
1835 char msgbuf[64];
1836 int msglen = sizeof msgbuf;
1837 int size;
1838 int i, nk;
1839
1840 /* get the number of keys first */
1841 *err = mpi_recv(&nk, 1, mpi_int, source, TAG_BUNDLE_SIZE,
1842 mpi_comm_world, MPI_STATUS_IGNORE);
1843
1844 if (!*err) {
1845 b = gretl_bundle_new();
1846 if (b == NULL) {
1847 *err = E_ALLOC;
1848 return NULL;
1849 }
1850 }
1851
1852 if (nk == 0) {
1853 /* the bundle is empty */
1854 return b;
1855 }
1856
1857 for (i=0; i<nk && !*err; i++) {
1858 /* loop across bundle keys */
1859 GretlType type = 0;
1860 void *data = NULL;
1861 double x;
1862
1863 memset(msgbuf, 0, msglen);
1864 size = 0;
1865
1866 /* get info on bundle member @i */
1867 *err = mpi_recv(msgbuf, msglen, mpi_byte, source, TAG_BMEMB_INFO,
1868 mpi_comm_world, MPI_STATUS_IGNORE);
1869 if (!*err) {
1870 *err = parse_msgbuf(msgbuf, key, &type, &size);
1871 }
1872 if (*err) {
1873 break;
1874 }
1875 if (type == GRETL_TYPE_DOUBLE) {
1876 x = gretl_scalar_receive(source, err);
1877 data = &x;
1878 } else if (type == GRETL_TYPE_INT) {
1879 x = gretl_int_receive(source, err);
1880 data = &x;
1881 } else if (type == GRETL_TYPE_MATRIX) {
1882 data = gretl_matrix_mpi_receive(source, err);
1883 } else if (type == GRETL_TYPE_ARRAY) {
1884 data = gretl_array_receive(source, err);
1885 } else if (type == GRETL_TYPE_BUNDLE) {
1886 data = gretl_bundle_receive(source, err);
1887 } else if (type == GRETL_TYPE_SERIES) {
1888 data = gretl_series_receive(source, err);
1889 } else if (type == GRETL_TYPE_LIST) {
1890 data = gretl_list_receive(source, err);
1891 } else if (type == GRETL_TYPE_STRING) {
1892 data = gretl_string_receive(source, err);
1893 }
1894 if (!*err) {
1895 *err = gretl_bundle_donate_data(b, key, data, type, size);
1896 }
1897 }
1898
1899 if (*err) {
1900 gretl_mpi_error(err);
1901 }
1902
1903 return b;
1904 }
1905
fill_tmp(double * restrict tmp,const gretl_matrix * m,int nr,int cmplx,int * offset)1906 static void fill_tmp (double * restrict tmp,
1907 const gretl_matrix *m,
1908 int nr, int cmplx,
1909 int *offset)
1910 {
1911 double x;
1912 double complex z;
1913 int imin = *offset;
1914 int imax = imin + nr;
1915 int i, j, k = 0;
1916
1917 for (j=0; j<m->cols; j++) {
1918 for (i=imin; i<imax; i++) {
1919 if (cmplx) {
1920 z = gretl_cmatrix_get(m, i, j);
1921 tmp[k++] = creal(z);
1922 tmp[k++] = cimag(z);
1923 } else {
1924 x = gretl_matrix_get(m, i, j);
1925 tmp[k++] = x;
1926 }
1927 }
1928 }
1929
1930 *offset += nr;
1931 }
1932
scatter_to_self(int * rc,double * val,gretl_matrix ** pm)1933 static int scatter_to_self (int *rc, double *val,
1934 gretl_matrix **pm)
1935 {
1936 int cmplx = rc[2];
1937 int err = 0;
1938
1939 if (cmplx) {
1940 *pm = gretl_cmatrix_new(rc[0], rc[1]);
1941 } else {
1942 *pm = gretl_matrix_alloc(rc[0], rc[1]);
1943 }
1944
1945 if (*pm == NULL) {
1946 err = E_ALLOC;
1947 } else {
1948 size_t n = rc[0] * rc[1] * sizeof *val;
1949
1950 if (cmplx) {
1951 n *= 2;
1952 }
1953 memcpy((*pm)->val, val, n);
1954 }
1955
1956 return err;
1957 }
1958
matsplit_rule(int n,int np,int * k1,int * n1,int * k2,int * n2)1959 static void matsplit_rule (int n, int np, int *k1, int *n1,
1960 int *k2, int *n2)
1961 {
1962 if (n <= np) {
1963 *k1 = n;
1964 *n1 = 1;
1965 *k2 = np - n;
1966 *n2 = 0;
1967 } else if (n % np == 0) {
1968 *k1 = np;
1969 *n1 = n / np;
1970 *k2 = 0;
1971 *n2 = 0;
1972 } else {
1973 int a = ceil(1.0e-12 + n / (double) np);
1974 int head = n % np;
1975
1976 *k1 = head;
1977 *n1 = a;
1978 *k2 = np - head;
1979 *n2 = a - 1;
1980 }
1981 }
1982
gretl_matrix_mpi_scatter(const gretl_matrix * m,gretl_matrix ** recvm,Gretl_MPI_Op op,int root)1983 int gretl_matrix_mpi_scatter (const gretl_matrix *m,
1984 gretl_matrix **recvm,
1985 Gretl_MPI_Op op,
1986 int root)
1987 {
1988 double *tmp = NULL;
1989 int id, np;
1990 int rc[MI_LEN] = {0};
1991 int err = 0;
1992
1993 mpi_comm_rank(mpi_comm_world, &id);
1994 mpi_comm_size(mpi_comm_world, &np);
1995
1996 if (root < 0 || root >= np) {
1997 return invalid_rank_error(root);
1998 }
1999
2000 if (id == root) {
2001 int cmplx = m->is_complex;
2002 int k1, n1, k2, n2;
2003 int i, n;
2004
2005 if (op == GRETL_MPI_VSPLIT) {
2006 /* scatter by rows */
2007 int offset = 0;
2008
2009 matsplit_rule(m->rows, np, &k1, &n1, &k2, &n2);
2010
2011 n = n1 * m->cols;
2012 rc[0] = n1;
2013 rc[1] = m->cols;
2014 if (cmplx) {
2015 rc[2] = 1;
2016 n *= 2;
2017 }
2018
2019 /* we'll need a working buffer */
2020 tmp = malloc(n * sizeof *tmp);
2021 if (tmp == NULL) {
2022 err = E_ALLOC;
2023 }
2024
2025 for (i=0; i<np; i++) {
2026 if (i == k1) {
2027 rc[0] = n2;
2028 n = n2 * m->cols;
2029 if (cmplx) n *= 2;
2030 }
2031 if (rc[0] > 0) {
2032 fill_tmp(tmp, m, rc[0], cmplx, &offset);
2033 }
2034 if (i == root) {
2035 err = scatter_to_self(rc, tmp, recvm);
2036 } else {
2037 err = mpi_send(rc, MI_LEN, mpi_int, i, TAG_MATRIX_INFO,
2038 mpi_comm_world);
2039 err = mpi_send(tmp, n, mpi_double, i, TAG_MATRIX_VAL,
2040 mpi_comm_world);
2041 }
2042 }
2043 } else {
2044 /* scatter by columns */
2045 double *val = m->val;
2046
2047 matsplit_rule(m->cols, np, &k1, &n1, &k2, &n2);
2048
2049 n = n1 * m->rows;
2050 if (cmplx) {
2051 n *= 2;
2052 }
2053 fill_matrix_info(rc, m);
2054 rc[1] = n1;
2055
2056 for (i=0; i<np; i++) {
2057 if (i == k1) {
2058 rc[1] = n2;
2059 n = n2 * m->rows;
2060 if (cmplx) n *= 2;
2061 }
2062 if (i == root) {
2063 err = scatter_to_self(rc, val, recvm);
2064 } else {
2065 err = mpi_send(rc, MI_LEN, mpi_int, i, TAG_MATRIX_INFO,
2066 mpi_comm_world);
2067 err = mpi_send(val, n, mpi_double, i, TAG_MATRIX_VAL,
2068 mpi_comm_world);
2069 }
2070 /* advance the read pointer */
2071 val += n;
2072 }
2073 }
2074 } else {
2075 /* non-root processes get their share-out */
2076 *recvm = gretl_matrix_mpi_receive(root, &err);
2077 }
2078
2079 if (id == root) {
2080 free(tmp);
2081 }
2082
2083 return err;
2084 }
2085
2086 /* MPI timer */
2087
gretl_mpi_time(void)2088 double gretl_mpi_time (void)
2089 {
2090 return mpi_wtime();
2091 }
2092
2093 /* end MPI timer */
2094
gretl_mpi_rank(void)2095 int gretl_mpi_rank (void)
2096 {
2097 int id = -1;
2098
2099 if (gretl_mpi_initialized()) {
2100 mpi_comm_rank(mpi_comm_world, &id);
2101 }
2102
2103 return id;
2104 }
2105
gretl_mpi_n_processes(void)2106 int gretl_mpi_n_processes (void)
2107 {
2108 int np = 0;
2109
2110 if (gretl_mpi_initialized()) {
2111 mpi_comm_size(mpi_comm_world, &np);
2112 }
2113
2114 return np;
2115 }
2116
gretl_mpi_initialized(void)2117 int gretl_mpi_initialized (void)
2118 {
2119 static int ret = -1;
2120
2121 if (!gretl_MPI_initted) {
2122 return 0;
2123 }
2124
2125 if (ret < 0) {
2126 mpi_initialized(&ret);
2127 ret = ret > 0;
2128 }
2129
2130 return ret;
2131 }
2132
2133 /* Experimental: implementation of the functions mwrite() and
2134 mread() via shared memory. This implementation is invoked
2135 when the matrix filename has suffix ".shm". The filename
2136 should not have any path component, as in "mymatrix.shm".
2137
2138 This is intended for fast transport between gretl or gretlcli
2139 and gretlmpi when an "mpi block" is used. In that context a
2140 single MPI process (presumably that with rank 0) must take
2141 exclusive charge of the shared-memory transaction. Unlike
2142 mwrite/mread using regular files, the matrix that is written
2143 into mapped memory by mwrite is cleared on the first access
2144 by mread so any subsequent attempts to read it will fail.
2145 */
2146
2147 #define MATRIX_ID "gretl_matrix"
2148 #define IDLEN 12 /* strlen(MATRIX_ID) */
2149
2150 #ifdef G_OS_WIN32
2151
2152 #include "gretl_win32.h"
2153
shm_write_matrix(const gretl_matrix * m,const char * fname)2154 int shm_write_matrix (const gretl_matrix *m,
2155 const char *fname)
2156 {
2157 int vsize = m->rows * m->cols * sizeof *m->val;
2158 int msize = IDLEN + 2 * sizeof(int) + vsize;
2159 HANDLE diskfile, mapfile;
2160 void *ptr = NULL;
2161 gchar *diskname;
2162 int err = 0;
2163
2164 diskname = g_strdup_printf("%s%s", gretl_dotdir(), fname);
2165
2166 diskfile = CreateFile(diskname,
2167 GENERIC_READ | GENERIC_WRITE,
2168 FILE_SHARE_DELETE | FILE_SHARE_READ,
2169 NULL,
2170 CREATE_ALWAYS,
2171 FILE_ATTRIBUTE_NORMAL,
2172 NULL);
2173
2174 if (diskfile == INVALID_HANDLE_VALUE) {
2175 fprintf(stderr, "mwrite: CreateFile failed for '%s'\n", diskname);
2176 win_print_last_error();
2177 err = E_FOPEN;
2178 }
2179
2180 if (!err) {
2181 mapfile = CreateFileMapping(diskfile, /* existing file handle */
2182 NULL, /* default security */
2183 PAGE_READWRITE, /* read/write access */
2184 0, /* max size (high DWORD) */
2185 msize, /* max size (low DWORD) */
2186 NULL); /* no name required */
2187 if (mapfile == NULL) {
2188 fprintf(stderr, "mwrite: CreateFileMapping failed for '%s'\n", diskname);
2189 win_print_last_error();
2190 err = E_FOPEN;
2191 }
2192 }
2193
2194 if (!err) {
2195 ptr = MapViewOfFile(mapfile, FILE_MAP_ALL_ACCESS, 0, 0, msize);
2196 if (ptr == NULL) {
2197 fprintf(stderr, "mwrite: MapViewOfFile failed\n");
2198 err = E_ALLOC;
2199 }
2200 }
2201
2202 if (!err) {
2203 void *pos = ptr;
2204
2205 memcpy(pos, MATRIX_ID, IDLEN);
2206 pos += IDLEN;
2207 memcpy(pos, &m->rows, sizeof m->rows);
2208 pos += sizeof m->rows;
2209 memcpy(pos, &m->cols, sizeof m->cols);
2210 pos += sizeof m->cols;
2211 memcpy(pos, m->val, vsize);
2212 }
2213
2214 if (ptr != NULL) {
2215 UnmapViewOfFile(ptr);
2216 }
2217 if (mapfile != NULL) {
2218 CloseHandle(mapfile);
2219 }
2220 if (diskfile != NULL) {
2221 CloseHandle(diskfile);
2222 }
2223
2224 g_free(diskname);
2225
2226 return err;
2227 }
2228
shm_read_matrix(const char * fname,int finalize,int * err)2229 gretl_matrix *shm_read_matrix (const char *fname,
2230 int finalize,
2231 int *err)
2232 {
2233 gretl_matrix *m = NULL;
2234 char buf[IDLEN] = {0};
2235 int isize;
2236 HANDLE diskfile, mapfile;
2237 void *ptr = NULL;
2238 gchar *diskname;
2239 int r, c, fd;
2240
2241 /* initial object size */
2242 isize = IDLEN + 2 * sizeof(int);
2243
2244 diskname = g_strdup_printf("%s%s", gretl_dotdir(), fname);
2245
2246 diskfile = CreateFile(diskname,
2247 GENERIC_READ | GENERIC_WRITE,
2248 0, /* no sharing here */
2249 NULL, /* security */
2250 OPEN_EXISTING, /* must already exist */
2251 FILE_ATTRIBUTE_NORMAL,
2252 NULL);
2253
2254 if (diskfile == INVALID_HANDLE_VALUE) {
2255 fprintf(stderr, "mread: CreateFile failed for '%s'\n", diskfile);
2256 win_print_last_error();
2257 *err = E_FOPEN;
2258 }
2259
2260 if (!*err) {
2261 mapfile = CreateFileMapping(diskfile,
2262 NULL, /* default security */
2263 PAGE_READWRITE, /* read/write access */
2264 0, /* accept existing size */
2265 0, /* ditto */
2266 NULL); /* no name required */
2267 if (mapfile == NULL) {
2268 fprintf(stderr, "mread: CreateFileMapping failed for '%s'\n", diskname);
2269 win_print_last_error();
2270 *err = E_FOPEN;
2271 }
2272 }
2273
2274 if (!*err) {
2275 ptr = MapViewOfFile(mapfile, FILE_MAP_ALL_ACCESS, 0, 0, isize);
2276 if (ptr == NULL) {
2277 fprintf(stderr, "mread: MapViewOfFile failed\n");
2278 *err = E_ALLOC;
2279 }
2280 }
2281
2282 if (ptr != NULL) {
2283 void *pos = ptr;
2284
2285 memcpy(buf, pos, IDLEN);
2286 if (memcmp(buf, MATRIX_ID, IDLEN) == 0) {
2287 pos += IDLEN;
2288 memcpy(&r, pos, sizeof r);
2289 pos += sizeof r;
2290 memcpy(&c, pos, sizeof c);
2291 m = gretl_matrix_alloc(r, c);
2292 if (m == NULL) {
2293 *err = E_ALLOC;
2294 }
2295 }
2296 }
2297
2298 if (m != NULL) {
2299 int vsize = r * c * sizeof *m->val;
2300 int msize = isize + vsize;
2301
2302 UnmapViewOfFile(ptr);
2303 ptr = MapViewOfFile(mapfile, FILE_MAP_ALL_ACCESS, 0, 0, msize);
2304 if (ptr == NULL) {
2305 *err = E_ALLOC;
2306 gretl_matrix_free(m);
2307 m = NULL;
2308 } else {
2309 memcpy(m->val, ptr + isize, vsize);
2310 }
2311 }
2312
2313 if (ptr != NULL) {
2314 UnmapViewOfFile(ptr);
2315 }
2316 if (mapfile != NULL) {
2317 CloseHandle(mapfile);
2318 }
2319 if (diskfile != NULL) {
2320 CloseHandle(diskfile);
2321 }
2322
2323 if (finalize) {
2324 gretl_remove(diskname);
2325 }
2326 g_free(diskname);
2327
2328 return m;
2329 }
2330
shm_finalize_matrix(const char * fname)2331 int shm_finalize_matrix (const char *fname)
2332 {
2333 gchar *diskname;
2334 int ret;
2335
2336 diskname = g_strdup_printf("%s%s", gretl_dotdir(), fname);
2337 ret = gretl_remove(diskname);
2338 g_free(diskname);
2339
2340 return ret;
2341 }
2342
2343 #else /* not MS Windows: Linux, OS X, etc. */
2344
2345 #include <unistd.h>
2346 #include <sys/mman.h>
2347 #include <fcntl.h>
2348 #include <errno.h>
2349
canonical_memname(const char * s)2350 static gchar *canonical_memname (const char *s)
2351 {
2352 gchar *name;
2353
2354 /* Ensure we have a suitable name for a memory-mapped
2355 file: should be something like "/foo.shm", with no
2356 path component other than "/").
2357 */
2358 if (*s != '/') {
2359 name = g_strdup_printf("/%s", s);
2360 } else {
2361 name = g_strdup(s);
2362 }
2363 gretl_charsub(name + 1, '/', '_');
2364
2365 return name;
2366 }
2367
shm_write_matrix(const gretl_matrix * m,const char * fname)2368 int shm_write_matrix (const gretl_matrix *m,
2369 const char *fname)
2370 {
2371 int vsize = m->rows * m->cols * sizeof *m->val;
2372 int msize = IDLEN + 2 * sizeof(int) + vsize;
2373 void *ptr = NULL;
2374 gchar *memname;
2375 int fd, err = 0;
2376
2377 memname = canonical_memname(fname);
2378
2379 fd = shm_open(memname, O_RDWR | O_CREAT, 0666);
2380 if (fd == -1) {
2381 fprintf(stderr, "mwrite: shm_open failed: %s\n", strerror(errno));
2382 err = E_FOPEN;
2383 }
2384
2385 if (ftruncate(fd, msize) == -1) {
2386 fprintf(stderr, "mwrite: ftruncate failed: %s\n", strerror(errno));
2387 err = E_ALLOC;
2388 }
2389
2390 if (!err) {
2391 ptr = mmap(NULL, msize, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
2392 if (ptr == MAP_FAILED) {
2393 fprintf(stderr, "mwrite: mmap failed: %s\n", strerror(errno));
2394 err = E_ALLOC;
2395 }
2396 }
2397
2398 if (ptr != NULL) {
2399 void *pos = ptr;
2400
2401 memcpy(pos, MATRIX_ID, IDLEN);
2402 pos += IDLEN;
2403 memcpy(pos, &m->rows, sizeof m->rows);
2404 pos += sizeof m->rows;
2405 memcpy(pos, &m->cols, sizeof m->cols);
2406 pos += sizeof m->cols;
2407 memcpy(pos, m->val, vsize);
2408
2409 munmap(ptr, msize);
2410 }
2411
2412 if (fd != -1) {
2413 close(fd);
2414 }
2415 if (err) {
2416 shm_unlink(memname);
2417 }
2418 g_free(memname);
2419
2420 return err;
2421 }
2422
shm_read_matrix(const char * fname,int finalize,int * err)2423 gretl_matrix *shm_read_matrix (const char *fname,
2424 int finalize,
2425 int *err)
2426 {
2427 gretl_matrix *m = NULL;
2428 char buf[IDLEN] = {0};
2429 int isize, msize = 0;
2430 void *ptr = NULL;
2431 gchar *memname;
2432 int r, c, fd;
2433
2434 memname = canonical_memname(fname);
2435
2436 fd = shm_open(memname, O_RDONLY, 0666);
2437 if (fd == -1) {
2438 fprintf(stderr, "mread: shm_open failed: %s\n", strerror(errno));
2439 *err = E_FOPEN;
2440 }
2441
2442 if (!*err) {
2443 msize = isize = IDLEN + 2 * sizeof(int);
2444 ptr = mmap(NULL, isize, PROT_READ, MAP_SHARED, fd, 0);
2445 if (ptr == MAP_FAILED) {
2446 fprintf(stderr, "mread: mmap failed: %s\n", strerror(errno));
2447 *err = E_ALLOC;
2448 }
2449 }
2450
2451 if (ptr != NULL) {
2452 void *pos = ptr;
2453
2454 memcpy(buf, pos, IDLEN);
2455 if (memcmp(buf, MATRIX_ID, IDLEN) == 0) {
2456 pos += IDLEN;
2457 memcpy(&r, pos, sizeof r);
2458 pos += sizeof r;
2459 memcpy(&c, pos, sizeof c);
2460 m = gretl_matrix_alloc(r, c);
2461 if (m == NULL) {
2462 *err = E_ALLOC;
2463 }
2464 }
2465 }
2466
2467 if (m != NULL) {
2468 int vsize = r * c * sizeof *m->val;
2469
2470 /* FIXME: use mremap() on Linux? */
2471 munmap(ptr, isize);
2472 msize = isize + vsize;
2473 ptr = mmap(NULL, msize, PROT_READ, MAP_SHARED, fd, 0);
2474 if (ptr == MAP_FAILED) {
2475 *err = E_ALLOC;
2476 gretl_matrix_free(m);
2477 m = NULL;
2478 } else {
2479 memcpy(m->val, ptr + isize, vsize);
2480 }
2481 }
2482
2483 if (ptr != NULL) {
2484 munmap(ptr, msize);
2485 }
2486 if (fd != -1) {
2487 close(fd);
2488 }
2489 if (finalize) {
2490 shm_unlink(memname);
2491 }
2492
2493 g_free(memname);
2494
2495 return m;
2496 }
2497
shm_finalize_matrix(const char * fname)2498 int shm_finalize_matrix (const char *fname)
2499 {
2500 gchar *memname;
2501 int ret;
2502
2503 memname = canonical_memname(fname);
2504 ret = shm_unlink(memname);
2505 g_free(memname);
2506
2507 return ret == -1 ? E_DATA: 0;
2508 }
2509
2510 #endif /* shared memory variants */
2511