Actual source code: ivec.c

petsc-3.3-p7 2013-05-11
  2: /**********************************ivec.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: ***********************************ivec.c*************************************/

 17: #include <../src/ksp/pc/impls/tfs/tfs.h>

 19: /* sorting args ivec.c ivec.c ... */
 20: #define   SORT_OPT        6     
 21: #define   SORT_STACK        50000


 24: /* allocate an address and size stack for sorter(s) */
 25: static void *offset_stack[2*SORT_STACK];
 26: static PetscInt   size_stack[SORT_STACK];

 28: /***********************************ivec.c*************************************/
 29: PetscInt *PCTFS_ivec_copy( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 30: {
 31:   while (n--)  {*arg1++ = *arg2++;}
 32:   return(arg1);
 33: }

 35: /***********************************ivec.c*************************************/
 36: PetscErrorCode PCTFS_ivec_zero( PetscInt *arg1,  PetscInt n)
 37: {
 39:   while (n--)  {*arg1++ = 0;}
 40:   return(0);
 41: }

 43: /***********************************ivec.c*************************************/
 44: PetscErrorCode PCTFS_ivec_set( PetscInt *arg1,  PetscInt arg2,  PetscInt n)
 45: {
 47:   while (n--)  {*arg1++ = arg2;}
 48:   return(0);
 49: }

 51: /***********************************ivec.c*************************************/
 52: PetscErrorCode PCTFS_ivec_max( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 53: {
 55:   while (n--)  {*arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++;}
 56:   return(0);
 57: }

 59: /***********************************ivec.c*************************************/
 60: PetscErrorCode PCTFS_ivec_min( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 61: {
 63:   while (n--)  {*(arg1) = PetscMin(*arg1,*arg2); arg1++; arg2++;}
 64:   return(0);
 65: }

 67: /***********************************ivec.c*************************************/
 68: PetscErrorCode PCTFS_ivec_mult( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 69: {
 71:   while (n--)  {*arg1++ *= *arg2++;}
 72:   return(0);
 73: }

 75: /***********************************ivec.c*************************************/
 76: PetscErrorCode PCTFS_ivec_add( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 77: {
 79:   while (n--)  {*arg1++ += *arg2++;}
 80:   return(0);
 81: }

 83: /***********************************ivec.c*************************************/
 84: PetscErrorCode PCTFS_ivec_lxor( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 85: {
 87:   while (n--) {*arg1=((*arg1 || *arg2) && !(*arg1 && *arg2)); arg1++; arg2++;}
 88:   return(0);
 89: }

 91: /***********************************ivec.c*************************************/
 92: PetscErrorCode PCTFS_ivec_xor( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
 93: {
 95:   while (n--)  {*arg1++ ^= *arg2++;}
 96:   return(0);
 97: }

 99: /***********************************ivec.c*************************************/
100: PetscErrorCode PCTFS_ivec_or( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
101: {
103:   while (n--)  {*arg1++ |= *arg2++;}
104:   return(0);
105: }

107: /***********************************ivec.c*************************************/
108: PetscErrorCode PCTFS_ivec_lor( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
109: {
111:   while (n--)  {*arg1 = (*arg1 || *arg2); arg1++; arg2++;}
112:   return(0);
113: }

115: /***********************************ivec.c*************************************/
116: PetscErrorCode PCTFS_ivec_and( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
117: {
119:   while (n--)  {*arg1++ &= *arg2++;}
120:   return(0);
121: }

123: /***********************************ivec.c*************************************/
124: PetscErrorCode PCTFS_ivec_land( PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
125: {
127:   while (n--) {*arg1 = (*arg1 && *arg2); arg1++; arg2++;}
128:   return(0);
129: }

131: /***********************************ivec.c*************************************/
132: PetscErrorCode PCTFS_ivec_and3( PetscInt *arg1,  PetscInt *arg2,  PetscInt *arg3, PetscInt n)
133: {
135:   while (n--)  {*arg1++ = (*arg2++ & *arg3++);}
136:   return(0);
137: }

139: /***********************************ivec.c*************************************/
140: PetscInt PCTFS_ivec_sum( PetscInt *arg1,  PetscInt n)
141: {
142:    PetscInt tmp = 0;


145:   while (n--) {tmp += *arg1++;}
146:   return(tmp);
147: }

149: /***********************************ivec.c*************************************/
150: PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2,  PetscInt n,  PetscInt *arg3)
151: {
152:    PetscInt i, j, type;


156:   /* LATER: if we're really motivated we can sort and then unsort */
157:   for (i=0;i<n;)
158:     {
159:       /* clump 'em for now */
160:       j=i+1;
161:       type = arg3[i];
162:       while ((j<n)&&(arg3[j]==type))
163:         {j++;}
164: 
165:       /* how many together */
166:       j -= i;

168:       /* call appropriate ivec function */
169:       if (type == GL_MAX)
170:         {PCTFS_ivec_max(arg1,arg2,j);}
171:       else if (type == GL_MIN)
172:         {PCTFS_ivec_min(arg1,arg2,j);}
173:       else if (type == GL_MULT)
174:         {PCTFS_ivec_mult(arg1,arg2,j);}
175:       else if (type == GL_ADD)
176:         {PCTFS_ivec_add(arg1,arg2,j);}
177:       else if (type == GL_B_XOR)
178:         {PCTFS_ivec_xor(arg1,arg2,j);}
179:       else if (type == GL_B_OR)
180:         {PCTFS_ivec_or(arg1,arg2,j);}
181:       else if (type == GL_B_AND)
182:         {PCTFS_ivec_and(arg1,arg2,j);}
183:       else if (type == GL_L_XOR)
184:         {PCTFS_ivec_lxor(arg1,arg2,j);}
185:       else if (type == GL_L_OR)
186:         {PCTFS_ivec_lor(arg1,arg2,j);}
187:       else if (type == GL_L_AND)
188:         {PCTFS_ivec_land(arg1,arg2,j);}
189:       else
190:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_ivec_non_uniform()!");

192:       arg1+=j; arg2+=j; i+=j;
193:     }
194:   return(0);
195: }

197: /***********************************ivec.c*************************************/
198: vfp PCTFS_ivec_fct_addr( PetscInt type)
199: {
201:   if (type == NON_UNIFORM)
202:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_non_uniform);}
203:   else if (type == GL_MAX)
204:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_max);}
205:   else if (type == GL_MIN)
206:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_min);}
207:   else if (type == GL_MULT)
208:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_mult);}
209:   else if (type == GL_ADD)
210:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_add);}
211:   else if (type == GL_B_XOR)
212:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_xor);}
213:   else if (type == GL_B_OR)
214:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_or);}
215:   else if (type == GL_B_AND)
216:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_and);}
217:   else if (type == GL_L_XOR)
218:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lxor);}
219:   else if (type == GL_L_OR)
220:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lor);}
221:   else if (type == GL_L_AND)
222:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_land);}

