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: {
 37:   while (n--) *arg1++ = 0;
 38:   return 0;
 39: }

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

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

 59: /***********************************ivec.c*************************************/
 60: PetscErrorCode PCTFS_ivec_min(PetscInt *arg1, PetscInt *arg2, PetscInt n)
 61: {
 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: {
 73:   while (n--) *arg1++ *= *arg2++;
 74:   return 0;
 75: }

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

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

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

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

109: /***********************************ivec.c*************************************/
110: PetscErrorCode PCTFS_ivec_lor(PetscInt *arg1, PetscInt *arg2, PetscInt n)
111: {
112:   while (n--) {
113:     *arg1 = (*arg1 || *arg2);
114:     arg1++;
115:     arg2++;
116:   }
117:   return 0;
118: }

120: /***********************************ivec.c*************************************/
121: PetscErrorCode PCTFS_ivec_and(PetscInt *arg1, PetscInt *arg2, PetscInt n)
122: {
123:   while (n--) *arg1++ &= *arg2++;
124:   return 0;
125: }

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

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

145: /***********************************ivec.c*************************************/
146: PetscInt PCTFS_ivec_sum(PetscInt *arg1, PetscInt n)
147: {
148:   PetscInt tmp = 0;
149:   while (n--) tmp += *arg1++;
150:   return (tmp);
151: }

153: /***********************************ivec.c*************************************/
154: PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2, PetscInt n, ...)
155: {
156:   PetscInt  i, j, type;
157:   PetscInt *arg3;
158:   va_list   ap;

160:   va_start(ap, n);
161:   arg3 = va_arg(ap, PetscInt *);
162:   va_end(ap);

164:   /* LATER: if we're really motivated we can sort and then unsort */
165:   for (i = 0; i < n;) {
166:     /* clump 'em for now */
167:     j    = i + 1;
168:     type = arg3[i];
169:     while ((j < n) && (arg3[j] == type)) j++;

171:     /* how many together */
172:     j -= i;

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

187:     arg1 += j;
188:     arg2 += j;
189:     i += j;
190:   }
191:   return 0;
192: }

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

209:   /* catch all ... not good if we get here */
210:   return (NULL);
211: }

213: /******************************************************************************/
214: PetscErrorCode PCTFS_ivec_sort(PetscInt *ar, PetscInt size)
215: {
216:   PetscInt  *pi, *pj, temp;
217:   PetscInt **top_a = (PetscInt **)offset_stack;
218:   PetscInt  *top_s = size_stack, *bottom_s = size_stack;

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

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

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

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

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

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)
272:         ; /* set up for next loop iff there is something to do */
273:       else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 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;

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

311:   /* do until we're done ... return when stack is exhausted */
312:   for (;;) {
313:     /* if list is large enough use quicksort partition exchange code */
314:     if (size > SORT_OPT) {
315:       /* start up pointer at element 1 and down at size     */
316:       mid = size >> 1;
317:       pi  = ar + 1;
318:       pj  = ar + mid;
319:       pi2 = ar2 + 1;
320:       pj2 = ar2 + mid;

322:       /* find middle element in list and swap w/ element 1 */
323:       SWAP(*pi, *pj)
324:       SWAP(*pi2, *pj2)

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

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 {
342:           pi++;
343:           pi2++;
344:         } while (*pi < *ar);
345:         do {
346:           pj--;
347:           pj2--;
348:         } while (*pj > *ar);

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

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

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

362:       /* test stack_size to see if we've exhausted our stack */

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

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

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

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

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

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

416:   /* do until we're done ... return when stack is exhausted */
417:   for (;;) {
418:     /* if list is large enough use quicksort partition exchange code */
419:     if (size > SORT_OPT) {
420:       /* start up pointer at element 1 and down at size     */
421:       mid = size >> 1;
422:       pi  = ar + 1;
423:       pj  = ar + mid;
424:       pi2 = ar2 + 1;
425:       pj2 = ar2 + mid;

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

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

442:       /* partition about pivot_value ...                              */
443:       /* note lists of length 2 are not guaranteed to be sorted */
444:       for (;;) {
445:         /* walk up ... and down ... swap if equal to pivot! */
446:         do {
447:           pi++;
448:           pi2++;
449:         } while (*pi < *ar);
450:         do {
451:           pj--;
452:           pj2--;
453:         } while (*pj > *ar);

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

458:         /* else swap */
459:         SWAP(*pi, *pj)
460:         P_SWAP(*pi2, *pj2)
461:       }

463:       /* place pivot_value in it's correct location */
464:       SWAP(*ar, *pj)
465:       P_SWAP(*ar2, *pj2)

467:       /* test stack_size to see if we've exhausted our stack */

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

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

499:       /* else pop another list from the stack */
500:       ar2  = (PetscInt **)*(--top_a);
501:       ar   = *(--top_a);
502:       size = *(--top_s);
503:     }
504:   }
505: }

507: /******************************************************************************/
508: PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type)
509: {
510:   if (type == SORT_INTEGER) {
511:     if (ar2) PCTFS_ivec_sort_companion((PetscInt *)ar1, (PetscInt *)ar2, size);
512:     else PCTFS_ivec_sort((PetscInt *)ar1, size);
513:   } else if (type == SORT_INT_PTR) {
514:     if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt *)ar1, (PetscInt **)ar2, size);
515:     else PCTFS_ivec_sort((PetscInt *)ar1, size);
516:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCTFS_SMI_sort only does SORT_INTEGER!");
517:   return 0;
518: }

520: /***********************************ivec.c*************************************/
521: PetscInt PCTFS_ivec_linear_search(PetscInt item, PetscInt *list, PetscInt n)
522: {
523:   PetscInt tmp = n - 1;

525:   while (n--) {
526:     if (*list++ == item) return (tmp - n);
527:   }
528:   return (-1);
529: }

531: /***********************************ivec.c*************************************/
532: PetscInt PCTFS_ivec_binary_search(PetscInt item, PetscInt *list, PetscInt rh)
533: {
534:   PetscInt mid, lh = 0;

536:   rh--;
537:   while (lh <= rh) {
538:     mid = (lh + rh) >> 1;
539:     if (*(list + mid) == item) return (mid);
540:     if (*(list + mid) > item) rh = mid - 1;
541:     else lh = mid + 1;
542:   }
543:   return (-1);
544: }

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

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

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

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

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

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

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

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

606: /*********************************ivec.c*************************************/
607: PetscErrorCode PCTFS_rvec_max_abs(PetscScalar *arg1, PetscScalar *arg2, PetscInt n)
608: {
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: {
620:   while (n--) {
621:     *arg1 = PetscMin(*arg1, *arg2);
622:     arg1++;
623:     arg2++;
624:   }
625:   return 0;
626: }

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

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

650: /***********************************ivec.c*************************************/
651: PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2, PetscInt n, PetscInt *arg3)
652: {
653:   PetscInt i, j, type;

655:   /* LATER: if we're really motivated we can sort and then unsort */
656:   for (i = 0; i < n;) {
657:     /* clump 'em for now */
658:     j    = i + 1;
659:     type = arg3[i];
660:     while ((j < n) && (arg3[j] == type)) j++;

662:     /* how many together */
663:     j -= i;

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

675:     arg1 += j;
676:     arg2 += j;
677:     i += j;
678:   }
679:   return 0;
680: }

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

694:   /* catch all ... not good if we get here */
695:   return (NULL);
696: }