Actual source code: ivec.c

      3: /**********************************ivec.c**************************************

      5: Author: Henry M. Tufo III

      7: e-mail: hmt@cs.brown.edu

      9: snail-mail:
     10: Division of Applied Mathematics
     11: Brown University
     12: Providence, RI 02912

     14: Last Modification:
     15: 6.21.97
     16: ***********************************ivec.c*************************************/

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

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

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

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

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

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

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

     60: /***********************************ivec.c*************************************/
     61: PetscErrorCode PCTFS_ivec_min(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
     62: {
     64:   while (n--) {
     65:     *(arg1) = PetscMin(*arg1,*arg2);
     66:     arg1++;
     67:     arg2++;
     68:   }
     69:   return(0);
     70: }

     72: /***********************************ivec.c*************************************/
     73: PetscErrorCode PCTFS_ivec_mult(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
     74: {
     76:   while (n--) *arg1++ *= *arg2++;
     77:   return(0);
     78: }

     80: /***********************************ivec.c*************************************/
     81: PetscErrorCode PCTFS_ivec_add(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
     82: {
     84:   while (n--) *arg1++ += *arg2++;
     85:   return(0);
     86: }

     88: /***********************************ivec.c*************************************/
     89: PetscErrorCode PCTFS_ivec_lxor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
     90: {
     92:   while (n--) {
     93:     *arg1=((*arg1 || *arg2) && !(*arg1 && *arg2));
     94:     arg1++;
     95:     arg2++;
     96:   }
     97:   return(0);
     98: }

    100: /***********************************ivec.c*************************************/
    101: PetscErrorCode PCTFS_ivec_xor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
    102: {
    104:   while (n--) *arg1++ ^= *arg2++;
    105:   return(0);
    106: }

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

    116: /***********************************ivec.c*************************************/
    117: PetscErrorCode PCTFS_ivec_lor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
    118: {
    120:   while (n--) {
    121:     *arg1 = (*arg1 || *arg2);
    122:     arg1++;
    123:     arg2++;
    124:   }
    125:   return(0);
    126: }

    128: /***********************************ivec.c*************************************/
    129: PetscErrorCode PCTFS_ivec_and(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
    130: {
    132:   while (n--) *arg1++ &= *arg2++;
    133:   return(0);
    134: }

    136: /***********************************ivec.c*************************************/
    137: PetscErrorCode PCTFS_ivec_land(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
    138: {
    140:   while (n--) {
    141:     *arg1 = (*arg1 && *arg2);
    142:     arg1++;
    143:     arg2++;
    144:   }
    145:   return(0);
    146: }

    148: /***********************************ivec.c*************************************/
    149: PetscErrorCode PCTFS_ivec_and3(PetscInt *arg1,  PetscInt *arg2,  PetscInt *arg3, PetscInt n)
    150: {
    152:   while (n--) *arg1++ = (*arg2++ & *arg3++);
    153:   return(0);
    154: }

    156: /***********************************ivec.c*************************************/
    157: PetscInt PCTFS_ivec_sum(PetscInt *arg1,  PetscInt n)
    158: {
    159:   PetscInt tmp = 0;
    160:   while (n--) tmp += *arg1++;
    161:   return(tmp);
    162: }

    164: /***********************************ivec.c*************************************/
    165: PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2,  PetscInt n,  PetscInt *arg3)
    166: {
    167:   PetscInt i, j, type;

    170:   /* LATER: if we're really motivated we can sort and then unsort */
    171:   for (i=0; i<n;) {
    172:     /* clump 'em for now */
    173:     j    =i+1;
    174:     type = arg3[i];
    175:     while ((j<n)&&(arg3[j]==type)) j++;

    177:     /* how many together */
    178:     j -= i;

    180:     /* call appropriate ivec function */
    181:     if (type == GL_MAX)        PCTFS_ivec_max(arg1,arg2,j);
    182:     else if (type == GL_MIN)   PCTFS_ivec_min(arg1,arg2,j);
    183:     else if (type == GL_MULT)  PCTFS_ivec_mult(arg1,arg2,j);
    184:     else if (type == GL_ADD)   PCTFS_ivec_add(arg1,arg2,j);
    185:     else if (type == GL_B_XOR) PCTFS_ivec_xor(arg1,arg2,j);
    186:     else if (type == GL_B_OR)  PCTFS_ivec_or(arg1,arg2,j);
    187:     else if (type == GL_B_AND) PCTFS_ivec_and(arg1,arg2,j);
    188:     else if (type == GL_L_XOR) PCTFS_ivec_lxor(arg1,arg2,j);
    189:     else if (type == GL_L_OR)  PCTFS_ivec_lor(arg1,arg2,j);
    190:     else if (type == GL_L_AND) PCTFS_ivec_land(arg1,arg2,j);
    191:     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_ivec_non_uniform()!");

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

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

    213:   /* catch all ... not good if we get here */
    214:   return(NULL);
    215: }

    217: /******************************************************************************/
    218: PetscErrorCode PCTFS_ivec_sort(PetscInt *ar,  PetscInt size)
    219: {
    220:   PetscInt *pi, *pj, temp;
    221:   PetscInt **top_a = (PetscInt**) offset_stack;
    222:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;

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

    229:   /* do until we're done ... return when stack is exhausted */
    230:   for (;;) {
    231:     /* if list is large enough use quicksort partition exchange code */
    232:     if (size > SORT_OPT) {
    233:       /* start up pointer at element 1 and down at size     */
    234:       pi = ar+1;
    235:       pj = ar+size;

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

    240:       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
    241:       /* note ==> pivot_value in index 0                   */
    242:       if (*pi > *pj) { SWAP(*pi,*pj) }
    243:       if (*ar > *pj) { SWAP(*ar,*pj) }
    244:       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) }

    246:       /* partition about pivot_value ...                              */
    247:       /* note lists of length 2 are not guaranteed to be sorted */
    248:       for (;;) {
    249:         /* walk up ... and down ... swap if equal to pivot! */
    250:         do pi++; while (*pi<*ar);
    251:         do pj--; while (*pj>*ar);

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

    256:         /* else swap */
    257:         SWAP(*pi,*pj)
    258:       }

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

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

    266:       /* push right hand child iff length > 1 */
    267:       if ((*top_s = size-((PetscInt) (pi-ar)))) {
    268:         *(top_a++) = pi;
    269:         size      -= *top_s+2;
    270:         top_s++;
    271:       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
    272:       else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */
    273:         ar   = *(--top_a);
    274:         size = *(--top_s);
    275:       }
    276:     } else { /* else sort small list directly then pop another off stack */

    278:       /* insertion sort for bottom */
    279:       for (pj=ar+1; pj<=ar+size; pj++) {
    280:         temp = *pj;
    281:         for (pi=pj-1; pi>=ar; pi--) {
    282:           if (*pi <= temp) break;
    283:           *(pi+1)=*pi;
    284:         }
    285:         *(pi+1)=temp;
    286:       }

    288:       /* check to see if stack is exhausted ==> DONE */
    289:       if (top_s==bottom_s) return(0);

    291:       /* else pop another list from the stack */
    292:       ar   = *(--top_a);
    293:       size = *(--top_s);
    294:     }
    295:   }
    296: }

    298: /******************************************************************************/
    299: PetscErrorCode PCTFS_ivec_sort_companion(PetscInt *ar,  PetscInt *ar2,  PetscInt size)
    300: {
    301:   PetscInt *pi, *pj, temp, temp2;
    302:   PetscInt **top_a = (PetscInt**)offset_stack;
    303:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
    304:   PetscInt *pi2, *pj2;
    305:   PetscInt mid;

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

    312:   /* do until we're done ... return when stack is exhausted */
    313:   for (;;) {

    315:     /* if list is large enough use quicksort partition exchange code */
    316:     if (size > SORT_OPT) {

    318:       /* start up pointer at element 1 and down at size     */
    319:       mid = size>>1;
    320:       pi  = ar+1;
    321:       pj  = ar+mid;
    322:       pi2 = ar2+1;
    323:       pj2 = ar2+mid;

    325:       /* find middle element in list and swap w/ element 1 */
    326:       SWAP(*pi,*pj)
    327:       SWAP(*pi2,*pj2)

    329:       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
    330:       /* note ==> pivot_value in index 0                   */
    331:       pj  = ar+size;
    332:       pj2 = ar2+size;
    333:       if (*pi > *pj) { SWAP(*pi,*pj) SWAP(*pi2,*pj2) }
    334:       if (*ar > *pj) { SWAP(*ar,*pj) SWAP(*ar2,*pj2) }
    335:       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1)) }

    337:       /* partition about pivot_value ...                              */
    338:       /* note lists of length 2 are not guaranteed to be sorted */
    339:       for (;;) {
    340:         /* walk up ... and down ... swap if equal to pivot! */
    341:         do { pi++; pi2++; } while (*pi<*ar);
    342:         do { pj--; pj2--; } while (*pj>*ar);

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

    347:         /* else swap */
    348:         SWAP(*pi,*pj)
    349:         SWAP(*pi2,*pj2)
    350:       }

    352:       /* place pivot_value in it's correct location */
    353:       SWAP(*ar,*pj)
    354:       SWAP(*ar2,*pj2)

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

    359:       /* push right hand child iff length > 1 */
    360:       if ((*top_s = size-((PetscInt) (pi-ar)))) {
    361:         *(top_a++) = pi;
    362:         *(top_a++) = pi2;
    363:         size      -= *top_s+2;
    364:         top_s++;
    365:       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
    366:       else {  /* might as well pop - note NR_OPT >=2 ==> we're ok! */
    367:         ar2  = *(--top_a);
    368:         ar   = *(--top_a);
    369:         size = *(--top_s);
    370:       }
    371:     } else { /* else sort small list directly then pop another off stack */

    373:       /* insertion sort for bottom */
    374:       for (pj=ar+1, pj2=ar2+1; pj<=ar+size; pj++,pj2++) {
    375:         temp  = *pj;
    376:         temp2 = *pj2;
    377:         for (pi=pj-1,pi2=pj2-1; pi>=ar; pi--,pi2--) {
    378:           if (*pi <= temp) break;
    379:           *(pi+1) =*pi;
    380:           *(pi2+1)=*pi2;
    381:         }
    382:         *(pi+1) =temp;
    383:         *(pi2+1)=temp2;
    384:       }

    386:       /* check to see if stack is exhausted ==> DONE */
    387:       if (top_s==bottom_s) return(0);

    389:       /* else pop another list from the stack */
    390:       ar2  = *(--top_a);
    391:       ar   = *(--top_a);
    392:       size = *(--top_s);
    393:     }
    394:   }
    395: }

    397: /******************************************************************************/
    398: PetscErrorCode PCTFS_ivec_sort_companion_hack(PetscInt *ar,  PetscInt **ar2, PetscInt size)
    399: {
    400:   PetscInt *pi, *pj, temp, *ptr;
    401:   PetscInt **top_a = (PetscInt**)offset_stack;
    402:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
    403:   PetscInt **pi2, **pj2;
    404:   PetscInt mid;

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

    411:   /* do until we're done ... return when stack is exhausted */
    412:   for (;;) {

    414:     /* if list is large enough use quicksort partition exchange code */
    415:     if (size > SORT_OPT) {

    417:       /* start up pointer at element 1 and down at size     */
    418:       mid = size>>1;
    419:       pi  = ar+1;
    420:       pj  = ar+mid;
    421:       pi2 = ar2+1;
    422:       pj2 = ar2+mid;

    424:       /* find middle element in list and swap w/ element 1 */
    425:       SWAP(*pi,*pj)
    426:       P_SWAP(*pi2,*pj2)

    428:       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
    429:       /* note ==> pivot_value in index 0                   */
    430:       pj  = ar+size;
    431:       pj2 = ar2+size;
    432:       if (*pi > *pj) { SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) }
    433:       if (*ar > *pj) { SWAP(*ar,*pj) P_SWAP(*ar2,*pj2) }
    434:       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1)) }

    436:       /* partition about pivot_value ...                              */
    437:       /* note lists of length 2 are not guaranteed to be sorted */
    438:       for (;;) {

    440:         /* walk up ... and down ... swap if equal to pivot! */
    441:         do {pi++; pi2++;} while (*pi<*ar);
    442:         do {pj--; pj2--;} while (*pj>*ar);

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

    447:         /* else swap */
    448:         SWAP(*pi,*pj)
    449:         P_SWAP(*pi2,*pj2)
    450:       }

    452:       /* place pivot_value in it's correct location */
    453:       SWAP(*ar,*pj)
    454:       P_SWAP(*ar2,*pj2)

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

    459:       /* push right hand child iff length > 1 */
    460:       if ((*top_s = size-((PetscInt) (pi-ar)))) {
    461:         *(top_a++) = pi;
    462:         *(top_a++) = (PetscInt*) pi2;
    463:         size      -= *top_s+2;
    464:         top_s++;
    465:       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
    466:       else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */
    467:         ar2  = (PetscInt**) *(--top_a);
    468:         ar   = *(--top_a);
    469:         size = *(--top_s);
    470:       }
    471:     } else  { /* else sort small list directly then pop another off stack */
    472:       /* insertion sort for bottom */
    473:       for (pj=ar+1, pj2=ar2+1; pj<=ar+size; pj++,pj2++) {
    474:         temp = *pj;
    475:         ptr  = *pj2;
    476:         for (pi=pj-1,pi2=pj2-1; pi>=ar; pi--,pi2--) {
    477:           if (*pi <= temp) break;
    478:           *(pi+1) =*pi;
    479:           *(pi2+1)=*pi2;
    480:         }
    481:         *(pi+1) =temp;
    482:         *(pi2+1)=ptr;
    483:       }

    485:       /* check to see if stack is exhausted ==> DONE */
    486:       if (top_s==bottom_s) return(0);

    488:       /* else pop another list from the stack */
    489:       ar2  = (PetscInt**)*(--top_a);
    490:       ar   = *(--top_a);
    491:       size = *(--top_s);
    492:     }
    493:   }
    494: }

    496: /******************************************************************************/
    497: PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type)
    498: {
    500:   if (type == SORT_INTEGER) {
    501:     if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size);
    502:     else PCTFS_ivec_sort((PetscInt*)ar1,size);
    503:   } else if (type == SORT_INT_PTR) {
    504:     if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt**)ar2,size);
    505:     else PCTFS_ivec_sort((PetscInt*)ar1,size);
    506:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!");
    507:   return(0);
    508: }

    510: /***********************************ivec.c*************************************/
    511: PetscInt PCTFS_ivec_linear_search(PetscInt item,  PetscInt *list,  PetscInt n)
    512: {
    513:   PetscInt tmp = n-1;

    516:   while (n--) {
    517:     if (*list++ == item) return(tmp-n);
    518:   }
    519:   return(-1);
    520: }

    522: /***********************************ivec.c*************************************/
    523: PetscInt PCTFS_ivec_binary_search(PetscInt item,  PetscInt *list,  PetscInt rh)
    524: {
    525:   PetscInt mid, lh=0;

    527:   rh--;
    528:   while (lh<=rh) {
    529:     mid = (lh+rh)>>1;
    530:     if (*(list+mid) == item) return(mid);
    531:     if (*(list+mid) > item) rh = mid-1;
    532:     else lh = mid+1;
    533:   }
    534:   return(-1);
    535: }

    537: /*********************************ivec.c*************************************/
    538: PetscErrorCode PCTFS_rvec_copy(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    539: {
    541:   while (n--) *arg1++ = *arg2++;
    542:   return(0);
    543: }

    545: /*********************************ivec.c*************************************/
    546: PetscErrorCode PCTFS_rvec_zero(PetscScalar *arg1,  PetscInt n)
    547: {
    549:   while (n--) *arg1++ = 0.0;
    550:   return(0);
    551: }

    553: /***********************************ivec.c*************************************/
    554: PetscErrorCode PCTFS_rvec_one(PetscScalar *arg1,  PetscInt n)
    555: {
    557:   while (n--) *arg1++ = 1.0;
    558:   return(0);
    559: }

    561: /***********************************ivec.c*************************************/
    562: PetscErrorCode PCTFS_rvec_set(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
    563: {
    565:   while (n--) *arg1++ = arg2;
    566:   return(0);
    567: }

    569: /***********************************ivec.c*************************************/
    570: PetscErrorCode PCTFS_rvec_scale(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
    571: {
    573:   while (n--) *arg1++ *= arg2;
    574:   return(0);
    575: }

    577: /*********************************ivec.c*************************************/
    578: PetscErrorCode PCTFS_rvec_add(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    579: {
    581:   while (n--) *arg1++ += *arg2++;
    582:   return(0);
    583: }

    585: /*********************************ivec.c*************************************/
    586: PetscErrorCode PCTFS_rvec_mult(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    587: {
    589:   while (n--) *arg1++ *= *arg2++;
    590:   return(0);
    591: }

    593: /*********************************ivec.c*************************************/
    594: PetscErrorCode PCTFS_rvec_max(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    595: {
    597:   while (n--) {
    598:     *arg1 = PetscMax(*arg1,*arg2);
    599:     arg1++;
    600:     arg2++;
    601:   }
    602:   return(0);
    603: }

    605: /*********************************ivec.c*************************************/
    606: PetscErrorCode PCTFS_rvec_max_abs(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    607: {
    609:   while (n--) {
    610:     *arg1 = MAX_FABS(*arg1,*arg2);
    611:     arg1++;
    612:     arg2++;
    613:   }
    614:   return(0);
    615: }

    617: /*********************************ivec.c*************************************/
    618: PetscErrorCode PCTFS_rvec_min(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    619: {
    621:   while (n--) {
    622:     *arg1 = PetscMin(*arg1,*arg2);
    623:     arg1++;
    624:     arg2++;
    625:   }
    626:   return(0);
    627: }

    629: /*********************************ivec.c*************************************/
    630: PetscErrorCode PCTFS_rvec_min_abs(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    631: {
    633:   while (n--) {
    634:     *arg1 = MIN_FABS(*arg1,*arg2);
    635:     arg1++;
    636:     arg2++;
    637:   }
    638:   return(0);
    639: }

    641: /*********************************ivec.c*************************************/
    642: PetscErrorCode PCTFS_rvec_exists(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
    643: {
    645:   while (n--) {
    646:     *arg1 = EXISTS(*arg1,*arg2);
    647:     arg1++;
    648:     arg2++;
    649:   }
    650:   return(0);
    651: }

    653: /***********************************ivec.c*************************************/
    654: PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2,  PetscInt n,  PetscInt *arg3)
    655: {
    656:   PetscInt i, j, type;

    659:   /* LATER: if we're really motivated we can sort and then unsort */
    660:   for (i=0; i<n;) {

    662:     /* clump 'em for now */
    663:     j    =i+1;
    664:     type = arg3[i];
    665:     while ((j<n)&&(arg3[j]==type)) j++;

    667:     /* how many together */
    668:     j -= i;

    670:     /* call appropriate ivec function */
    671:     if (type == GL_MAX)          PCTFS_rvec_max(arg1,arg2,j);
    672:     else if (type == GL_MIN)     PCTFS_rvec_min(arg1,arg2,j);
    673:     else if (type == GL_MULT)    PCTFS_rvec_mult(arg1,arg2,j);
    674:     else if (type == GL_ADD)     PCTFS_rvec_add(arg1,arg2,j);
    675:     else if (type == GL_MAX_ABS) PCTFS_rvec_max_abs(arg1,arg2,j);
    676:     else if (type == GL_MIN_ABS) PCTFS_rvec_min_abs(arg1,arg2,j);
    677:     else if (type == GL_EXISTS)  PCTFS_rvec_exists(arg1,arg2,j);
    678:     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_rvec_non_uniform()!");

    680:     arg1+=j; arg2+=j; i+=j;
    681:   }
    682:   return(0);
    683: }

    685: /***********************************ivec.c*************************************/
    686: vfp PCTFS_rvec_fct_addr(PetscInt type)
    687: {
    688:   if (type == NON_UNIFORM)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_non_uniform);
    689:   else if (type == GL_MAX)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_max);
    690:   else if (type == GL_MIN)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_min);
    691:   else if (type == GL_MULT)    return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_mult);
    692:   else if (type == GL_ADD)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_add);
    693:   else if (type == GL_MAX_ABS) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_max_abs);
    694:   else if (type == GL_MIN_ABS) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_min_abs);
    695:   else if (type == GL_EXISTS)  return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_exists);

    697:   /* catch all ... not good if we get here */
    698:   return(NULL);
    699: }