224:   /* catch all ... not good if we get here */
225:   return(NULL);
226: }

228: /******************************************************************************/
229: PetscErrorCode PCTFS_ivec_sort( PetscInt *ar,  PetscInt size)
230: {
231:    PetscInt *pi, *pj, temp;
232:    PetscInt **top_a = (PetscInt **) offset_stack;
233:    PetscInt *top_s = size_stack, *bottom_s = size_stack;


236:   /* we're really interested in the offset of the last element */
237:   /* ==> length of the list is now size + 1                    */
238:   size--;

240:   /* do until we're done ... return when stack is exhausted */
241:   for (;;)
242:     {
243:       /* if list is large enough use quicksort partition exchange code */
244:       if (size > SORT_OPT)
245:         {
246:           /* start up pointer at element 1 and down at size     */
247:           pi = ar+1;
248:           pj = ar+size;

250:           /* find middle element in list and swap w/ element 1 */
251:           SWAP(*(ar+(size>>1)),*pi)

253:           /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
254:           /* note ==> pivot_value in index 0                   */
255:           if (*pi > *pj)
256:             {SWAP(*pi,*pj)}
257:           if (*ar > *pj)
258:             {SWAP(*ar,*pj)}
259:           else if (*pi > *ar)
260:             {SWAP(*(ar),*(ar+1))}

262:           /* partition about pivot_value ...                              */
263:           /* note lists of length 2 are not guaranteed to be sorted */
264:           for(;;)
265:             {
266:               /* walk up ... and down ... swap if equal to pivot! */
267:               do pi++; while (*pi<*ar);
268:               do pj--; while (*pj>*ar);

270:               /* if we've crossed we're done */
271:               if (pj<pi) break;

273:               /* else swap */
274:               SWAP(*pi,*pj)
275:             }

277:           /* place pivot_value in it's correct location */
278:           SWAP(*ar,*pj)

280:           /* test stack_size to see if we've exhausted our stack */
281:           if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort() :: STACK EXHAUSTED!!!");

283:           /* push right hand child iff length > 1 */
284:           if ((*top_s = size-((PetscInt) (pi-ar))))
285:             {
286:               *(top_a++) = pi;
287:               size -= *top_s+2;
288:               top_s++;
289:             }
290:           /* set up for next loop iff there is something to do */
291:           else if (size -= *top_s+2)
292:             {;}
293:           /* might as well pop - note NR_OPT >=2 ==> we're ok! */
294:           else
295:             {
296:               ar = *(--top_a);
297:               size = *(--top_s);
298:             }
299:         }

301:       /* else sort small list directly then pop another off stack */
302:       else
303:         {
304:           /* insertion sort for bottom */
305:           for (pj=ar+1;pj<=ar+size;pj++)
306:             {
307:               temp = *pj;
308:               for (pi=pj-1;pi>=ar;pi--)
309:                 {
310:                   if (*pi <= temp) break;
311:                   *(pi+1)=*pi;
312:                 }
313:               *(pi+1)=temp;
314:             }

316:           /* check to see if stack is exhausted ==> DONE */
317:           if (top_s==bottom_s)   return(0);
318: 
319:           /* else pop another list from the stack */
320:           ar = *(--top_a);
321:           size = *(--top_s);
322:         }
323:     }
324:   return(0);
325: }

