Actual source code: gs.c
petsc-3.9.4 2018-09-11
2: /***********************************gs.c***************************************
4: Author: Henry M. Tufo III
6: e-mail: hmt@cs.brown.edu
8: snail-mail:
9: Division of Applied Mathematics
10: Brown University
11: Providence, RI 02912
13: Last Modification:
14: 6.21.97
15: ************************************gs.c**************************************/
17: /***********************************gs.c***************************************
18: File Description:
19: -----------------
21: ************************************gs.c**************************************/
23: #include <../src/ksp/pc/impls/tfs/tfs.h>
25: /* default length of number of items via tree - doubles if exceeded */
26: #define TREE_BUF_SZ 2048;
27: #define GS_VEC_SZ 1
31: /***********************************gs.c***************************************
32: Type: struct gather_scatter_id
33: ------------------------------
35: ************************************gs.c**************************************/
36: typedef struct gather_scatter_id {
37: PetscInt id;
38: PetscInt nel_min;
39: PetscInt nel_max;
40: PetscInt nel_sum;
41: PetscInt negl;
42: PetscInt gl_max;
43: PetscInt gl_min;
44: PetscInt repeats;
45: PetscInt ordered;
46: PetscInt positive;
47: PetscScalar *vals;
49: /* bit mask info */
50: PetscInt *my_proc_mask;
51: PetscInt mask_sz;
52: PetscInt *ngh_buf;
53: PetscInt ngh_buf_sz;
54: PetscInt *nghs;
55: PetscInt num_nghs;
56: PetscInt max_nghs;
57: PetscInt *pw_nghs;
58: PetscInt num_pw_nghs;
59: PetscInt *tree_nghs;
60: PetscInt num_tree_nghs;
62: PetscInt num_loads;
64: /* repeats == true -> local info */
65: PetscInt nel; /* number of unique elememts */
66: PetscInt *elms; /* of size nel */
67: PetscInt nel_total;
68: PetscInt *local_elms; /* of size nel_total */
69: PetscInt *companion; /* of size nel_total */
71: /* local info */
72: PetscInt num_local_total;
73: PetscInt local_strength;
74: PetscInt num_local;
75: PetscInt *num_local_reduce;
76: PetscInt **local_reduce;
77: PetscInt num_local_gop;
78: PetscInt *num_gop_local_reduce;
79: PetscInt **gop_local_reduce;
81: /* pairwise info */
82: PetscInt level;
83: PetscInt num_pairs;
84: PetscInt max_pairs;
85: PetscInt loc_node_pairs;
86: PetscInt max_node_pairs;
87: PetscInt min_node_pairs;
88: PetscInt avg_node_pairs;
89: PetscInt *pair_list;
90: PetscInt *msg_sizes;
91: PetscInt **node_list;
92: PetscInt len_pw_list;
93: PetscInt *pw_elm_list;
94: PetscScalar *pw_vals;
96: MPI_Request *msg_ids_in;
97: MPI_Request *msg_ids_out;
99: PetscScalar *out;
100: PetscScalar *in;
101: PetscInt msg_total;
103: /* tree - crystal accumulator info */
104: PetscInt max_left_over;
105: PetscInt *pre;
106: PetscInt *in_num;
107: PetscInt *out_num;
108: PetscInt **in_list;
109: PetscInt **out_list;
111: /* new tree work*/
112: PetscInt tree_nel;
113: PetscInt *tree_elms;
114: PetscScalar *tree_buf;
115: PetscScalar *tree_work;
117: PetscInt tree_map_sz;
118: PetscInt *tree_map_in;
119: PetscInt *tree_map_out;
121: /* current memory status */
122: PetscInt gl_bss_min;
123: PetscInt gl_perm_min;
125: /* max segment size for PCTFS_gs_gop_vec() */
126: PetscInt vec_sz;
128: /* hack to make paul happy */
129: MPI_Comm PCTFS_gs_comm;
131: } PCTFS_gs_id;
133: static PCTFS_gs_id *gsi_check_args(PetscInt *elms, PetscInt nel, PetscInt level);
134: static PetscErrorCode gsi_via_bit_mask(PCTFS_gs_id *gs);
135: static PetscErrorCode get_ngh_buf(PCTFS_gs_id *gs);
136: static PetscErrorCode set_pairwise(PCTFS_gs_id *gs);
137: static PCTFS_gs_id *gsi_new(void);
138: static PetscErrorCode set_tree(PCTFS_gs_id *gs);
140: /* same for all but vector flavor */
141: static PetscErrorCode PCTFS_gs_gop_local_out(PCTFS_gs_id *gs, PetscScalar *vals);
142: /* vector flavor */
143: static PetscErrorCode PCTFS_gs_gop_vec_local_out(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step);
145: static PetscErrorCode PCTFS_gs_gop_vec_plus(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt step);
146: static PetscErrorCode PCTFS_gs_gop_vec_pairwise_plus(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt step);
147: static PetscErrorCode PCTFS_gs_gop_vec_local_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step);
148: static PetscErrorCode PCTFS_gs_gop_vec_local_in_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step);
149: static PetscErrorCode PCTFS_gs_gop_vec_tree_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step);
152: static PetscErrorCode PCTFS_gs_gop_local_plus(PCTFS_gs_id *gs, PetscScalar *vals);
153: static PetscErrorCode PCTFS_gs_gop_local_in_plus(PCTFS_gs_id *gs, PetscScalar *vals);
155: static PetscErrorCode PCTFS_gs_gop_plus_hc(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt dim);
156: static PetscErrorCode PCTFS_gs_gop_pairwise_plus_hc(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt dim);
157: static PetscErrorCode PCTFS_gs_gop_tree_plus_hc(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt dim);
159: /* global vars */
160: /* from comm.c module */
162: static PetscInt num_gs_ids = 0;
164: /* should make this dynamic ... later */
165: static PetscInt msg_buf =MAX_MSG_BUF;
166: static PetscInt vec_sz =GS_VEC_SZ;
167: static PetscInt *tree_buf =NULL;
168: static PetscInt tree_buf_sz=0;
169: static PetscInt ntree =0;
171: /***************************************************************************/
172: PetscErrorCode PCTFS_gs_init_vec_sz(PetscInt size)
173: {
175: vec_sz = size;
176: return(0);
177: }
179: /******************************************************************************/
180: PetscErrorCode PCTFS_gs_init_msg_buf_sz(PetscInt buf_size)
181: {
183: msg_buf = buf_size;
184: return(0);
185: }
187: /******************************************************************************/
188: PCTFS_gs_id *PCTFS_gs_init(PetscInt *elms, PetscInt nel, PetscInt level)
189: {
190: PCTFS_gs_id *gs;
191: MPI_Group PCTFS_gs_group;
192: MPI_Comm PCTFS_gs_comm;
196: /* ensure that communication package has been initialized */
197: PCTFS_comm_init();
200: /* determines if we have enough dynamic/semi-static memory */
201: /* checks input, allocs and sets gd_id template */
202: gs = gsi_check_args(elms,nel,level);
204: /* only bit mask version up and working for the moment */
205: /* LATER :: get int list version working for sparse pblms */
206: gsi_via_bit_mask(gs);CHKERRABORT(PETSC_COMM_WORLD,ierr);
209: MPI_Comm_group(MPI_COMM_WORLD,&PCTFS_gs_group);CHKERRABORT(PETSC_COMM_WORLD,ierr);
210: MPI_Comm_create(MPI_COMM_WORLD,PCTFS_gs_group,&PCTFS_gs_comm);CHKERRABORT(PETSC_COMM_WORLD,ierr);
211: MPI_Group_free(&PCTFS_gs_group);CHKERRABORT(PETSC_COMM_WORLD,ierr);
213: gs->PCTFS_gs_comm=PCTFS_gs_comm;
215: return(gs);
216: }
218: /******************************************************************************/
219: static PCTFS_gs_id *gsi_new(void)
220: {
222: PCTFS_gs_id *gs;
223: gs = (PCTFS_gs_id*) malloc(sizeof(PCTFS_gs_id));
224: PetscMemzero(gs,sizeof(PCTFS_gs_id));CHKERRABORT(PETSC_COMM_WORLD,ierr);
225: return(gs);
226: }
228: /******************************************************************************/
229: static PCTFS_gs_id *gsi_check_args(PetscInt *in_elms, PetscInt nel, PetscInt level)
230: {
231: PetscInt i, j, k, t2;
232: PetscInt *companion, *elms, *unique, *iptr;
233: PetscInt num_local=0, *num_to_reduce, **local_reduce;
234: PetscInt oprs[] = {NON_UNIFORM,GL_MIN,GL_MAX,GL_ADD,GL_MIN,GL_MAX,GL_MIN,GL_B_AND};
235: PetscInt vals[sizeof(oprs)/sizeof(oprs[0])-1];
236: PetscInt work[sizeof(oprs)/sizeof(oprs[0])-1];
237: PCTFS_gs_id *gs;
241: if (!in_elms) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"elms point to nothing!!!\n");
242: if (nel<0) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"can't have fewer than 0 elms!!!\n");
244: if (nel==0) { PetscInfo(0,"I don't have any elements!!!\n");CHKERRABORT(PETSC_COMM_WORLD,ierr); }
246: /* get space for gs template */
247: gs = gsi_new();
248: gs->id = ++num_gs_ids;
250: /* hmt 6.4.99 */
251: /* caller can set global ids that don't participate to 0 */
252: /* PCTFS_gs_init ignores all zeros in elm list */
253: /* negative global ids are still invalid */
254: for (i=j=0; i<nel; i++) {
255: if (in_elms[i]!=0) j++;
256: }
258: k=nel; nel=j;
260: /* copy over in_elms list and create inverse map */
261: elms = (PetscInt*) malloc((nel+1)*sizeof(PetscInt));
262: companion = (PetscInt*) malloc(nel*sizeof(PetscInt));
264: for (i=j=0; i<k; i++) {
265: if (in_elms[i]!=0) { elms[j] = in_elms[i]; companion[j++] = i; }
266: }
268: if (j!=nel) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"nel j mismatch!\n");
270: /* pre-pass ... check to see if sorted */
271: elms[nel] = INT_MAX;
272: iptr = elms;
273: unique = elms+1;
274: j =0;
275: while (*iptr!=INT_MAX) {
276: if (*iptr++>*unique++) { j=1; break; }
277: }
279: /* set up inverse map */
280: if (j) {
281: PetscInfo(0,"gsi_check_args() :: elm list *not* sorted!\n");CHKERRABORT(PETSC_COMM_WORLD,ierr);
282: PCTFS_SMI_sort((void*)elms, (void*)companion, nel, SORT_INTEGER);CHKERRABORT(PETSC_COMM_WORLD,ierr);
283: } else { PetscInfo(0,"gsi_check_args() :: elm list sorted!\n");CHKERRABORT(PETSC_COMM_WORLD,ierr); }
284: elms[nel] = INT_MIN;
286: /* first pass */
287: /* determine number of unique elements, check pd */
288: for (i=k=0; i<nel; i+=j) {
289: t2 = elms[i];
290: j = ++i;
292: /* clump 'em for now */
293: while (elms[j]==t2) j++;
295: /* how many together and num local */
296: if (j-=i) { num_local++; k+=j; }
297: }
299: /* how many unique elements? */
300: gs->repeats = k;
301: gs->nel = nel-k;
304: /* number of repeats? */
305: gs->num_local = num_local;
306: num_local += 2;
307: gs->local_reduce = local_reduce=(PetscInt**)malloc(num_local*sizeof(PetscInt*));
308: gs->num_local_reduce = num_to_reduce=(PetscInt*) malloc(num_local*sizeof(PetscInt));
310: unique = (PetscInt*) malloc((gs->nel+1)*sizeof(PetscInt));
311: gs->elms = unique;
312: gs->nel_total = nel;
313: gs->local_elms = elms;
314: gs->companion = companion;
316: /* compess map as well as keep track of local ops */
317: for (num_local=i=j=0; i<gs->nel; i++) {
318: k = j;
319: t2 = unique[i] = elms[j];
320: companion[i] = companion[j];
322: while (elms[j]==t2) j++;
324: if ((t2=(j-k))>1) {
325: /* number together */
326: num_to_reduce[num_local] = t2++;
328: iptr = local_reduce[num_local++] = (PetscInt*)malloc(t2*sizeof(PetscInt));
330: /* to use binary searching don't remap until we check intersection */
331: *iptr++ = i;
333: /* note that we're skipping the first one */
334: while (++k<j) *(iptr++) = companion[k];
335: *iptr = -1;
336: }
337: }
339: /* sentinel for ngh_buf */
340: unique[gs->nel]=INT_MAX;
342: /* for two partition sort hack */
343: num_to_reduce[num_local] = 0;
344: local_reduce[num_local] = NULL;
345: num_to_reduce[++num_local] = 0;
346: local_reduce[num_local] = NULL;
348: /* load 'em up */
349: /* note one extra to hold NON_UNIFORM flag!!! */
350: vals[2] = vals[1] = vals[0] = nel;
351: if (gs->nel>0) {
352: vals[3] = unique[0];
353: vals[4] = unique[gs->nel-1];
354: } else {
355: vals[3] = INT_MAX;
356: vals[4] = INT_MIN;
357: }
358: vals[5] = level;
359: vals[6] = num_gs_ids;
361: /* GLOBAL: send 'em out */
362: PCTFS_giop(vals,work,sizeof(oprs)/sizeof(oprs[0])-1,oprs);CHKERRABORT(PETSC_COMM_WORLD,ierr);
364: /* must be semi-pos def - only pairwise depends on this */
365: /* LATER - remove this restriction */
366: if (vals[3]<0) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"gsi_check_args() :: system not semi-pos def \n");
367: if (vals[4]==INT_MAX) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"gsi_check_args() :: system ub too large !\n");
369: gs->nel_min = vals[0];
370: gs->nel_max = vals[1];
371: gs->nel_sum = vals[2];
372: gs->gl_min = vals[3];
373: gs->gl_max = vals[4];
374: gs->negl = vals[4]-vals[3]+1;
376: if (gs->negl<=0) SETERRABORT(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"gsi_check_args() :: system empty or neg :: %d\n");
378: /* LATER :: add level == -1 -> program selects level */
379: if (vals[5]<0) vals[5]=0;
380: else if (vals[5]>PCTFS_num_nodes) vals[5]=PCTFS_num_nodes;
381: gs->level = vals[5];
383: return(gs);
384: }
386: /******************************************************************************/
387: static PetscErrorCode gsi_via_bit_mask(PCTFS_gs_id *gs)
388: {
389: PetscInt i, nel, *elms;
390: PetscInt t1;
391: PetscInt **reduce;
392: PetscInt *map;
396: /* totally local removes ... PCTFS_ct_bits == 0 */
397: get_ngh_buf(gs);
399: if (gs->level) set_pairwise(gs);
400: if (gs->max_left_over) set_tree(gs);
402: /* intersection local and pairwise/tree? */
403: gs->num_local_total = gs->num_local;
404: gs->gop_local_reduce = gs->local_reduce;
405: gs->num_gop_local_reduce = gs->num_local_reduce;
407: map = gs->companion;
409: /* is there any local compression */
410: if (!gs->num_local) {
411: gs->local_strength = NONE;
412: gs->num_local_gop = 0;
413: } else {
414: /* ok find intersection */
415: map = gs->companion;
416: reduce = gs->local_reduce;
417: for (i=0, t1=0; i<gs->num_local; i++, reduce++) {
418: if ((PCTFS_ivec_binary_search(**reduce,gs->pw_elm_list,gs->len_pw_list)>=0) || PCTFS_ivec_binary_search(**reduce,gs->tree_map_in,gs->tree_map_sz)>=0) {
419: t1++;
420: if (gs->num_local_reduce[i]<=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nobody in list?");
421: gs->num_local_reduce[i] *= -1;
422: }
423: **reduce=map[**reduce];
424: }
426: /* intersection is empty */
427: if (!t1) {
428: gs->local_strength = FULL;
429: gs->num_local_gop = 0;
430: } else { /* intersection not empty */
431: gs->local_strength = PARTIAL;
433: PCTFS_SMI_sort((void*)gs->num_local_reduce, (void*)gs->local_reduce, gs->num_local + 1, SORT_INT_PTR);
435: gs->num_local_gop = t1;
436: gs->num_local_total = gs->num_local;
437: gs->num_local -= t1;
438: gs->gop_local_reduce = gs->local_reduce;
439: gs->num_gop_local_reduce = gs->num_local_reduce;
441: for (i=0; i<t1; i++) {
442: if (gs->num_gop_local_reduce[i]>=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"they aren't negative?");
443: gs->num_gop_local_reduce[i] *= -1;
444: gs->local_reduce++;
445: gs->num_local_reduce++;
446: }
447: gs->local_reduce++;
448: gs->num_local_reduce++;
449: }
450: }
452: elms = gs->pw_elm_list;
453: nel = gs->len_pw_list;
454: for (i=0; i<nel; i++) elms[i] = map[elms[i]];
456: elms = gs->tree_map_in;
457: nel = gs->tree_map_sz;
458: for (i=0; i<nel; i++) elms[i] = map[elms[i]];
460: /* clean up */
461: free((void*) gs->local_elms);
462: free((void*) gs->companion);
463: free((void*) gs->elms);
464: free((void*) gs->ngh_buf);
465: gs->local_elms = gs->companion = gs->elms = gs->ngh_buf = NULL;
466: return(0);
467: }
469: /******************************************************************************/
470: static PetscErrorCode place_in_tree(PetscInt elm)
471: {
472: PetscInt *tp, n;
475: if (ntree==tree_buf_sz) {
476: if (tree_buf_sz) {
477: tp = tree_buf;
478: n = tree_buf_sz;
479: tree_buf_sz<<=1;
480: tree_buf = (PetscInt*)malloc(tree_buf_sz*sizeof(PetscInt));
481: PCTFS_ivec_copy(tree_buf,tp,n);
482: free(tp);
483: } else {
484: tree_buf_sz = TREE_BUF_SZ;
485: tree_buf = (PetscInt*)malloc(tree_buf_sz*sizeof(PetscInt));
486: }
487: }
489: tree_buf[ntree++] = elm;
490: return(0);
491: }
493: /******************************************************************************/
494: static PetscErrorCode get_ngh_buf(PCTFS_gs_id *gs)
495: {
496: PetscInt i, j, npw=0, ntree_map=0;
497: PetscInt p_mask_size, ngh_buf_size, buf_size;
498: PetscInt *p_mask, *sh_proc_mask, *pw_sh_proc_mask;
499: PetscInt *ngh_buf, *buf1, *buf2;
500: PetscInt offset, per_load, num_loads, or_ct, start, end;
501: PetscInt *ptr1, *ptr2, i_start, negl, nel, *elms;
502: PetscInt oper=GL_B_OR;
503: PetscInt *ptr3, *t_mask, level, ct1, ct2;
507: /* to make life easier */
508: nel = gs->nel;
509: elms = gs->elms;
510: level = gs->level;
512: /* det #bytes needed for processor bit masks and init w/mask cor. to PCTFS_my_id */
513: p_mask = (PetscInt*) malloc(p_mask_size=PCTFS_len_bit_mask(PCTFS_num_nodes));
514: PCTFS_set_bit_mask(p_mask,p_mask_size,PCTFS_my_id);
516: /* allocate space for masks and info bufs */
517: gs->nghs = sh_proc_mask = (PetscInt*) malloc(p_mask_size);
518: gs->pw_nghs = pw_sh_proc_mask = (PetscInt*) malloc(p_mask_size);
519: gs->ngh_buf_sz = ngh_buf_size = p_mask_size*nel;
520: t_mask = (PetscInt*) malloc(p_mask_size);
521: gs->ngh_buf = ngh_buf = (PetscInt*) malloc(ngh_buf_size);
523: /* comm buffer size ... memory usage bounded by ~2*msg_buf */
524: /* had thought I could exploit rendezvous threshold */
526: /* default is one pass */
527: per_load = negl = gs->negl;
528: gs->num_loads = num_loads = 1;
529: i = p_mask_size*negl;
531: /* possible overflow on buffer size */
532: /* overflow hack */
533: if (i<0) i=INT_MAX;
535: buf_size = PetscMin(msg_buf,i);
537: /* can we do it? */
538: if (p_mask_size>buf_size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"get_ngh_buf() :: buf<pms :: %d>%d\n",p_mask_size,buf_size);
540: /* get PCTFS_giop buf space ... make *only* one malloc */
541: buf1 = (PetscInt*) malloc(buf_size<<1);
543: /* more than one gior exchange needed? */
544: if (buf_size!=i) {
545: per_load = buf_size/p_mask_size;
546: buf_size = per_load*p_mask_size;
547: gs->num_loads = num_loads = negl/per_load + (negl%per_load>0);
548: }
551: /* convert buf sizes from #bytes to #ints - 32 bit only! */
552: p_mask_size/=sizeof(PetscInt); ngh_buf_size/=sizeof(PetscInt); buf_size/=sizeof(PetscInt);
554: /* find PCTFS_giop work space */
555: buf2 = buf1+buf_size;
557: /* hold #ints needed for processor masks */
558: gs->mask_sz=p_mask_size;
560: /* init buffers */
561: PCTFS_ivec_zero(sh_proc_mask,p_mask_size);
562: PCTFS_ivec_zero(pw_sh_proc_mask,p_mask_size);
563: PCTFS_ivec_zero(ngh_buf,ngh_buf_size);
565: /* HACK reset tree info */
566: tree_buf = NULL;
567: tree_buf_sz = ntree = 0;
569: /* ok do it */
570: for (ptr1=ngh_buf,ptr2=elms,end=gs->gl_min,or_ct=i=0; or_ct<num_loads; or_ct++) {
571: /* identity for bitwise or is 000...000 */
572: PCTFS_ivec_zero(buf1,buf_size);
574: /* load msg buffer */
575: for (start=end,end+=per_load,i_start=i; (offset=*ptr2)<end; i++, ptr2++) {
576: offset = (offset-start)*p_mask_size;
577: PCTFS_ivec_copy(buf1+offset,p_mask,p_mask_size);
578: }
580: /* GLOBAL: pass buffer */
581: PCTFS_giop(buf1,buf2,buf_size,&oper);
584: /* unload buffer into ngh_buf */
585: ptr2=(elms+i_start);
586: for (ptr3=buf1,j=start; j<end; ptr3+=p_mask_size,j++) {
587: /* I own it ... may have to pairwise it */
588: if (j==*ptr2) {
589: /* do i share it w/anyone? */
590: ct1 = PCTFS_ct_bits((char*)ptr3,p_mask_size*sizeof(PetscInt));
591: /* guess not */
592: if (ct1<2) { ptr2++; ptr1+=p_mask_size; continue; }
594: /* i do ... so keep info and turn off my bit */
595: PCTFS_ivec_copy(ptr1,ptr3,p_mask_size);
596: PCTFS_ivec_xor(ptr1,p_mask,p_mask_size);
597: PCTFS_ivec_or(sh_proc_mask,ptr1,p_mask_size);
599: /* is it to be done pairwise? */
600: if (--ct1<=level) {
601: npw++;
603: /* turn on high bit to indicate pw need to process */
604: *ptr2++ |= TOP_BIT;
605: PCTFS_ivec_or(pw_sh_proc_mask,ptr1,p_mask_size);
606: ptr1 += p_mask_size;
607: continue;
608: }
610: /* get set for next and note that I have a tree contribution */
611: /* could save exact elm index for tree here -> save a search */
612: ptr2++; ptr1+=p_mask_size; ntree_map++;
613: } else { /* i don't but still might be involved in tree */
615: /* shared by how many? */
616: ct1 = PCTFS_ct_bits((char*)ptr3,p_mask_size*sizeof(PetscInt));
618: /* none! */
619: if (ct1<2) continue;
621: /* is it going to be done pairwise? but not by me of course!*/
622: if (--ct1<=level) continue;
623: }
624: /* LATER we're going to have to process it NOW */
625: /* nope ... tree it */
626: place_in_tree(j);
627: }
628: }
630: free((void*)t_mask);
631: free((void*)buf1);
633: gs->len_pw_list = npw;
634: gs->num_nghs = PCTFS_ct_bits((char*)sh_proc_mask,p_mask_size*sizeof(PetscInt));
636: /* expand from bit mask list to int list and save ngh list */
637: gs->nghs = (PetscInt*) malloc(gs->num_nghs * sizeof(PetscInt));
638: PCTFS_bm_to_proc((char*)sh_proc_mask,p_mask_size*sizeof(PetscInt),gs->nghs);
640: gs->num_pw_nghs = PCTFS_ct_bits((char*)pw_sh_proc_mask,p_mask_size*sizeof(PetscInt));
642: oper = GL_MAX;
643: ct1 = gs->num_nghs;
644: PCTFS_giop(&ct1,&ct2,1,&oper);
645: gs->max_nghs = ct1;
647: gs->tree_map_sz = ntree_map;
648: gs->max_left_over=ntree;
650: free((void*)p_mask);
651: free((void*)sh_proc_mask);
652: return(0);
653: }
655: /******************************************************************************/
656: static PetscErrorCode set_pairwise(PCTFS_gs_id *gs)
657: {
658: PetscInt i, j;
659: PetscInt p_mask_size;
660: PetscInt *p_mask, *sh_proc_mask, *tmp_proc_mask;
661: PetscInt *ngh_buf, *buf2;
662: PetscInt offset;
663: PetscInt *msg_list, *msg_size, **msg_nodes, nprs;
664: PetscInt *pairwise_elm_list, len_pair_list=0;
665: PetscInt *iptr, t1, i_start, nel, *elms;
666: PetscInt ct;
670: /* to make life easier */
671: nel = gs->nel;
672: elms = gs->elms;
673: ngh_buf = gs->ngh_buf;
674: sh_proc_mask = gs->pw_nghs;
676: /* need a few temp masks */
677: p_mask_size = PCTFS_len_bit_mask(PCTFS_num_nodes);
678: p_mask = (PetscInt*) malloc(p_mask_size);
679: tmp_proc_mask = (PetscInt*) malloc(p_mask_size);
681: /* set mask to my PCTFS_my_id's bit mask */
682: PCTFS_set_bit_mask(p_mask,p_mask_size,PCTFS_my_id);
684: p_mask_size /= sizeof(PetscInt);
686: len_pair_list = gs->len_pw_list;
687: gs->pw_elm_list = pairwise_elm_list=(PetscInt*)malloc((len_pair_list+1)*sizeof(PetscInt));
689: /* how many processors (nghs) do we have to exchange with? */
690: nprs = gs->num_pairs = PCTFS_ct_bits((char*)sh_proc_mask,p_mask_size*sizeof(PetscInt));
693: /* allocate space for PCTFS_gs_gop() info */
694: gs->pair_list = msg_list = (PetscInt*) malloc(sizeof(PetscInt)*nprs);
695: gs->msg_sizes = msg_size = (PetscInt*) malloc(sizeof(PetscInt)*nprs);
696: gs->node_list = msg_nodes = (PetscInt**) malloc(sizeof(PetscInt*)*(nprs+1));
698: /* init msg_size list */
699: PCTFS_ivec_zero(msg_size,nprs);
701: /* expand from bit mask list to int list */
702: PCTFS_bm_to_proc((char*)sh_proc_mask,p_mask_size*sizeof(PetscInt),msg_list);
704: /* keep list of elements being handled pairwise */
705: for (i=j=0; i<nel; i++) {
706: if (elms[i] & TOP_BIT) { elms[i] ^= TOP_BIT; pairwise_elm_list[j++] = i; }
707: }
708: pairwise_elm_list[j] = -1;
710: gs->msg_ids_out = (MPI_Request*) malloc(sizeof(MPI_Request)*(nprs+1));
711: gs->msg_ids_out[nprs] = MPI_REQUEST_NULL;
712: gs->msg_ids_in = (MPI_Request*) malloc(sizeof(MPI_Request)*(nprs+1));
713: gs->msg_ids_in[nprs] = MPI_REQUEST_NULL;
714: gs->pw_vals = (PetscScalar*) malloc(sizeof(PetscScalar)*len_pair_list*vec_sz);
716: /* find who goes to each processor */
717: for (i_start=i=0; i<nprs; i++) {
718: /* processor i's mask */
719: PCTFS_set_bit_mask(p_mask,p_mask_size*sizeof(PetscInt),msg_list[i]);
721: /* det # going to processor i */
722: for (ct=j=0; j<len_pair_list; j++) {
723: buf2 = ngh_buf+(pairwise_elm_list[j]*p_mask_size);
724: PCTFS_ivec_and3(tmp_proc_mask,p_mask,buf2,p_mask_size);
725: if (PCTFS_ct_bits((char*)tmp_proc_mask,p_mask_size*sizeof(PetscInt))) ct++;
726: }
727: msg_size[i] = ct;
728: i_start = PetscMax(i_start,ct);
730: /*space to hold nodes in message to first neighbor */
731: msg_nodes[i] = iptr = (PetscInt*) malloc(sizeof(PetscInt)*(ct+1));
733: for (j=0;j<len_pair_list;j++) {
734: buf2 = ngh_buf+(pairwise_elm_list[j]*p_mask_size);
735: PCTFS_ivec_and3(tmp_proc_mask,p_mask,buf2,p_mask_size);
736: if (PCTFS_ct_bits((char*)tmp_proc_mask,p_mask_size*sizeof(PetscInt))) *iptr++ = j;
737: }
738: *iptr = -1;
739: }
740: msg_nodes[nprs] = NULL;
742: j = gs->loc_node_pairs=i_start;
743: t1 = GL_MAX;
744: PCTFS_giop(&i_start,&offset,1,&t1);
745: gs->max_node_pairs = i_start;
747: i_start = j;
748: t1 = GL_MIN;
749: PCTFS_giop(&i_start,&offset,1,&t1);
750: gs->min_node_pairs = i_start;
752: i_start = j;
753: t1 = GL_ADD;
754: PCTFS_giop(&i_start,&offset,1,&t1);
755: gs->avg_node_pairs = i_start/PCTFS_num_nodes + 1;
757: i_start = nprs;
758: t1 = GL_MAX;
759: PCTFS_giop(&i_start,&offset,1,&t1);
760: gs->max_pairs = i_start;
763: /* remap pairwise in tail of gsi_via_bit_mask() */
764: gs->msg_total = PCTFS_ivec_sum(gs->msg_sizes,nprs);
765: gs->out = (PetscScalar*) malloc(sizeof(PetscScalar)*gs->msg_total*vec_sz);
766: gs->in = (PetscScalar*) malloc(sizeof(PetscScalar)*gs->msg_total*vec_sz);
768: /* reset malloc pool */
769: free((void*)p_mask);
770: free((void*)tmp_proc_mask);
771: return(0);
772: }
774: /* to do pruned tree just save ngh buf copy for each one and decode here!
775: ******************************************************************************/
776: static PetscErrorCode set_tree(PCTFS_gs_id *gs)
777: {
778: PetscInt i, j, n, nel;
779: PetscInt *iptr_in, *iptr_out, *tree_elms, *elms;
782: /* local work ptrs */
783: elms = gs->elms;
784: nel = gs->nel;
786: /* how many via tree */
787: gs->tree_nel = n = ntree;
788: gs->tree_elms = tree_elms = iptr_in = tree_buf;
789: gs->tree_buf = (PetscScalar*) malloc(sizeof(PetscScalar)*n*vec_sz);
790: gs->tree_work = (PetscScalar*) malloc(sizeof(PetscScalar)*n*vec_sz);
791: j = gs->tree_map_sz;
792: gs->tree_map_in = iptr_in = (PetscInt*) malloc(sizeof(PetscInt)*(j+1));
793: gs->tree_map_out = iptr_out = (PetscInt*) malloc(sizeof(PetscInt)*(j+1));
795: /* search the longer of the two lists */
796: /* note ... could save this info in get_ngh_buf and save searches */
797: if (n<=nel) {
798: /* bijective fct w/remap - search elm list */
799: for (i=0; i<n; i++) {
800: if ((j=PCTFS_ivec_binary_search(*tree_elms++,elms,nel))>=0) {*iptr_in++ = j; *iptr_out++ = i;}
801: }
802: } else {
803: for (i=0; i<nel; i++) {
804: if ((j=PCTFS_ivec_binary_search(*elms++,tree_elms,n))>=0) {*iptr_in++ = i; *iptr_out++ = j;}
805: }
806: }
808: /* sentinel */
809: *iptr_in = *iptr_out = -1;
810: return(0);
811: }
813: /******************************************************************************/
814: static PetscErrorCode PCTFS_gs_gop_local_out(PCTFS_gs_id *gs, PetscScalar *vals)
815: {
816: PetscInt *num, *map, **reduce;
817: PetscScalar tmp;
820: num = gs->num_gop_local_reduce;
821: reduce = gs->gop_local_reduce;
822: while ((map = *reduce++)) {
823: /* wall */
824: if (*num == 2) {
825: num++;
826: vals[map[1]] = vals[map[0]];
827: } else if (*num == 3) { /* corner shared by three elements */
828: num++;
829: vals[map[2]] = vals[map[1]] = vals[map[0]];
830: } else if (*num == 4) { /* corner shared by four elements */
831: num++;
832: vals[map[3]] = vals[map[2]] = vals[map[1]] = vals[map[0]];
833: } else { /* general case ... odd geoms ... 3D*/
834: num++;
835: tmp = *(vals + *map++);
836: while (*map >= 0) *(vals + *map++) = tmp;
837: }
838: }
839: return(0);
840: }
842: /******************************************************************************/
843: static PetscErrorCode PCTFS_gs_gop_local_plus(PCTFS_gs_id *gs, PetscScalar *vals)
844: {
845: PetscInt *num, *map, **reduce;
846: PetscScalar tmp;
849: num = gs->num_local_reduce;
850: reduce = gs->local_reduce;
851: while ((map = *reduce)) {
852: /* wall */
853: if (*num == 2) {
854: num++; reduce++;
855: vals[map[1]] = vals[map[0]] += vals[map[1]];
856: } else if (*num == 3) { /* corner shared by three elements */
857: num++; reduce++;
858: vals[map[2]]=vals[map[1]]=vals[map[0]]+=(vals[map[1]]+vals[map[2]]);
859: } else if (*num == 4) { /* corner shared by four elements */
860: num++; reduce++;
861: vals[map[1]]=vals[map[2]]=vals[map[3]]=vals[map[0]] += (vals[map[1]] + vals[map[2]] + vals[map[3]]);
862: } else { /* general case ... odd geoms ... 3D*/
863: num++;
864: tmp = 0.0;
865: while (*map >= 0) tmp += *(vals + *map++);
867: map = *reduce++;
868: while (*map >= 0) *(vals + *map++) = tmp;
869: }
870: }
871: return(0);
872: }
874: /******************************************************************************/
875: static PetscErrorCode PCTFS_gs_gop_local_in_plus(PCTFS_gs_id *gs, PetscScalar *vals)
876: {
877: PetscInt *num, *map, **reduce;
878: PetscScalar *base;
881: num = gs->num_gop_local_reduce;
882: reduce = gs->gop_local_reduce;
883: while ((map = *reduce++)) {
884: /* wall */
885: if (*num == 2) {
886: num++;
887: vals[map[0]] += vals[map[1]];
888: } else if (*num == 3) { /* corner shared by three elements */
889: num++;
890: vals[map[0]] += (vals[map[1]] + vals[map[2]]);
891: } else if (*num == 4) { /* corner shared by four elements */
892: num++;
893: vals[map[0]] += (vals[map[1]] + vals[map[2]] + vals[map[3]]);
894: } else { /* general case ... odd geoms ... 3D*/
895: num++;
896: base = vals + *map++;
897: while (*map >= 0) *base += *(vals + *map++);
898: }
899: }
900: return(0);
901: }
903: /******************************************************************************/
904: PetscErrorCode PCTFS_gs_free(PCTFS_gs_id *gs)
905: {
906: PetscInt i;
910: MPI_Comm_free(&gs->PCTFS_gs_comm);CHKERRABORT(PETSC_COMM_WORLD,ierr);
911: if (gs->nghs) free((void*) gs->nghs);
912: if (gs->pw_nghs) free((void*) gs->pw_nghs);
914: /* tree */
915: if (gs->max_left_over) {
916: if (gs->tree_elms) free((void*) gs->tree_elms);
917: if (gs->tree_buf) free((void*) gs->tree_buf);
918: if (gs->tree_work) free((void*) gs->tree_work);
919: if (gs->tree_map_in) free((void*) gs->tree_map_in);
920: if (gs->tree_map_out) free((void*) gs->tree_map_out);
921: }
923: /* pairwise info */
924: if (gs->num_pairs) {
925: /* should be NULL already */
926: if (gs->ngh_buf) free((void*) gs->ngh_buf);
927: if (gs->elms) free((void*) gs->elms);
928: if (gs->local_elms) free((void*) gs->local_elms);
929: if (gs->companion) free((void*) gs->companion);
931: /* only set if pairwise */
932: if (gs->vals) free((void*) gs->vals);
933: if (gs->in) free((void*) gs->in);
934: if (gs->out) free((void*) gs->out);
935: if (gs->msg_ids_in) free((void*) gs->msg_ids_in);
936: if (gs->msg_ids_out) free((void*) gs->msg_ids_out);
937: if (gs->pw_vals) free((void*) gs->pw_vals);
938: if (gs->pw_elm_list) free((void*) gs->pw_elm_list);
939: if (gs->node_list) {
940: for (i=0;i<gs->num_pairs;i++) {
941: if (gs->node_list[i]) {
942: free((void*) gs->node_list[i]);
943: }
944: }
945: free((void*) gs->node_list);
946: }
947: if (gs->msg_sizes) free((void*) gs->msg_sizes);
948: if (gs->pair_list) free((void*) gs->pair_list);
949: }
951: /* local info */
952: if (gs->num_local_total>=0) {
953: for (i=0;i<gs->num_local_total+1;i++) {
954: if (gs->num_gop_local_reduce[i]) free((void*) gs->gop_local_reduce[i]);
955: }
956: }
958: /* if intersection tree/pairwise and local isn't empty */
959: if (gs->gop_local_reduce) free((void*) gs->gop_local_reduce);
960: if (gs->num_gop_local_reduce) free((void*) gs->num_gop_local_reduce);
962: free((void*) gs);
963: return(0);
964: }
966: /******************************************************************************/
967: PetscErrorCode PCTFS_gs_gop_vec(PCTFS_gs_id *gs, PetscScalar *vals, const char *op, PetscInt step)
968: {
972: switch (*op) {
973: case '+':
974: PCTFS_gs_gop_vec_plus(gs,vals,step);
975: break;
976: default:
977: PetscInfo1(0,"PCTFS_gs_gop_vec() :: %c is not a valid op\n",op[0]);
978: PetscInfo(0,"PCTFS_gs_gop_vec() :: default :: plus\n");
979: PCTFS_gs_gop_vec_plus(gs,vals,step);
980: break;
981: }
982: return(0);
983: }
985: /******************************************************************************/
986: static PetscErrorCode PCTFS_gs_gop_vec_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step)
987: {
989: if (!gs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_gs_gop_vec() passed NULL gs handle!!!");
991: /* local only operations!!! */
992: if (gs->num_local) PCTFS_gs_gop_vec_local_plus(gs,vals,step);
994: /* if intersection tree/pairwise and local isn't empty */
995: if (gs->num_local_gop) {
996: PCTFS_gs_gop_vec_local_in_plus(gs,vals,step);
998: /* pairwise */
999: if (gs->num_pairs) PCTFS_gs_gop_vec_pairwise_plus(gs,vals,step);
1001: /* tree */
1002: else if (gs->max_left_over) PCTFS_gs_gop_vec_tree_plus(gs,vals,step);
1004: PCTFS_gs_gop_vec_local_out(gs,vals,step);
1005: } else { /* if intersection tree/pairwise and local is empty */
1006: /* pairwise */
1007: if (gs->num_pairs) PCTFS_gs_gop_vec_pairwise_plus(gs,vals,step);
1009: /* tree */
1010: else if (gs->max_left_over) PCTFS_gs_gop_vec_tree_plus(gs,vals,step);
1011: }
1012: return(0);
1013: }
1015: /******************************************************************************/
1016: static PetscErrorCode PCTFS_gs_gop_vec_local_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step)
1017: {
1018: PetscInt *num, *map, **reduce;
1019: PetscScalar *base;
1022: num = gs->num_local_reduce;
1023: reduce = gs->local_reduce;
1024: while ((map = *reduce)) {
1025: base = vals + map[0] * step;
1027: /* wall */
1028: if (*num == 2) {
1029: num++; reduce++;
1030: PCTFS_rvec_add (base,vals+map[1]*step,step);
1031: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1032: } else if (*num == 3) { /* corner shared by three elements */
1033: num++; reduce++;
1034: PCTFS_rvec_add (base,vals+map[1]*step,step);
1035: PCTFS_rvec_add (base,vals+map[2]*step,step);
1036: PCTFS_rvec_copy(vals+map[2]*step,base,step);
1037: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1038: } else if (*num == 4) { /* corner shared by four elements */
1039: num++; reduce++;
1040: PCTFS_rvec_add (base,vals+map[1]*step,step);
1041: PCTFS_rvec_add (base,vals+map[2]*step,step);
1042: PCTFS_rvec_add (base,vals+map[3]*step,step);
1043: PCTFS_rvec_copy(vals+map[3]*step,base,step);
1044: PCTFS_rvec_copy(vals+map[2]*step,base,step);
1045: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1046: } else { /* general case ... odd geoms ... 3D */
1047: num++;
1048: while (*++map >= 0) PCTFS_rvec_add (base,vals+*map*step,step);
1050: map = *reduce;
1051: while (*++map >= 0) PCTFS_rvec_copy(vals+*map*step,base,step);
1053: reduce++;
1054: }
1055: }
1056: return(0);
1057: }
1059: /******************************************************************************/
1060: static PetscErrorCode PCTFS_gs_gop_vec_local_in_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step)
1061: {
1062: PetscInt *num, *map, **reduce;
1063: PetscScalar *base;
1066: num = gs->num_gop_local_reduce;
1067: reduce = gs->gop_local_reduce;
1068: while ((map = *reduce++)) {
1069: base = vals + map[0] * step;
1071: /* wall */
1072: if (*num == 2) {
1073: num++;
1074: PCTFS_rvec_add(base,vals+map[1]*step,step);
1075: } else if (*num == 3) { /* corner shared by three elements */
1076: num++;
1077: PCTFS_rvec_add(base,vals+map[1]*step,step);
1078: PCTFS_rvec_add(base,vals+map[2]*step,step);
1079: } else if (*num == 4) { /* corner shared by four elements */
1080: num++;
1081: PCTFS_rvec_add(base,vals+map[1]*step,step);
1082: PCTFS_rvec_add(base,vals+map[2]*step,step);
1083: PCTFS_rvec_add(base,vals+map[3]*step,step);
1084: } else { /* general case ... odd geoms ... 3D*/
1085: num++;
1086: while (*++map >= 0) PCTFS_rvec_add(base,vals+*map*step,step);
1087: }
1088: }
1089: return(0);
1090: }
1092: /******************************************************************************/
1093: static PetscErrorCode PCTFS_gs_gop_vec_local_out(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step)
1094: {
1095: PetscInt *num, *map, **reduce;
1096: PetscScalar *base;
1099: num = gs->num_gop_local_reduce;
1100: reduce = gs->gop_local_reduce;
1101: while ((map = *reduce++)) {
1102: base = vals + map[0] * step;
1104: /* wall */
1105: if (*num == 2) {
1106: num++;
1107: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1108: } else if (*num == 3) { /* corner shared by three elements */
1109: num++;
1110: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1111: PCTFS_rvec_copy(vals+map[2]*step,base,step);
1112: } else if (*num == 4) { /* corner shared by four elements */
1113: num++;
1114: PCTFS_rvec_copy(vals+map[1]*step,base,step);
1115: PCTFS_rvec_copy(vals+map[2]*step,base,step);
1116: PCTFS_rvec_copy(vals+map[3]*step,base,step);
1117: } else { /* general case ... odd geoms ... 3D*/
1118: num++;
1119: while (*++map >= 0) PCTFS_rvec_copy(vals+*map*step,base,step);
1120: }
1121: }
1122: return(0);
1123: }
1125: /******************************************************************************/
1126: static PetscErrorCode PCTFS_gs_gop_vec_pairwise_plus(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt step)
1127: {
1128: PetscScalar *dptr1, *dptr2, *dptr3, *in1, *in2;
1129: PetscInt *iptr, *msg_list, *msg_size, **msg_nodes;
1130: PetscInt *pw, *list, *size, **nodes;
1131: MPI_Request *msg_ids_in, *msg_ids_out, *ids_in, *ids_out;
1132: MPI_Status status;
1133: PetscBLASInt i1 = 1,dstep;
1137: /* strip and load s */
1138: msg_list = list = gs->pair_list;
1139: msg_size = size = gs->msg_sizes;
1140: msg_nodes = nodes = gs->node_list;
1141: iptr = pw = gs->pw_elm_list;
1142: dptr1 = dptr3 = gs->pw_vals;
1143: msg_ids_in = ids_in = gs->msg_ids_in;
1144: msg_ids_out = ids_out = gs->msg_ids_out;
1145: dptr2 = gs->out;
1146: in1=in2 = gs->in;
1148: /* post the receives */
1149: /* msg_nodes=nodes; */
1150: do {
1151: /* Should MPI_ANY_SOURCE be replaced by *list ? In that case do the
1152: second one *list and do list++ afterwards */
1153: MPI_Irecv(in1, *size *step, MPIU_SCALAR, MPI_ANY_SOURCE, MSGTAG1 + *list, gs->PCTFS_gs_comm, msg_ids_in);
1154: list++;msg_ids_in++;
1155: in1 += *size++ *step;
1156: } while (*++msg_nodes);
1157: msg_nodes=nodes;
1159: /* load gs values into in out gs buffers */
1160: while (*iptr >= 0) {
1161: PCTFS_rvec_copy(dptr3,in_vals + *iptr*step,step);
1162: dptr3+=step;
1163: iptr++;
1164: }
1166: /* load out buffers and post the sends */
1167: while ((iptr = *msg_nodes++)) {
1168: dptr3 = dptr2;
1169: while (*iptr >= 0) {
1170: PCTFS_rvec_copy(dptr2,dptr1 + *iptr*step,step);
1171: dptr2+=step;
1172: iptr++;
1173: }
1174: MPI_Isend(dptr3, *msg_size *step, MPIU_SCALAR, *msg_list, MSGTAG1+PCTFS_my_id, gs->PCTFS_gs_comm, msg_ids_out);
1175: msg_size++; msg_list++;msg_ids_out++;
1176: }
1178: /* tree */
1179: if (gs->max_left_over) PCTFS_gs_gop_vec_tree_plus(gs,in_vals,step);
1181: /* process the received data */
1182: msg_nodes=nodes;
1183: while ((iptr = *nodes++)) {
1184: PetscScalar d1 = 1.0;
1186: /* Should I check the return value of MPI_Wait() or status? */
1187: /* Can this loop be replaced by a call to MPI_Waitall()? */
1188: MPI_Wait(ids_in, &status);
1189: ids_in++;
1190: while (*iptr >= 0) {
1191: PetscBLASIntCast(step,&dstep);
1192: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&dstep,&d1,in2,&i1,dptr1 + *iptr*step,&i1));
1193: in2+=step;
1194: iptr++;
1195: }
1196: }
1198: /* replace vals */
1199: while (*pw >= 0) {
1200: PCTFS_rvec_copy(in_vals + *pw*step,dptr1,step);
1201: dptr1+=step;
1202: pw++;
1203: }
1205: /* clear isend message handles */
1206: /* This changed for clarity though it could be the same */
1208: /* Should I check the return value of MPI_Wait() or status? */
1209: /* Can this loop be replaced by a call to MPI_Waitall()? */
1210: while (*msg_nodes++) {
1211: MPI_Wait(ids_out, &status);
1212: ids_out++;
1213: }
1214: return(0);
1215: }
1217: /******************************************************************************/
1218: static PetscErrorCode PCTFS_gs_gop_vec_tree_plus(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt step)
1219: {
1220: PetscInt size, *in, *out;
1221: PetscScalar *buf, *work;
1222: PetscInt op[] = {GL_ADD,0};
1223: PetscBLASInt i1 = 1;
1225: PetscBLASInt dstep;
1228: /* copy over to local variables */
1229: in = gs->tree_map_in;
1230: out = gs->tree_map_out;
1231: buf = gs->tree_buf;
1232: work = gs->tree_work;
1233: size = gs->tree_nel*step;
1235: /* zero out collection buffer */
1236: PCTFS_rvec_zero(buf,size);
1239: /* copy over my contributions */
1240: while (*in >= 0) {
1241: PetscBLASIntCast(step,&dstep);
1242: PetscStackCallBLAS("BLAScopy",BLAScopy_(&dstep,vals + *in++ * step,&i1,buf + *out++ * step,&i1));
1243: }
1245: /* perform fan in/out on full buffer */
1246: /* must change PCTFS_grop to handle the blas */
1247: PCTFS_grop(buf,work,size,op);
1249: /* reset */
1250: in = gs->tree_map_in;
1251: out = gs->tree_map_out;
1253: /* get the portion of the results I need */
1254: while (*in >= 0) {
1255: PetscBLASIntCast(step,&dstep);
1256: PetscStackCallBLAS("BLAScopy",BLAScopy_(&dstep,buf + *out++ * step,&i1,vals + *in++ * step,&i1));
1257: }
1258: return(0);
1259: }
1261: /******************************************************************************/
1262: PetscErrorCode PCTFS_gs_gop_hc(PCTFS_gs_id *gs, PetscScalar *vals, const char *op, PetscInt dim)
1263: {
1267: switch (*op) {
1268: case '+':
1269: PCTFS_gs_gop_plus_hc(gs,vals,dim);
1270: break;
1271: default:
1272: PetscInfo1(0,"PCTFS_gs_gop_hc() :: %c is not a valid op\n",op[0]);
1273: PetscInfo(0,"PCTFS_gs_gop_hc() :: default :: plus\n");
1274: PCTFS_gs_gop_plus_hc(gs,vals,dim);
1275: break;
1276: }
1277: return(0);
1278: }
1280: /******************************************************************************/
1281: static PetscErrorCode PCTFS_gs_gop_plus_hc(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt dim)
1282: {
1284: /* if there's nothing to do return */
1285: if (dim<=0) return(0);
1287: /* can't do more dimensions then exist */
1288: dim = PetscMin(dim,PCTFS_i_log2_num_nodes);
1290: /* local only operations!!! */
1291: if (gs->num_local) PCTFS_gs_gop_local_plus(gs,vals);
1293: /* if intersection tree/pairwise and local isn't empty */
1294: if (gs->num_local_gop) {
1295: PCTFS_gs_gop_local_in_plus(gs,vals);
1297: /* pairwise will do tree inside ... */
1298: if (gs->num_pairs) PCTFS_gs_gop_pairwise_plus_hc(gs,vals,dim); /* tree only */
1299: else if (gs->max_left_over) PCTFS_gs_gop_tree_plus_hc(gs,vals,dim);
1301: PCTFS_gs_gop_local_out(gs,vals);
1302: } else { /* if intersection tree/pairwise and local is empty */
1303: /* pairwise will do tree inside */
1304: if (gs->num_pairs) PCTFS_gs_gop_pairwise_plus_hc(gs,vals,dim); /* tree */
1305: else if (gs->max_left_over) PCTFS_gs_gop_tree_plus_hc(gs,vals,dim);
1306: }
1307: return(0);
1308: }
1310: /******************************************************************************/
1311: static PetscErrorCode PCTFS_gs_gop_pairwise_plus_hc(PCTFS_gs_id *gs, PetscScalar *in_vals, PetscInt dim)
1312: {
1313: PetscScalar *dptr1, *dptr2, *dptr3, *in1, *in2;
1314: PetscInt *iptr, *msg_list, *msg_size, **msg_nodes;
1315: PetscInt *pw, *list, *size, **nodes;
1316: MPI_Request *msg_ids_in, *msg_ids_out, *ids_in, *ids_out;
1317: MPI_Status status;
1318: PetscInt i, mask=1;
1322: for (i=1; i<dim; i++) { mask<<=1; mask++; }
1324: /* strip and load s */
1325: msg_list = list = gs->pair_list;
1326: msg_size = size = gs->msg_sizes;
1327: msg_nodes = nodes = gs->node_list;
1328: iptr = pw = gs->pw_elm_list;
1329: dptr1 = dptr3 = gs->pw_vals;
1330: msg_ids_in = ids_in = gs->msg_ids_in;
1331: msg_ids_out = ids_out = gs->msg_ids_out;
1332: dptr2 = gs->out;
1333: in1 = in2 = gs->in;
1335: /* post the receives */
1336: /* msg_nodes=nodes; */
1337: do {
1338: /* Should MPI_ANY_SOURCE be replaced by *list ? In that case do the
1339: second one *list and do list++ afterwards */
1340: if ((PCTFS_my_id|mask)==(*list|mask)) {
1341: MPI_Irecv(in1, *size, MPIU_SCALAR, MPI_ANY_SOURCE, MSGTAG1 + *list, gs->PCTFS_gs_comm, msg_ids_in);
1342: list++; msg_ids_in++;in1 += *size++;
1343: } else { list++; size++; }
1344: } while (*++msg_nodes);
1346: /* load gs values into in out gs buffers */
1347: while (*iptr >= 0) *dptr3++ = *(in_vals + *iptr++);
1349: /* load out buffers and post the sends */
1350: msg_nodes=nodes;
1351: list = msg_list;
1352: while ((iptr = *msg_nodes++)) {
1353: if ((PCTFS_my_id|mask)==(*list|mask)) {
1354: dptr3 = dptr2;
1355: while (*iptr >= 0) *dptr2++ = *(dptr1 + *iptr++);
1356: /* CHECK PERSISTENT COMMS MODE FOR ALL THIS STUFF */
1357: /* is msg_ids_out++ correct? */
1358: MPI_Isend(dptr3, *msg_size, MPIU_SCALAR, *list, MSGTAG1+PCTFS_my_id, gs->PCTFS_gs_comm, msg_ids_out);
1359: msg_size++;list++;msg_ids_out++;
1360: } else {list++; msg_size++;}
1361: }
1363: /* do the tree while we're waiting */
1364: if (gs->max_left_over) PCTFS_gs_gop_tree_plus_hc(gs,in_vals,dim);
1366: /* process the received data */
1367: msg_nodes=nodes;
1368: list = msg_list;
1369: while ((iptr = *nodes++)) {
1370: if ((PCTFS_my_id|mask)==(*list|mask)) {
1371: /* Should I check the return value of MPI_Wait() or status? */
1372: /* Can this loop be replaced by a call to MPI_Waitall()? */
1373: MPI_Wait(ids_in, &status);
1374: ids_in++;
1375: while (*iptr >= 0) *(dptr1 + *iptr++) += *in2++;
1376: }
1377: list++;
1378: }
1380: /* replace vals */
1381: while (*pw >= 0) *(in_vals + *pw++) = *dptr1++;
1383: /* clear isend message handles */
1384: /* This changed for clarity though it could be the same */
1385: while (*msg_nodes++) {
1386: if ((PCTFS_my_id|mask)==(*msg_list|mask)) {
1387: /* Should I check the return value of MPI_Wait() or status? */
1388: /* Can this loop be replaced by a call to MPI_Waitall()? */
1389: MPI_Wait(ids_out, &status);
1390: ids_out++;
1391: }
1392: msg_list++;
1393: }
1394: return(0);
1395: }
1397: /******************************************************************************/
1398: static PetscErrorCode PCTFS_gs_gop_tree_plus_hc(PCTFS_gs_id *gs, PetscScalar *vals, PetscInt dim)
1399: {
1400: PetscInt size;
1401: PetscInt *in, *out;
1402: PetscScalar *buf, *work;
1403: PetscInt op[] = {GL_ADD,0};
1406: in = gs->tree_map_in;
1407: out = gs->tree_map_out;
1408: buf = gs->tree_buf;
1409: work = gs->tree_work;
1410: size = gs->tree_nel;
1412: PCTFS_rvec_zero(buf,size);
1414: while (*in >= 0) *(buf + *out++) = *(vals + *in++);
1416: in = gs->tree_map_in;
1417: out = gs->tree_map_out;
1419: PCTFS_grop_hc(buf,work,size,op,dim);
1421: while (*in >= 0) *(vals + *in++) = *(buf + *out++);
1422: return(0);
1423: }