Actual source code: gs.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

  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);
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: }