327: /******************************************************************************/
328: PetscErrorCode PCTFS_ivec_sort_companion( PetscInt *ar,  PetscInt *ar2,  PetscInt size)
329: {
330:    PetscInt *pi, *pj, temp, temp2;
331:    PetscInt **top_a = (PetscInt **)offset_stack;
332:    PetscInt *top_s = size_stack, *bottom_s = size_stack;
333:    PetscInt *pi2, *pj2;
334:    PetscInt mid;

337:   /* we're really interested in the offset of the last element */
338:   /* ==> length of the list is now size + 1                    */
339:   size--;

341:   /* do until we're done ... return when stack is exhausted */
342:   for (;;)
343:     {
344:       /* if list is large enough use quicksort partition exchange code */
345:       if (size > SORT_OPT)
346:         {
347:           /* start up pointer at element 1 and down at size     */
348:           mid = size>>1;
349:           pi = ar+1;
350:           pj = ar+mid;
351:           pi2 = ar2+1;
352:           pj2 = ar2+mid;

354:           /* find middle element in list and swap w/ element 1 */
355:           SWAP(*pi,*pj)
356:           SWAP(*pi2,*pj2)

358:           /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
359:           /* note ==> pivot_value in index 0                   */
360:           pj = ar+size;
361:           pj2 = ar2+size;
362:           if (*pi > *pj)
363:             {SWAP(*pi,*pj) SWAP(*pi2,*pj2)}
364:           if (*ar > *pj)
365:             {SWAP(*ar,*pj) SWAP(*ar2,*pj2)}
366:           else if (*pi > *ar)
367:             {SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1))}

369:           /* partition about pivot_value ...                              */
370:           /* note lists of length 2 are not guaranteed to be sorted */
371:           for(;;)
372:             {
373:               /* walk up ... and down ... swap if equal to pivot! */
374:               do {pi++; pi2++;} while (*pi<*ar);
375:               do {pj--; pj2--;} while (*pj>*ar);

377:               /* if we've crossed we're done */
378:               if (pj<pi) break;

380:               /* else swap */
381:               SWAP(*pi,*pj)
382:               SWAP(*pi2,*pj2)
383:             }

385:           /* place pivot_value in it's correct location */
386:           SWAP(*ar,*pj)
387:           SWAP(*ar2,*pj2)

389:           /* test stack_size to see if we've exhausted our stack */
390:           if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion() :: STACK EXHAUSTED!!!");

392:           /* push right hand child iff length > 1 */
393:           if ((*top_s = size-((PetscInt) (pi-ar))))
394:             {
395:               *(top_a++) = pi;
396:               *(top_a++) = pi2;
397:               size -= *top_s+2;
398:               top_s++;
399:             }
400:           /* set up for next loop iff there is something to do */
401:           else if (size -= *top_s+2)
402:             {;}
403:           /* might as well pop - note NR_OPT >=2 ==> we're ok! */
404:           else
405:             {
406:               ar2 = *(--top_a);
407:               ar  = *(--top_a);
408:               size = *(--top_s);
409:             }
410:         }

