xref: /petsc/src/ksp/pc/impls/tfs/ivec.c (revision d57bb9db46f99bff4aaa690e4e9a7e55805ae415)
1 
2 /**********************************ivec.c**************************************
3 
4 Author: Henry M. Tufo III
5 
6 e-mail: hmt@cs.brown.edu
7 
8 snail-mail:
9 Division of Applied Mathematics
10 Brown University
11 Providence, RI 02912
12 
13 Last Modification:
14 6.21.97
15 ***********************************ivec.c*************************************/
16 
17 #include <../src/ksp/pc/impls/tfs/tfs.h>
18 
19 /* sorting args ivec.c ivec.c ... */
20 #define   SORT_OPT     6
21 #define   SORT_STACK   50000
22 
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];
26 
27 /***********************************ivec.c*************************************/
28 PetscInt *PCTFS_ivec_copy(PetscInt *arg1, PetscInt *arg2, PetscInt n)
29 {
30   while (n--) *arg1++ = *arg2++;
31   return(arg1);
32 }
33 
34 /***********************************ivec.c*************************************/
35 PetscErrorCode PCTFS_ivec_zero(PetscInt *arg1, PetscInt n)
36 {
37   PetscFunctionBegin;
38   while (n--) *arg1++ = 0;
39   PetscFunctionReturn(0);
40 }
41 
42 /***********************************ivec.c*************************************/
43 PetscErrorCode PCTFS_ivec_set(PetscInt *arg1, PetscInt arg2, PetscInt n)
44 {
45   PetscFunctionBegin;
46   while (n--) *arg1++ = arg2;
47   PetscFunctionReturn(0);
48 }
49 
50 /***********************************ivec.c*************************************/
51 PetscErrorCode PCTFS_ivec_max(PetscInt *arg1, PetscInt *arg2, PetscInt n)
52 {
53   PetscFunctionBegin;
54   while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; }
55   PetscFunctionReturn(0);
56 }
57 
58 /***********************************ivec.c*************************************/
59 PetscErrorCode PCTFS_ivec_min(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
60 {
61   PetscFunctionBegin;
62   while (n--) {
63     *(arg1) = PetscMin(*arg1,*arg2);
64     arg1++;
65     arg2++;
66   }
67   PetscFunctionReturn(0);
68 }
69 
70 /***********************************ivec.c*************************************/
71 PetscErrorCode PCTFS_ivec_mult(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
72 {
73   PetscFunctionBegin;
74   while (n--) *arg1++ *= *arg2++;
75   PetscFunctionReturn(0);
76 }
77 
78 /***********************************ivec.c*************************************/
79 PetscErrorCode PCTFS_ivec_add(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
80 {
81   PetscFunctionBegin;
82   while (n--) *arg1++ += *arg2++;
83   PetscFunctionReturn(0);
84 }
85 
86 /***********************************ivec.c*************************************/
87 PetscErrorCode PCTFS_ivec_lxor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
88 {
89   PetscFunctionBegin;
90   while (n--) {
91     *arg1=((*arg1 || *arg2) && !(*arg1 && *arg2));
92     arg1++;
93     arg2++;
94   }
95   PetscFunctionReturn(0);
96 }
97 
98 /***********************************ivec.c*************************************/
99 PetscErrorCode PCTFS_ivec_xor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
100 {
101   PetscFunctionBegin;
102   while (n--) *arg1++ ^= *arg2++;
103   PetscFunctionReturn(0);
104 }
105 
106 /***********************************ivec.c*************************************/
107 PetscErrorCode PCTFS_ivec_or(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
108 {
109   PetscFunctionBegin;
110   while (n--) *arg1++ |= *arg2++;
111   PetscFunctionReturn(0);
112 }
113 
114 /***********************************ivec.c*************************************/
115 PetscErrorCode PCTFS_ivec_lor(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
116 {
117   PetscFunctionBegin;
118   while (n--) {
119     *arg1 = (*arg1 || *arg2);
120     arg1++;
121     arg2++;
122   }
123   PetscFunctionReturn(0);
124 }
125 
126 /***********************************ivec.c*************************************/
127 PetscErrorCode PCTFS_ivec_and(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
128 {
129   PetscFunctionBegin;
130   while (n--) *arg1++ &= *arg2++;
131   PetscFunctionReturn(0);
132 }
133 
134 /***********************************ivec.c*************************************/
135 PetscErrorCode PCTFS_ivec_land(PetscInt *arg1,  PetscInt *arg2,  PetscInt n)
136 {
137   PetscFunctionBegin;
138   while (n--) {
139     *arg1 = (*arg1 && *arg2);
140     arg1++;
141     arg2++;
142   }
143   PetscFunctionReturn(0);
144 }
145 
146 /***********************************ivec.c*************************************/
147 PetscErrorCode PCTFS_ivec_and3(PetscInt *arg1,  PetscInt *arg2,  PetscInt *arg3, PetscInt n)
148 {
149   PetscFunctionBegin;
150   while (n--) *arg1++ = (*arg2++ & *arg3++);
151   PetscFunctionReturn(0);
152 }
153 
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 }
161 
162 /***********************************ivec.c*************************************/
163 PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2,  PetscInt n,  PetscInt *arg3)
164 {
165   PetscInt i, j, type;
166 
167   PetscFunctionBegin;
168   /* LATER: if we're really motivated we can sort and then unsort */
169   for (i=0; i<n;) {
170     /* clump 'em for now */
171     j    =i+1;
172     type = arg3[i];
173     while ((j<n)&&(arg3[j]==type)) j++;
174 
175     /* how many together */
176     j -= i;
177 
178     /* call appropriate ivec function */
179     if (type == GL_MAX)        PCTFS_ivec_max(arg1,arg2,j);
180     else if (type == GL_MIN)   PCTFS_ivec_min(arg1,arg2,j);
181     else if (type == GL_MULT)  PCTFS_ivec_mult(arg1,arg2,j);
182     else if (type == GL_ADD)   PCTFS_ivec_add(arg1,arg2,j);
183     else if (type == GL_B_XOR) PCTFS_ivec_xor(arg1,arg2,j);
184     else if (type == GL_B_OR)  PCTFS_ivec_or(arg1,arg2,j);
185     else if (type == GL_B_AND) PCTFS_ivec_and(arg1,arg2,j);
186     else if (type == GL_L_XOR) PCTFS_ivec_lxor(arg1,arg2,j);
187     else if (type == GL_L_OR)  PCTFS_ivec_lor(arg1,arg2,j);
188     else if (type == GL_L_AND) PCTFS_ivec_land(arg1,arg2,j);
189     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_ivec_non_uniform()!");
190 
191     arg1+=j; arg2+=j; i+=j;
192   }
193   PetscFunctionReturn(0);
194 }
195 
196 /***********************************ivec.c*************************************/
197 vfp PCTFS_ivec_fct_addr(PetscInt type)
198 {
199   if (type == NON_UNIFORM)   return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_non_uniform);
200   else if (type == GL_MAX)   return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_max);
201   else if (type == GL_MIN)   return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_min);
202   else if (type == GL_MULT)  return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_mult);
203   else if (type == GL_ADD)   return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_add);
204   else if (type == GL_B_XOR) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_xor);
205   else if (type == GL_B_OR)  return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_or);
206   else if (type == GL_B_AND) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_and);
207   else if (type == GL_L_XOR) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_lxor);
208   else if (type == GL_L_OR)  return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_lor);
209   else if (type == GL_L_AND) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_ivec_land);
210 
211   /* catch all ... not good if we get here */
212   return(NULL);
213 }
214 
215 /******************************************************************************/
216 PetscErrorCode PCTFS_ivec_sort(PetscInt *ar,  PetscInt size)
217 {
218   PetscInt *pi, *pj, temp;
219   PetscInt **top_a = (PetscInt**) offset_stack;
220   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
221 
222   PetscFunctionBegin;
223   /* we're really interested in the offset of the last element */
224   /* ==> length of the list is now size + 1                    */
225   size--;
226 
227   /* do until we're done ... return when stack is exhausted */
228   for (;;) {
229     /* if list is large enough use quicksort partition exchange code */
230     if (size > SORT_OPT) {
231       /* start up pointer at element 1 and down at size     */
232       pi = ar+1;
233       pj = ar+size;
234 
235       /* find middle element in list and swap w/ element 1 */
236       SWAP(*(ar+(size>>1)),*pi)
237 
238       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
239       /* note ==> pivot_value in index 0                   */
240       if (*pi > *pj) { SWAP(*pi,*pj) }
241       if (*ar > *pj) { SWAP(*ar,*pj) }
242       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) }
243 
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++; while (*pi<*ar);
249         do pj--; while (*pj>*ar);
250 
251         /* if we've crossed we're done */
252         if (pj<pi) break;
253 
254         /* else swap */
255         SWAP(*pi,*pj)
256       }
257 
258       /* place pivot_value in it's correct location */
259       SWAP(*ar,*pj)
260 
261       /* test stack_size to see if we've exhausted our stack */
262       if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort() :: STACK EXHAUSTED!!!");
263 
264       /* push right hand child iff length > 1 */
265       if ((*top_s = size-((PetscInt) (pi-ar)))) {
266         *(top_a++) = pi;
267         size      -= *top_s+2;
268         top_s++;
269       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
270       else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */
271         ar   = *(--top_a);
272         size = *(--top_s);
273       }
274     } else { /* else sort small list directly then pop another off stack */
275 
276       /* insertion sort for bottom */
277       for (pj=ar+1; pj<=ar+size; pj++) {
278         temp = *pj;
279         for (pi=pj-1; pi>=ar; pi--) {
280           if (*pi <= temp) break;
281           *(pi+1)=*pi;
282         }
283         *(pi+1)=temp;
284       }
285 
286       /* check to see if stack is exhausted ==> DONE */
287       if (top_s==bottom_s) PetscFunctionReturn(0);
288 
289       /* else pop another list from the stack */
290       ar   = *(--top_a);
291       size = *(--top_s);
292     }
293   }
294 }
295 
296 /******************************************************************************/
297 PetscErrorCode PCTFS_ivec_sort_companion(PetscInt *ar,  PetscInt *ar2,  PetscInt size)
298 {
299   PetscInt *pi, *pj, temp, temp2;
300   PetscInt **top_a = (PetscInt**)offset_stack;
301   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
302   PetscInt *pi2, *pj2;
303   PetscInt mid;
304 
305   PetscFunctionBegin;
306   /* we're really interested in the offset of the last element */
307   /* ==> length of the list is now size + 1                    */
308   size--;
309 
310   /* do until we're done ... return when stack is exhausted */
311   for (;;) {
312 
313     /* if list is large enough use quicksort partition exchange code */
314     if (size > SORT_OPT) {
315 
316       /* start up pointer at element 1 and down at size     */
317       mid = size>>1;
318       pi  = ar+1;
319       pj  = ar+mid;
320       pi2 = ar2+1;
321       pj2 = ar2+mid;
322 
323       /* find middle element in list and swap w/ element 1 */
324       SWAP(*pi,*pj)
325       SWAP(*pi2,*pj2)
326 
327       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
328       /* note ==> pivot_value in index 0                   */
329       pj  = ar+size;
330       pj2 = ar2+size;
331       if (*pi > *pj) { SWAP(*pi,*pj) SWAP(*pi2,*pj2) }
332       if (*ar > *pj) { SWAP(*ar,*pj) SWAP(*ar2,*pj2) }
333       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1)) }
334 
335       /* partition about pivot_value ...                              */
336       /* note lists of length 2 are not guaranteed to be sorted */
337       for (;;) {
338         /* walk up ... and down ... swap if equal to pivot! */
339         do { pi++; pi2++; } while (*pi<*ar);
340         do { pj--; pj2--; } while (*pj>*ar);
341 
342         /* if we've crossed we're done */
343         if (pj<pi) break;
344 
345         /* else swap */
346         SWAP(*pi,*pj)
347         SWAP(*pi2,*pj2)
348       }
349 
350       /* place pivot_value in it's correct location */
351       SWAP(*ar,*pj)
352       SWAP(*ar2,*pj2)
353 
354       /* test stack_size to see if we've exhausted our stack */
355       if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion() :: STACK EXHAUSTED!!!");
356 
357       /* push right hand child iff length > 1 */
358       if ((*top_s = size-((PetscInt) (pi-ar)))) {
359         *(top_a++) = pi;
360         *(top_a++) = pi2;
361         size      -= *top_s+2;
362         top_s++;
363       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
364       else {  /* might as well pop - note NR_OPT >=2 ==> we're ok! */
365         ar2  = *(--top_a);
366         ar   = *(--top_a);
367         size = *(--top_s);
368       }
369     } else { /* else sort small list directly then pop another off stack */
370 
371       /* insertion sort for bottom */
372       for (pj=ar+1, pj2=ar2+1; pj<=ar+size; pj++,pj2++) {
373         temp  = *pj;
374         temp2 = *pj2;
375         for (pi=pj-1,pi2=pj2-1; pi>=ar; pi--,pi2--) {
376           if (*pi <= temp) break;
377           *(pi+1) =*pi;
378           *(pi2+1)=*pi2;
379         }
380         *(pi+1) =temp;
381         *(pi2+1)=temp2;
382       }
383 
384       /* check to see if stack is exhausted ==> DONE */
385       if (top_s==bottom_s) PetscFunctionReturn(0);
386 
387       /* else pop another list from the stack */
388       ar2  = *(--top_a);
389       ar   = *(--top_a);
390       size = *(--top_s);
391     }
392   }
393 }
394 
395 /******************************************************************************/
396 PetscErrorCode PCTFS_ivec_sort_companion_hack(PetscInt *ar,  PetscInt **ar2, PetscInt size)
397 {
398   PetscInt *pi, *pj, temp, *ptr;
399   PetscInt **top_a = (PetscInt**)offset_stack;
400   PetscInt *top_s  = size_stack, *bottom_s = size_stack;
401   PetscInt **pi2, **pj2;
402   PetscInt mid;
403 
404   PetscFunctionBegin;
405   /* we're really interested in the offset of the last element */
406   /* ==> length of the list is now size + 1                    */
407   size--;
408 
409   /* do until we're done ... return when stack is exhausted */
410   for (;;) {
411 
412     /* if list is large enough use quicksort partition exchange code */
413     if (size > SORT_OPT) {
414 
415       /* start up pointer at element 1 and down at size     */
416       mid = size>>1;
417       pi  = ar+1;
418       pj  = ar+mid;
419       pi2 = ar2+1;
420       pj2 = ar2+mid;
421 
422       /* find middle element in list and swap w/ element 1 */
423       SWAP(*pi,*pj)
424       P_SWAP(*pi2,*pj2)
425 
426       /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
427       /* note ==> pivot_value in index 0                   */
428       pj  = ar+size;
429       pj2 = ar2+size;
430       if (*pi > *pj) { SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) }
431       if (*ar > *pj) { SWAP(*ar,*pj) P_SWAP(*ar2,*pj2) }
432       else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1)) }
433 
434       /* partition about pivot_value ...                              */
435       /* note lists of length 2 are not guaranteed to be sorted */
436       for (;;) {
437 
438         /* walk up ... and down ... swap if equal to pivot! */
439         do {pi++; pi2++;} while (*pi<*ar);
440         do {pj--; pj2--;} while (*pj>*ar);
441 
442         /* if we've crossed we're done */
443         if (pj<pi) break;
444 
445         /* else swap */
446         SWAP(*pi,*pj)
447         P_SWAP(*pi2,*pj2)
448       }
449 
450       /* place pivot_value in it's correct location */
451       SWAP(*ar,*pj)
452       P_SWAP(*ar2,*pj2)
453 
454       /* test stack_size to see if we've exhausted our stack */
455       if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion_hack() :: STACK EXHAUSTED!!!");
456 
457       /* push right hand child iff length > 1 */
458       if ((*top_s = size-((PetscInt) (pi-ar)))) {
459         *(top_a++) = pi;
460         *(top_a++) = (PetscInt*) pi2;
461         size      -= *top_s+2;
462         top_s++;
463       } else if (size -= *top_s+2) ;   /* set up for next loop iff there is something to do */
464       else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */
465         ar2  = (PetscInt**) *(--top_a);
466         ar   = *(--top_a);
467         size = *(--top_s);
468       }
469     } else  { /* else sort small list directly then pop another off stack */
470       /* insertion sort for bottom */
471       for (pj=ar+1, pj2=ar2+1; pj<=ar+size; pj++,pj2++) {
472         temp = *pj;
473         ptr  = *pj2;
474         for (pi=pj-1,pi2=pj2-1; pi>=ar; pi--,pi2--) {
475           if (*pi <= temp) break;
476           *(pi+1) =*pi;
477           *(pi2+1)=*pi2;
478         }
479         *(pi+1) =temp;
480         *(pi2+1)=ptr;
481       }
482 
483       /* check to see if stack is exhausted ==> DONE */
484       if (top_s==bottom_s) PetscFunctionReturn(0);
485 
486       /* else pop another list from the stack */
487       ar2  = (PetscInt**)*(--top_a);
488       ar   = *(--top_a);
489       size = *(--top_s);
490     }
491   }
492 }
493 
494 /******************************************************************************/
495 PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type)
496 {
497   PetscFunctionBegin;
498   if (type == SORT_INTEGER) {
499     if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size);
500     else PCTFS_ivec_sort((PetscInt*)ar1,size);
501   } else if (type == SORT_INT_PTR) {
502     if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt**)ar2,size);
503     else PCTFS_ivec_sort((PetscInt*)ar1,size);
504   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!");
505   PetscFunctionReturn(0);
506 }
507 
508 /***********************************ivec.c*************************************/
509 PetscInt PCTFS_ivec_linear_search(PetscInt item,  PetscInt *list,  PetscInt n)
510 {
511   PetscInt tmp = n-1;
512 
513   PetscFunctionBegin;
514   while (n--) {
515     if (*list++ == item) return(tmp-n);
516   }
517   return(-1);
518 }
519 
520 /***********************************ivec.c*************************************/
521 PetscInt PCTFS_ivec_binary_search(PetscInt item,  PetscInt *list,  PetscInt rh)
522 {
523   PetscInt mid, lh=0;
524 
525   rh--;
526   while (lh<=rh) {
527     mid = (lh+rh)>>1;
528     if (*(list+mid) == item) return(mid);
529     if (*(list+mid) > item) rh = mid-1;
530     else lh = mid+1;
531   }
532   return(-1);
533 }
534 
535 /*********************************ivec.c*************************************/
536 PetscErrorCode PCTFS_rvec_copy(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
537 {
538   PetscFunctionBegin;
539   while (n--) *arg1++ = *arg2++;
540   PetscFunctionReturn(0);
541 }
542 
543 /*********************************ivec.c*************************************/
544 PetscErrorCode PCTFS_rvec_zero(PetscScalar *arg1,  PetscInt n)
545 {
546   PetscFunctionBegin;
547   while (n--) *arg1++ = 0.0;
548   PetscFunctionReturn(0);
549 }
550 
551 /***********************************ivec.c*************************************/
552 PetscErrorCode PCTFS_rvec_one(PetscScalar *arg1,  PetscInt n)
553 {
554   PetscFunctionBegin;
555   while (n--) *arg1++ = 1.0;
556   PetscFunctionReturn(0);
557 }
558 
559 /***********************************ivec.c*************************************/
560 PetscErrorCode PCTFS_rvec_set(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
561 {
562   PetscFunctionBegin;
563   while (n--) *arg1++ = arg2;
564   PetscFunctionReturn(0);
565 }
566 
567 /***********************************ivec.c*************************************/
568 PetscErrorCode PCTFS_rvec_scale(PetscScalar *arg1,  PetscScalar arg2,  PetscInt n)
569 {
570   PetscFunctionBegin;
571   while (n--) *arg1++ *= arg2;
572   PetscFunctionReturn(0);
573 }
574 
575 /*********************************ivec.c*************************************/
576 PetscErrorCode PCTFS_rvec_add(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
577 {
578   PetscFunctionBegin;
579   while (n--) *arg1++ += *arg2++;
580   PetscFunctionReturn(0);
581 }
582 
583 /*********************************ivec.c*************************************/
584 PetscErrorCode PCTFS_rvec_mult(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
585 {
586   PetscFunctionBegin;
587   while (n--) *arg1++ *= *arg2++;
588   PetscFunctionReturn(0);
589 }
590 
591 /*********************************ivec.c*************************************/
592 PetscErrorCode PCTFS_rvec_max(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
593 {
594   PetscFunctionBegin;
595   while (n--) {
596     *arg1 = PetscMax(*arg1,*arg2);
597     arg1++;
598     arg2++;
599   }
600   PetscFunctionReturn(0);
601 }
602 
603 /*********************************ivec.c*************************************/
604 PetscErrorCode PCTFS_rvec_max_abs(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
605 {
606   PetscFunctionBegin;
607   while (n--) {
608     *arg1 = MAX_FABS(*arg1,*arg2);
609     arg1++;
610     arg2++;
611   }
612   PetscFunctionReturn(0);
613 }
614 
615 /*********************************ivec.c*************************************/
616 PetscErrorCode PCTFS_rvec_min(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
617 {
618   PetscFunctionBegin;
619   while (n--) {
620     *arg1 = PetscMin(*arg1,*arg2);
621     arg1++;
622     arg2++;
623   }
624   PetscFunctionReturn(0);
625 }
626 
627 /*********************************ivec.c*************************************/
628 PetscErrorCode PCTFS_rvec_min_abs(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
629 {
630   PetscFunctionBegin;
631   while (n--) {
632     *arg1 = MIN_FABS(*arg1,*arg2);
633     arg1++;
634     arg2++;
635   }
636   PetscFunctionReturn(0);
637 }
638 
639 /*********************************ivec.c*************************************/
640 PetscErrorCode PCTFS_rvec_exists(PetscScalar *arg1,  PetscScalar *arg2,  PetscInt n)
641 {
642   PetscFunctionBegin;
643   while (n--) {
644     *arg1 = EXISTS(*arg1,*arg2);
645     arg1++;
646     arg2++;
647   }
648   PetscFunctionReturn(0);
649 }
650 
651 /***********************************ivec.c*************************************/
652 PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2,  PetscInt n,  PetscInt *arg3)
653 {
654   PetscInt i, j, type;
655 
656   PetscFunctionBegin;
657   /* LATER: if we're really motivated we can sort and then unsort */
658   for (i=0; i<n;) {
659 
660     /* clump 'em for now */
661     j    =i+1;
662     type = arg3[i];
663     while ((j<n)&&(arg3[j]==type)) j++;
664 
665     /* how many together */
666     j -= i;
667 
668     /* call appropriate ivec function */
669     if (type == GL_MAX)          PCTFS_rvec_max(arg1,arg2,j);
670     else if (type == GL_MIN)     PCTFS_rvec_min(arg1,arg2,j);
671     else if (type == GL_MULT)    PCTFS_rvec_mult(arg1,arg2,j);
672     else if (type == GL_ADD)     PCTFS_rvec_add(arg1,arg2,j);
673     else if (type == GL_MAX_ABS) PCTFS_rvec_max_abs(arg1,arg2,j);
674     else if (type == GL_MIN_ABS) PCTFS_rvec_min_abs(arg1,arg2,j);
675     else if (type == GL_EXISTS)  PCTFS_rvec_exists(arg1,arg2,j);
676     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_rvec_non_uniform()!");
677 
678     arg1+=j; arg2+=j; i+=j;
679   }
680   PetscFunctionReturn(0);
681 }
682 
683 /***********************************ivec.c*************************************/
684 vfp PCTFS_rvec_fct_addr(PetscInt type)
685 {
686   if (type == NON_UNIFORM)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_non_uniform);
687   else if (type == GL_MAX)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_max);
688   else if (type == GL_MIN)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_min);
689   else if (type == GL_MULT)    return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_mult);
690   else if (type == GL_ADD)     return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_add);
691   else if (type == GL_MAX_ABS) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_max_abs);
692   else if (type == GL_MIN_ABS) return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_min_abs);
693   else if (type == GL_EXISTS)  return((PetscErrorCode (*)(void*, void*, PetscInt, ...))&PCTFS_rvec_exists);
694 
695   /* catch all ... not good if we get here */
696   return(NULL);
697 }
698 
699