Actual source code: ivec.c


  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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

170:   va_start(ap, n);
171:   arg3 = va_arg(ap, PetscInt*);
172:   va_end(ap);

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

181:     /* how many together */
182:     j -= i;

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

197:     arg1+=j; arg2+=j; i+=j;
198:   }
199:   return(0);
200: }

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

217:   /* catch all ... not good if we get here */
218:   return(NULL);
219: }

221: /******************************************************************************/
222: PetscErrorCode PCTFS_ivec_sort(PetscInt *ar,  PetscInt size)
223: {
224:   PetscInt *pi, *pj, temp;
225:   PetscInt **top_a = (PetscInt**) offset_stack;
226:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;

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

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

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

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

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

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

260:         /* else swap */
261:         SWAP(*pi,*pj)
262:       }

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

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

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

282:       /* insertion sort for bottom */
283:       for (pj=ar+1; pj<=ar+size; pj++) {
284:         temp = *pj;
285:         for (pi=pj-1; pi>=ar; pi--) {
286:           if (*pi <= temp) break;
287:           *(pi+1)=*pi;
288:         }
289:         *(pi+1)=temp;
290:       }

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

295:       /* else pop another list from the stack */
296:       ar   = *(--top_a);
297:       size = *(--top_s);
298:     }
299:   }
300: }

302: /******************************************************************************/
303: PetscErrorCode PCTFS_ivec_sort_companion(PetscInt *ar,  PetscInt *ar2,  PetscInt size)
304: {
305:   PetscInt *pi, *pj, temp, temp2;
306:   PetscInt **top_a = (PetscInt**)offset_stack;
307:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
308:   PetscInt *pi2, *pj2;
309:   PetscInt mid;

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

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

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

322:       /* start up pointer at element 1 and down at size     */
323:       mid = size>>1;
324:       pi  = ar+1;
325:       pj  = ar+mid;
326:       pi2 = ar2+1;
327:       pj2 = ar2+mid;

329:       /* find middle element in list and swap w/ element 1 */
330:       SWAP(*pi,*pj)
331:       SWAP(*pi2,*pj2)

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

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

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

351:         /* else swap */
352:         SWAP(*pi,*pj)
353:         SWAP(*pi2,*pj2)
354:       }

356:       /* place pivot_value in it's correct location */
357:       SWAP(*ar,*pj)
358:       SWAP(*ar2,*pj2)

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

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

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

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

393:       /* else pop another list from the stack */
394:       ar2  = *(--top_a);
395:       ar   = *(--top_a);
396:       size = *(--top_s);
397:     }
398:   }
399: }

401: /******************************************************************************/
402: PetscErrorCode PCTFS_ivec_sort_companion_hack(PetscInt *ar,  PetscInt **ar2, PetscInt size)
403: {
404:   PetscInt *pi, *pj, temp, *ptr;
405:   PetscInt **top_a = (PetscInt**)offset_stack;
406:   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
407:   PetscInt **pi2, **pj2;
408:   PetscInt mid;

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

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

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

421:       /* start up pointer at element 1 and down at size     */
422:       mid = size>>1;
423:       pi  = ar+1;
424:       pj  = ar+mid;
425:       pi2 = ar2+1;
426:       pj2 = ar2+mid;

428:       /* find middle element in list and swap w/ element 1 */
429:       SWAP(*pi,*pj)
430:       P_SWAP(*pi2,*pj2)

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

440:       /* partition about pivot_value ...                              */
441:       /* note lists of length 2 are not guaranteed to be sorted */
442:       for (;;) {

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

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

451:         /* else swap */
452:         SWAP(*pi,*pj)
453:         P_SWAP(*pi2,*pj2)
454:       }

456:       /* place pivot_value in it's correct location */
457:       SWAP(*ar,*pj)
458:       P_SWAP(*ar2,*pj2)

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

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

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

492:       /* else pop another list from the stack */
493:       ar2  = (PetscInt**)*(--top_a);
494:       ar   = *(--top_a);
495:       size = *(--top_s);
496:     }
497:   }
498: }

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

514: /***********************************ivec.c*************************************/
515: PetscInt PCTFS_ivec_linear_search(PetscInt item,  PetscInt *list,  PetscInt n)
516: {
517:   PetscInt tmp = n-1;

519:   while (n--) {
520:     if (*list++ == item) return(tmp-n);
521:   }
522:   return(-1);
523: }

525: /***********************************ivec.c*************************************/
526: PetscInt PCTFS_ivec_binary_search(PetscInt item,  PetscInt *list,  PetscInt rh)
527: {
528:   PetscInt mid, lh=0;

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

540: /*********************************ivec.c*************************************/
541: PetscErrorCode PCTFS_rvec_copy(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
542: {
544:   while (n--) *arg1++ = *arg2++;
545:   return(0);
546: }

548: /*********************************ivec.c*************************************/
549: PetscErrorCode PCTFS_rvec_zero(PetscScalar *arg1,  PetscInt n)
550: {
552:   while (n--) *arg1++ = 0.0;
553:   return(0);
554: }

556: /***********************************ivec.c*************************************/
557: PetscErrorCode PCTFS_rvec_one(PetscScalar *arg1,  PetscInt n)
558: {
560:   while (n--) *arg1++ = 1.0;
561:   return(0);
562: }

564: /***********************************ivec.c*************************************/
565: PetscErrorCode PCTFS_rvec_set(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
566: {
568:   while (n--) *arg1++ = arg2;
569:   return(0);
570: }

572: /***********************************ivec.c*************************************/
573: PetscErrorCode PCTFS_rvec_scale(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
574: {
576:   while (n--) *arg1++ *= arg2;
577:   return(0);
578: }

580: /*********************************ivec.c*************************************/
581: PetscErrorCode PCTFS_rvec_add(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
582: {
584:   while (n--) *arg1++ += *arg2++;
585:   return(0);
586: }

588: /*********************************ivec.c*************************************/
589: PetscErrorCode PCTFS_rvec_mult(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
590: {
592:   while (n--) *arg1++ *= *arg2++;
593:   return(0);
594: }

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

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

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

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

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

656: /***********************************ivec.c*************************************/
657: PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2,  PetscInt n,  PetscInt *arg3)
658: {
659:   PetscInt i, j, type;

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

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

670:     /* how many together */
671:     j -= i;

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

683:     arg1+=j; arg2+=j; i+=j;
684:   }
685:   return(0);
686: }

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

700:   /* catch all ... not good if we get here */
701:   return(NULL);
702: }