412:       /* else sort small list directly then pop another off stack */
413:       else
414:         {
415:           /* insertion sort for bottom */
416:           for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++)
417:             {
418:               temp = *pj;
419:               temp2 = *pj2;
420:               for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--)
421:                 {
422:                   if (*pi <= temp) break;
423:                   *(pi+1)=*pi;
424:                   *(pi2+1)=*pi2;
425:                 }
426:               *(pi+1)=temp;
427:               *(pi2+1)=temp2;
428:             }

430:           /* check to see if stack is exhausted ==> DONE */
431:           if (top_s==bottom_s)   return(0);
432: 
433:           /* else pop another list from the stack */
434:           ar2 = *(--top_a);
435:           ar  = *(--top_a);
436:           size = *(--top_s);
437:         }
438:     }
439:   return(0);
440: }

442: /******************************************************************************/
443: PetscErrorCode PCTFS_ivec_sort_companion_hack( PetscInt *ar,  PetscInt **ar2, PetscInt size)
444: {
445:    PetscInt *pi, *pj, temp, *ptr;
446:    PetscInt **top_a = (PetscInt **)offset_stack;
447:    PetscInt *top_s = size_stack, *bottom_s = size_stack;
448:    PetscInt **pi2, **pj2;
449:    PetscInt mid;

452:   /* we're really interested in the offset of the last element */
453:   /* ==> length of the list is now size + 1                    */
454:   size--;

456:   /* do until we're done ... return when stack is exhausted */
457:   for (;;)
458:     {
459:       /* if list is large enough use quicksort partition exchange code */
460:       if (size > SORT_OPT)
461:         {
462:           /* start up pointer at element 1 and down at size     */
463:           mid = size>>1;
464:           pi = ar+1;
465:           pj = ar+mid;
466:           pi2 = ar2+1;
467:           pj2 = ar2+mid;

469:           /* find middle element in list and swap w/ element 1 */
470:           SWAP(*pi,*pj)
471:           P_SWAP(*pi2,*pj2)

473:           /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
474:           /* note ==> pivot_value in index 0                   */
475:           pj = ar+size;
476:           pj2 = ar2+size;
477:           if (*pi > *pj)
478:             {SWAP(*pi,*pj) P_SWAP(*pi2,*pj2)}
479:           if (*ar > *pj)
480:             {SWAP(*ar,*pj) P_SWAP(*ar2,*pj2)}
481:           else if (*pi > *ar)
482:             {SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1))}

484:           /* partition about pivot_value ...                              */
485:           /* note lists of length 2 are not guaranteed to be sorted */
486:           for(;;)
487:             {
488:               /* walk up ... and down ... swap if equal to pivot! */
489:               do {pi++; pi2++;} while (*pi<*ar);
490:               do {pj--; pj2--;} while (*pj>*ar);

492:               /* if we've crossed we're done */
493:               if (pj<pi) break;

495:               /* else swap */
496:               SWAP(*pi,*pj)
497:               P_SWAP(*pi2,*pj2)
498:             }

500:           /* place pivot_value in it's correct location */
501:           SWAP(*ar,*pj)
502:           P_SWAP(*ar2,*pj2)

504:           /* test stack_size to see if we've exhausted our stack */
505:           if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion_hack() :: STACK EXHAUSTED!!!");

507:           /* push right hand child iff length > 1 */
508:           if ((*top_s = size-((PetscInt) (pi-ar))))
509:             {
510:               *(top_a++) = pi;
511:               *(top_a++) = (PetscInt*) pi2;
512:               size -= *top_s+2;
513:               top_s++;
514:             }
515:           /* set up for next loop iff there is something to do */
516:           else if (size -= *top_s+2)
517:             {;}
518:           /* might as well pop - note NR_OPT >=2 ==> we're ok! */
519:           else
520:             {
521:               ar2 = (PetscInt **) *(--top_a);
522:               ar  = *(--top_a);
523:               size = *(--top_s);
524:             }
525:         }

527:       /* else sort small list directly then pop another off stack */
528:       else
529:         {
530:           /* insertion sort for bottom */
531:           for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++)
532:             {
533:               temp = *pj;
534:               ptr = *pj2;
535:               for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--)
536:                 {
537:                   if (*pi <= temp) break;
538:                   *(pi+1)=*pi;
539:                   *(pi2+1)=*pi2;
540:                 }
541:               *(pi+1)=temp;
542:               *(pi2+1)=ptr;
543:             }

545:           /* check to see if stack is exhausted ==> DONE */
546:           if (top_s==bottom_s)   return(0);
547: 
548:           /* else pop another list from the stack */
549:           ar2 = (PetscInt **)*(--top_a);
550:           ar  = *(--top_a);
551:           size = *(--top_s);
552:         }
553:     }
554:   return(0);
555: }

557: /******************************************************************************/
558: PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type)
559: {
561:   if (type == SORT_INTEGER) {
562:     if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size);
563:     else PCTFS_ivec_sort((PetscInt*)ar1,size);
564:   } else if (type == SORT_INT_PTR) {
565:     if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt **)ar2,size);
566:     else PCTFS_ivec_sort((PetscInt*)ar1,size);
567:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!");
568:   return(0);
569: }

571: /***********************************ivec.c*************************************/
572: PetscInt PCTFS_ivec_linear_search( PetscInt item,  PetscInt *list,  PetscInt n)
573: {
574:   PetscInt tmp = n-1;
576:   while (n--)  {if (*list++ == item) {return(tmp-n);}}
577:   return(-1);
578: }

580: /***********************************ivec.c*************************************/
581: PetscInt PCTFS_ivec_binary_search( PetscInt item,  PetscInt *list,  PetscInt rh)
582: {
583:    PetscInt mid, lh=0;

585:   rh--;
586:   while (lh<=rh)
587:     {
588:       mid = (lh+rh)>>1;
589:       if (*(list+mid) == item)
590:         {return(mid);}
591:       if (*(list+mid) > item)
592:         {rh = mid-1;}
593:       else
594:         {lh = mid+1;}
595:     }
596:   return(-1);
597: }

599: /*********************************ivec.c*************************************/
600: PetscErrorCode PCTFS_rvec_copy( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
601: {
603:   while (n--)  {*arg1++ = *arg2++;}
604:   return(0);
605: }

607: /*********************************ivec.c*************************************/
608: PetscErrorCode PCTFS_rvec_zero( PetscScalar *arg1,  PetscInt n)
609: {
611:   while (n--)  {*arg1++ = 0.0;}
612:   return(0);
613: }

615: /***********************************ivec.c*************************************/
616: PetscErrorCode PCTFS_rvec_one( PetscScalar *arg1,  PetscInt n)
617: {
619:   while (n--)  {*arg1++ = 1.0;}
620:   return(0);
621: }

623: /***********************************ivec.c*************************************/
624: PetscErrorCode PCTFS_rvec_set( PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
625: {
627:   while (n--)  {*arg1++ = arg2;}
628:   return(0);
629: }

631: /***********************************ivec.c*************************************/
632: PetscErrorCode PCTFS_rvec_scale( PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
633: {
635:   while (n--)  {*arg1++ *= arg2;}
636:   return(0);
637: }

639: /*********************************ivec.c*************************************/
640: PetscErrorCode PCTFS_rvec_add( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
641: {
643:   while (n--)  {*arg1++ += *arg2++;}
644:   return(0);
645: }

647: /*********************************ivec.c*************************************/
648: PetscErrorCode PCTFS_rvec_mult( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
649: {
651:   while (n--)  {*arg1++ *= *arg2++;}
652:   return(0);
653: }

655: /*********************************ivec.c*************************************/
656: PetscErrorCode PCTFS_rvec_max( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
657: {
659:   while (n--)  {*arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++;}
660:   return(0);
661: }

663: /*********************************ivec.c*************************************/
664: PetscErrorCode PCTFS_rvec_max_abs( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
665: {
667:   while (n--)  {*arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++;}
668:   return(0);
669: }

671: /*********************************ivec.c*************************************/
672: PetscErrorCode PCTFS_rvec_min( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
673: {
675:   while (n--)  {*arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++;}
676:   return(0);
677: }

679: /*********************************ivec.c*************************************/
680: PetscErrorCode PCTFS_rvec_min_abs( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
681: {
683:   while (n--)  {*arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++;}
684:   return(0);
685: }

687: /*********************************ivec.c*************************************/
688: PetscErrorCode PCTFS_rvec_exists( PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
689: {
691:   while (n--)  {*arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++;}
692:   return(0);
693: }

695: /***********************************ivec.c*************************************/
696: PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2,  PetscInt n,  PetscInt *arg3)
697: {
698:    PetscInt i, j, type;

701:   /* LATER: if we're really motivated we can sort and then unsort */
702:   for (i=0;i<n;)
703:     {
704:       /* clump 'em for now */
705:       j=i+1;
706:       type = arg3[i];
707:       while ((j<n)&&(arg3[j]==type))
708:         {j++;}
709: 
710:       /* how many together */
711:       j -= i;

713:       /* call appropriate ivec function */
714:       if (type == GL_MAX)
715:         {PCTFS_rvec_max(arg1,arg2,j);}
716:       else if (type == GL_MIN)
717:         {PCTFS_rvec_min(arg1,arg2,j);}
718:       else if (type == GL_MULT)
719:         {PCTFS_rvec_mult(arg1,arg2,j);}
720:       else if (type == GL_ADD)
721:         {PCTFS_rvec_add(arg1,arg2,j);}
722:       else if (type == GL_MAX_ABS)
723:         {PCTFS_rvec_max_abs(arg1,arg2,j);}
724:       else if (type == GL_MIN_ABS)
725:         {PCTFS_rvec_min_abs(arg1,arg2,j);}
726:       else if (type == GL_EXISTS)
727:         {PCTFS_rvec_exists(arg1,arg2,j);}
728:       else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_rvec_non_uniform()!");

730:       arg1+=j; arg2+=j; i+=j;
731:     }
732:   return(0);
733: }

735: /***********************************ivec.c*************************************/
736: vfp PCTFS_rvec_fct_addr( PetscInt type)
737: {
738:   if (type == NON_UNIFORM)
739:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_non_uniform);}
740:   else if (type == GL_MAX)
741:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max);}
742:   else if (type == GL_MIN)
743:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min);}
744:   else if (type == GL_MULT)
745:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_mult);}
746:   else if (type == GL_ADD)
747:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_add);}
748:   else if (type == GL_MAX_ABS)
749:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max_abs);}
750:   else if (type == GL_MIN_ABS)
751:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min_abs);}
752:   else if (type == GL_EXISTS)
753:     {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_exists);}

755:   /* catch all ... not good if we get here */
756:   return(NULL);
757: }