xref: /petsc/src/ksp/pc/impls/tfs/ivec.c (revision 1d7d09051fe52a3a6e6d1decb2e711a4db60a9c8) 
1 #define PETSCKSP_DLL
2 
3 /**********************************ivec.c**************************************
4 
5 Author: Henry M. Tufo III
6 
7 e-mail: hmt@cs.brown.edu
8 
9 snail-mail:
10 Division of Applied Mathematics
11 Brown University
12 Providence, RI 02912
13 
14 Last Modification:
15 6.21.97
16 ***********************************ivec.c*************************************/
17 
18 /**********************************ivec.c**************************************
19 File Description:
20 -----------------
21 
22 ***********************************ivec.c*************************************/
23 #include "src/ksp/pc/impls/tfs/tfs.h"
24 
25 /* sorting args ivec.c ivec.c ... */
26 #define   SORT_OPT	6
27 #define   SORT_STACK	50000
28 
29 
30 /* allocate an address and size stack for sorter(s) */
31 static void *offset_stack[2*SORT_STACK];
32 static int   size_stack[SORT_STACK];
33 
34 /**********************************ivec.c**************************************
35 Function ivec_copy()
36 
37 Input :
38 Output:
39 Return:
40 Description:
41 ***********************************ivec.c*************************************/
42 int *ivec_copy( int *arg1,  int *arg2,  int n)
43 {
44   while (n--)  {*arg1++ = *arg2++;}
45   return(arg1);
46 }
47 
48 
49 
50 /**********************************ivec.c**************************************
51 Function ivec_zero()
52 
53 Input :
54 Output:
55 Return:
56 Description:
57 ***********************************ivec.c*************************************/
58 PetscErrorCode ivec_zero( int *arg1,  int n)
59 {
60   PetscFunctionBegin;
61   while (n--)  {*arg1++ = 0;}
62   PetscFunctionReturn(0);
63 }
64 
65 /**********************************ivec.c**************************************
66 Function ivec_set()
67 
68 Input :
69 Output:
70 Return:
71 Description:
72 ***********************************ivec.c*************************************/
73 PetscErrorCode ivec_set( int *arg1,  int arg2,  int n)
74 {
75   PetscFunctionBegin;
76   while (n--)  {*arg1++ = arg2;}
77   PetscFunctionReturn(0);
78 }
79 
80 /**********************************ivec.c**************************************
81 Function ivec_max()
82 
83 Input :
84 Output:
85 Return:
86 Description:
87 ***********************************ivec.c*************************************/
88 PetscErrorCode ivec_max( int *arg1,  int *arg2,  int n)
89 {
90   PetscFunctionBegin;
91   while (n--)  {*arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++;}
92   PetscFunctionReturn(0);
93 }
94 
95 
96 
97 /**********************************ivec.c**************************************
98 Function ivec_min()
99 
100 Input :
101 Output:
102 Return:
103 Description:
104 ***********************************ivec.c*************************************/
105 PetscErrorCode ivec_min( int *arg1,  int *arg2,  int n)
106 {
107   PetscFunctionBegin;
108   while (n--)  {*(arg1) = PetscMin(*arg1,*arg2); arg1++; arg2++;}
109   PetscFunctionReturn(0);
110 }
111 
112 
113 
114 /**********************************ivec.c**************************************
115 Function ivec_mult()
116 
117 Input :
118 Output:
119 Return:
120 Description:
121 ***********************************ivec.c*************************************/
122 PetscErrorCode ivec_mult( int *arg1,  int *arg2,  int n)
123 {
124   PetscFunctionBegin;
125   while (n--)  {*arg1++ *= *arg2++;}
126   PetscFunctionReturn(0);
127 }
128 
129 
130 
131 /**********************************ivec.c**************************************
132 Function ivec_add()
133 
134 Input :
135 Output:
136 Return:
137 Description:
138 ***********************************ivec.c*************************************/
139 PetscErrorCode ivec_add( int *arg1,  int *arg2,  int n)
140 {
141   PetscFunctionBegin;
142   while (n--)  {*arg1++ += *arg2++;}
143   PetscFunctionReturn(0);
144 }
145 
146 
147 
148 /**********************************ivec.c**************************************
149 Function ivec_lxor()
150 
151 Input :
152 Output:
153 Return:
154 Description:
155 ***********************************ivec.c*************************************/
156 PetscErrorCode ivec_lxor( int *arg1,  int *arg2,  int n)
157 {
158   PetscFunctionBegin;
159   while (n--) {*arg1=((*arg1 || *arg2) && !(*arg1 && *arg2)); arg1++; arg2++;}
160   PetscFunctionReturn(0);
161 }
162 
163 
164 
165 /**********************************ivec.c**************************************
166 Function ivec_xor()
167 
168 Input :
169 Output:
170 Return:
171 Description:
172 ***********************************ivec.c*************************************/
173 PetscErrorCode ivec_xor( int *arg1,  int *arg2,  int n)
174 {
175   PetscFunctionBegin;
176   while (n--)  {*arg1++ ^= *arg2++;}
177   PetscFunctionReturn(0);
178 }
179 
180 
181 
182 /**********************************ivec.c**************************************
183 Function ivec_or()
184 
185 Input :
186 Output:
187 Return:
188 Description:
189 ***********************************ivec.c*************************************/
190 PetscErrorCode ivec_or( int *arg1,  int *arg2,  int n)
191 {
192   PetscFunctionBegin;
193   while (n--)  {*arg1++ |= *arg2++;}
194   PetscFunctionReturn(0);
195 }
196 
197 
198 
199 /**********************************ivec.c**************************************
200 Function ivec_lor()
201 
202 Input :
203 Output:
204 Return:
205 Description:
206 ***********************************ivec.c*************************************/
207 PetscErrorCode ivec_lor( int *arg1,  int *arg2,  int n)
208 {
209   PetscFunctionBegin;
210   while (n--)  {*arg1 = (*arg1 || *arg2); arg1++; arg2++;}
211   PetscFunctionReturn(0);
212 }
213 
214 /**********************************ivec.c**************************************
215 Function ivec_and()
216 
217 Input :
218 Output:
219 Return:
220 Description:
221 ***********************************ivec.c*************************************/
222 PetscErrorCode ivec_and( int *arg1,  int *arg2,  int n)
223 {
224   PetscFunctionBegin;
225   while (n--)  {*arg1++ &= *arg2++;}
226   PetscFunctionReturn(0);
227 }
228 
229 
230 
231 /**********************************ivec.c**************************************
232 Function ivec_land()
233 
234 Input :
235 Output:
236 Return:
237 Description:
238 ***********************************ivec.c*************************************/
239 PetscErrorCode ivec_land( int *arg1,  int *arg2,  int n)
240 {
241   PetscFunctionBegin;
242   while (n--) {*arg1 = (*arg1 && *arg2); arg1++; arg2++;}
243   PetscFunctionReturn(0);
244 }
245 
246 
247 
248 /**********************************ivec.c**************************************
249 Function ivec_and3()
250 
251 Input :
252 Output:
253 Return:
254 Description:
255 ***********************************ivec.c*************************************/
256 PetscErrorCode ivec_and3( int *arg1,  int *arg2,  int *arg3, int n)
257 {
258   PetscFunctionBegin;
259   while (n--)  {*arg1++ = (*arg2++ & *arg3++);}
260   PetscFunctionReturn(0);
261 }
262 
263 
264 
265 /**********************************ivec.c**************************************
266 Function ivec_sum
267 
268 Input :
269 Output:
270 Return:
271 Description:
272 ***********************************ivec.c*************************************/
273 int ivec_sum( int *arg1,  int n)
274 {
275    int tmp = 0;
276 
277 
278   while (n--) {tmp += *arg1++;}
279   return(tmp);
280 }
281 
282 /**********************************ivec.c**************************************
283 Function ivec_non_uniform()
284 
285 Input :
286 Output:
287 Return:
288 Description:
289 ***********************************ivec.c*************************************/
290 PetscErrorCode ivec_non_uniform(int *arg1, int *arg2,  int n,  int *arg3)
291 {
292    int i, j, type;
293 
294 
295   PetscFunctionBegin;
296   /* LATER: if we're really motivated we can sort and then unsort */
297   for (i=0;i<n;)
298     {
299       /* clump 'em for now */
300       j=i+1;
301       type = arg3[i];
302       while ((j<n)&&(arg3[j]==type))
303 	{j++;}
304 
305       /* how many together */
306       j -= i;
307 
308       /* call appropriate ivec function */
309       if (type == GL_MAX)
310 	{ivec_max(arg1,arg2,j);}
311       else if (type == GL_MIN)
312 	{ivec_min(arg1,arg2,j);}
313       else if (type == GL_MULT)
314 	{ivec_mult(arg1,arg2,j);}
315       else if (type == GL_ADD)
316 	{ivec_add(arg1,arg2,j);}
317       else if (type == GL_B_XOR)
318 	{ivec_xor(arg1,arg2,j);}
319       else if (type == GL_B_OR)
320 	{ivec_or(arg1,arg2,j);}
321       else if (type == GL_B_AND)
322 	{ivec_and(arg1,arg2,j);}
323       else if (type == GL_L_XOR)
324 	{ivec_lxor(arg1,arg2,j);}
325       else if (type == GL_L_OR)
326 	{ivec_lor(arg1,arg2,j);}
327       else if (type == GL_L_AND)
328 	{ivec_land(arg1,arg2,j);}
329       else
330 	{error_msg_fatal("unrecognized type passed to ivec_non_uniform()!");}
331 
332       arg1+=j; arg2+=j; i+=j;
333     }
334   PetscFunctionReturn(0);
335 }
336 
337 
338 
339 /**********************************ivec.c**************************************
340 Function ivec_addr()
341 
342 Input :
343 Output:
344 Return:
345 Description:
346 ***********************************ivec.c*************************************/
347 vfp ivec_fct_addr( int type)
348 {
349   PetscFunctionBegin;
350   if (type == NON_UNIFORM)
351     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_non_uniform);}
352   else if (type == GL_MAX)
353     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_max);}
354   else if (type == GL_MIN)
355     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_min);}
356   else if (type == GL_MULT)
357     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_mult);}
358   else if (type == GL_ADD)
359     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_add);}
360   else if (type == GL_B_XOR)
361     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_xor);}
362   else if (type == GL_B_OR)
363     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_or);}
364   else if (type == GL_B_AND)
365     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_and);}
366   else if (type == GL_L_XOR)
367     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_lxor);}
368   else if (type == GL_L_OR)
369     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_lor);}
370   else if (type == GL_L_AND)
371     {return((PetscErrorCode (*)(void*, void *, int, ...))&ivec_land);}
372 
373   /* catch all ... not good if we get here */
374   return(NULL);
375 }
376 
377 
378 /******************************************************************************
379 Function: ivec_sort().
380 
381 Input : offset of list to be sorted, number of elements to be sorted.
382 Output: sorted list (in ascending order).
383 Return: none.
384 Description: stack based (nonrecursive) quicksort w/brute-shell bottom.
385 ******************************************************************************/
386 PetscErrorCode ivec_sort( int *ar,  int size)
387 {
388    int *pi, *pj, temp;
389    int **top_a = (int **) offset_stack;
390    int *top_s = size_stack, *bottom_s = size_stack;
391 
392 
393   /* we're really interested in the offset of the last element */
394   /* ==> length of the list is now size + 1                    */
395   size--;
396 
397   /* do until we're done ... return when stack is exhausted */
398   for (;;)
399     {
400       /* if list is large enough use quicksort partition exchange code */
401       if (size > SORT_OPT)
402 	{
403 	  /* start up pointer at element 1 and down at size     */
404 	  pi = ar+1;
405 	  pj = ar+size;
406 
407 	  /* find middle element in list and swap w/ element 1 */
408 	  SWAP(*(ar+(size>>1)),*pi)
409 
410 	  /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
411 	  /* note ==> pivot_value in index 0                   */
412 	  if (*pi > *pj)
413 	    {SWAP(*pi,*pj)}
414 	  if (*ar > *pj)
415 	    {SWAP(*ar,*pj)}
416 	  else if (*pi > *ar)
417 	    {SWAP(*(ar),*(ar+1))}
418 
419 	  /* partition about pivot_value ...  	                    */
420 	  /* note lists of length 2 are not guaranteed to be sorted */
421 	  for(;;)
422 	    {
423 	      /* walk up ... and down ... swap if equal to pivot! */
424 	      do pi++; while (*pi<*ar);
425 	      do pj--; while (*pj>*ar);
426 
427 	      /* if we've crossed we're done */
428 	      if (pj<pi) break;
429 
430 	      /* else swap */
431 	      SWAP(*pi,*pj)
432 	    }
433 
434 	  /* place pivot_value in it's correct location */
435 	  SWAP(*ar,*pj)
436 
437 	  /* test stack_size to see if we've exhausted our stack */
438 	  if (top_s-bottom_s >= SORT_STACK)
439 	    {error_msg_fatal("ivec_sort() :: STACK EXHAUSTED!!!");}
440 
441 	  /* push right hand child iff length > 1 */
442 	  if ((*top_s = size-((int) (pi-ar))))
443 	    {
444 	      *(top_a++) = pi;
445 	      size -= *top_s+2;
446 	      top_s++;
447 	    }
448 	  /* set up for next loop iff there is something to do */
449 	  else if (size -= *top_s+2)
450 	    {;}
451 	  /* might as well pop - note NR_OPT >=2 ==> we're ok! */
452 	  else
453 	    {
454 	      ar = *(--top_a);
455 	      size = *(--top_s);
456 	    }
457 	}
458 
459       /* else sort small list directly then pop another off stack */
460       else
461 	{
462 	  /* insertion sort for bottom */
463           for (pj=ar+1;pj<=ar+size;pj++)
464             {
465               temp = *pj;
466               for (pi=pj-1;pi>=ar;pi--)
467                 {
468                   if (*pi <= temp) break;
469                   *(pi+1)=*pi;
470                 }
471               *(pi+1)=temp;
472 	    }
473 
474 	  /* check to see if stack is exhausted ==> DONE */
475 	  if (top_s==bottom_s)   PetscFunctionReturn(0);
476 
477 	  /* else pop another list from the stack */
478 	  ar = *(--top_a);
479 	  size = *(--top_s);
480 	}
481     }
482   PetscFunctionReturn(0);
483 }
484 
485 
486 
487 /******************************************************************************
488 Function: ivec_sort_companion().
489 
490 Input : offset of list to be sorted, number of elements to be sorted.
491 Output: sorted list (in ascending order).
492 Return: none.
493 Description: stack based (nonrecursive) quicksort w/brute-shell bottom.
494 ******************************************************************************/
495 PetscErrorCode ivec_sort_companion( int *ar,  int *ar2,  int size)
496 {
497    int *pi, *pj, temp, temp2;
498    int **top_a = (int **)offset_stack;
499    int *top_s = size_stack, *bottom_s = size_stack;
500    int *pi2, *pj2;
501    int mid;
502 
503   PetscFunctionBegin;
504   /* we're really interested in the offset of the last element */
505   /* ==> length of the list is now size + 1                    */
506   size--;
507 
508   /* do until we're done ... return when stack is exhausted */
509   for (;;)
510     {
511       /* if list is large enough use quicksort partition exchange code */
512       if (size > SORT_OPT)
513 	{
514 	  /* start up pointer at element 1 and down at size     */
515 	  mid = size>>1;
516 	  pi = ar+1;
517 	  pj = ar+mid;
518 	  pi2 = ar2+1;
519 	  pj2 = ar2+mid;
520 
521 	  /* find middle element in list and swap w/ element 1 */
522 	  SWAP(*pi,*pj)
523 	  SWAP(*pi2,*pj2)
524 
525 	  /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
526 	  /* note ==> pivot_value in index 0                   */
527 	  pj = ar+size;
528 	  pj2 = ar2+size;
529 	  if (*pi > *pj)
530 	    {SWAP(*pi,*pj) SWAP(*pi2,*pj2)}
531 	  if (*ar > *pj)
532 	    {SWAP(*ar,*pj) SWAP(*ar2,*pj2)}
533 	  else if (*pi > *ar)
534 	    {SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1))}
535 
536 	  /* partition about pivot_value ...  	                    */
537 	  /* note lists of length 2 are not guaranteed to be sorted */
538 	  for(;;)
539 	    {
540 	      /* walk up ... and down ... swap if equal to pivot! */
541 	      do {pi++; pi2++;} while (*pi<*ar);
542 	      do {pj--; pj2--;} while (*pj>*ar);
543 
544 	      /* if we've crossed we're done */
545 	      if (pj<pi) break;
546 
547 	      /* else swap */
548 	      SWAP(*pi,*pj)
549 	      SWAP(*pi2,*pj2)
550 	    }
551 
552 	  /* place pivot_value in it's correct location */
553 	  SWAP(*ar,*pj)
554 	  SWAP(*ar2,*pj2)
555 
556 	  /* test stack_size to see if we've exhausted our stack */
557 	  if (top_s-bottom_s >= SORT_STACK)
558 	    {error_msg_fatal("ivec_sort_companion() :: STACK EXHAUSTED!!!");}
559 
560 	  /* push right hand child iff length > 1 */
561 	  if ((*top_s = size-((int) (pi-ar))))
562 	    {
563 	      *(top_a++) = pi;
564 	      *(top_a++) = pi2;
565 	      size -= *top_s+2;
566 	      top_s++;
567 	    }
568 	  /* set up for next loop iff there is something to do */
569 	  else if (size -= *top_s+2)
570 	    {;}
571 	  /* might as well pop - note NR_OPT >=2 ==> we're ok! */
572 	  else
573 	    {
574 	      ar2 = *(--top_a);
575 	      ar  = *(--top_a);
576 	      size = *(--top_s);
577 	    }
578 	}
579 
580       /* else sort small list directly then pop another off stack */
581       else
582 	{
583 	  /* insertion sort for bottom */
584           for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++)
585             {
586               temp = *pj;
587               temp2 = *pj2;
588               for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--)
589                 {
590                   if (*pi <= temp) break;
591                   *(pi+1)=*pi;
592                   *(pi2+1)=*pi2;
593                 }
594               *(pi+1)=temp;
595               *(pi2+1)=temp2;
596 	    }
597 
598 	  /* check to see if stack is exhausted ==> DONE */
599 	  if (top_s==bottom_s)   PetscFunctionReturn(0);
600 
601 	  /* else pop another list from the stack */
602 	  ar2 = *(--top_a);
603 	  ar  = *(--top_a);
604 	  size = *(--top_s);
605 	}
606     }
607   PetscFunctionReturn(0);
608 }
609 
610 
611 
612 /******************************************************************************
613 Function: ivec_sort_companion_hack().
614 
615 Input : offset of list to be sorted, number of elements to be sorted.
616 Output: sorted list (in ascending order).
617 Return: none.
618 Description: stack based (nonrecursive) quicksort w/brute-shell bottom.
619 ******************************************************************************/
620 PetscErrorCode ivec_sort_companion_hack( int *ar,  int **ar2, int size)
621 {
622    int *pi, *pj, temp, *ptr;
623    int **top_a = (int **)offset_stack;
624    int *top_s = size_stack, *bottom_s = size_stack;
625    int **pi2, **pj2;
626    int mid;
627 
628   PetscFunctionBegin;
629   /* we're really interested in the offset of the last element */
630   /* ==> length of the list is now size + 1                    */
631   size--;
632 
633   /* do until we're done ... return when stack is exhausted */
634   for (;;)
635     {
636       /* if list is large enough use quicksort partition exchange code */
637       if (size > SORT_OPT)
638 	{
639 	  /* start up pointer at element 1 and down at size     */
640 	  mid = size>>1;
641 	  pi = ar+1;
642 	  pj = ar+mid;
643 	  pi2 = ar2+1;
644 	  pj2 = ar2+mid;
645 
646 	  /* find middle element in list and swap w/ element 1 */
647 	  SWAP(*pi,*pj)
648 	  P_SWAP(*pi2,*pj2)
649 
650 	  /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */
651 	  /* note ==> pivot_value in index 0                   */
652 	  pj = ar+size;
653 	  pj2 = ar2+size;
654 	  if (*pi > *pj)
655 	    {SWAP(*pi,*pj) P_SWAP(*pi2,*pj2)}
656 	  if (*ar > *pj)
657 	    {SWAP(*ar,*pj) P_SWAP(*ar2,*pj2)}
658 	  else if (*pi > *ar)
659 	    {SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1))}
660 
661 	  /* partition about pivot_value ...  	                    */
662 	  /* note lists of length 2 are not guaranteed to be sorted */
663 	  for(;;)
664 	    {
665 	      /* walk up ... and down ... swap if equal to pivot! */
666 	      do {pi++; pi2++;} while (*pi<*ar);
667 	      do {pj--; pj2--;} while (*pj>*ar);
668 
669 	      /* if we've crossed we're done */
670 	      if (pj<pi) break;
671 
672 	      /* else swap */
673 	      SWAP(*pi,*pj)
674 	      P_SWAP(*pi2,*pj2)
675 	    }
676 
677 	  /* place pivot_value in it's correct location */
678 	  SWAP(*ar,*pj)
679 	  P_SWAP(*ar2,*pj2)
680 
681 	  /* test stack_size to see if we've exhausted our stack */
682 	  if (top_s-bottom_s >= SORT_STACK)
683          {error_msg_fatal("ivec_sort_companion_hack() :: STACK EXHAUSTED!!!");}
684 
685 	  /* push right hand child iff length > 1 */
686 	  if ((*top_s = size-((int) (pi-ar))))
687 	    {
688 	      *(top_a++) = pi;
689 	      *(top_a++) = (int*) pi2;
690 	      size -= *top_s+2;
691 	      top_s++;
692 	    }
693 	  /* set up for next loop iff there is something to do */
694 	  else if (size -= *top_s+2)
695 	    {;}
696 	  /* might as well pop - note NR_OPT >=2 ==> we're ok! */
697 	  else
698 	    {
699 	      ar2 = (int **) *(--top_a);
700 	      ar  = *(--top_a);
701 	      size = *(--top_s);
702 	    }
703 	}
704 
705       /* else sort small list directly then pop another off stack */
706       else
707 	{
708 	  /* insertion sort for bottom */
709           for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++)
710             {
711               temp = *pj;
712               ptr = *pj2;
713               for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--)
714                 {
715                   if (*pi <= temp) break;
716                   *(pi+1)=*pi;
717                   *(pi2+1)=*pi2;
718                 }
719               *(pi+1)=temp;
720               *(pi2+1)=ptr;
721 	    }
722 
723 	  /* check to see if stack is exhausted ==> DONE */
724 	  if (top_s==bottom_s)   PetscFunctionReturn(0);
725 
726 	  /* else pop another list from the stack */
727 	  ar2 = (int **)*(--top_a);
728 	  ar  = *(--top_a);
729 	  size = *(--top_s);
730 	}
731     }
732   PetscFunctionReturn(0);
733 }
734 
735 
736 
737 /******************************************************************************
738 Function: SMI_sort().
739 Input : offset of list to be sorted, number of elements to be sorted.
740 Output: sorted list (in ascending order).
741 Return: none.
742 Description: stack based (nonrecursive) quicksort w/brute-shell bottom.
743 ******************************************************************************/
744 PetscErrorCode SMI_sort(void *ar1, void *ar2, int size, int type)
745 {
746   PetscFunctionBegin;
747   if (type == SORT_INTEGER)
748     {
749       if (ar2)
750 	{ivec_sort_companion((int*)ar1,(int*)ar2,size);}
751       else
752 	{ivec_sort((int*)ar1,size);}
753     }
754   else if (type == SORT_INT_PTR)
755     {
756       if (ar2)
757 	{ivec_sort_companion_hack((int*)ar1,(int **)ar2,size);}
758       else
759 	{ivec_sort((int*)ar1,size);}
760     }
761 
762   else
763     {
764       error_msg_fatal("SMI_sort only does SORT_INTEGER!");
765     }
766   PetscFunctionReturn(0);
767 }
768 
769 
770 
771 /**********************************ivec.c**************************************
772 Function ivec_linear_search()
773 
774 Input :
775 Output:
776 Return:
777 Description:
778 ***********************************ivec.c*************************************/
779 int
780 ivec_linear_search( int item,  int *list,  int n)
781 {
782    int tmp = n-1;
783   PetscFunctionBegin;
784   while (n--)  {if (*list++ == item) {return(tmp-n);}}
785   return(-1);
786 }
787 
788 
789 
790 /**********************************ivec.c**************************************
791 Function ivec_binary_search()
792 
793 Input :
794 Output:
795 Return:
796 Description:
797 ***********************************ivec.c*************************************/
798 int
799 ivec_binary_search( int item,  int *list,  int rh)
800 {
801    int mid, lh=0;
802 
803   rh--;
804   while (lh<=rh)
805     {
806       mid = (lh+rh)>>1;
807       if (*(list+mid) == item)
808 	{return(mid);}
809       if (*(list+mid) > item)
810 	{rh = mid-1;}
811       else
812 	{lh = mid+1;}
813     }
814   return(-1);
815 }
816 
817 
818 /********************************ivec.c**************************************
819 Function rvec_copy()
820 
821 Input :
822 Output:
823 Return:
824 Description:
825 *********************************ivec.c*************************************/
826 PetscErrorCode rvec_copy( PetscScalar *arg1,  PetscScalar *arg2,  int n)
827 {
828   PetscFunctionBegin;
829   while (n--)  {*arg1++ = *arg2++;}
830   PetscFunctionReturn(0);
831 }
832 
833 
834 
835 /********************************ivec.c**************************************
836 Function rvec_zero()
837 
838 Input :
839 Output:
840 Return:
841 Description:
842 *********************************ivec.c*************************************/
843 PetscErrorCode rvec_zero( PetscScalar *arg1,  int n)
844 {
845   PetscFunctionBegin;
846   while (n--)  {*arg1++ = 0.0;}
847   PetscFunctionReturn(0);
848 }
849 
850 
851 
852 /**********************************ivec.c**************************************
853 Function rvec_one()
854 
855 Input :
856 Output:
857 Return:
858 Description:
859 ***********************************ivec.c*************************************/
860 PetscErrorCode rvec_one( PetscScalar *arg1,  int n)
861 {
862   PetscFunctionBegin;
863   while (n--)  {*arg1++ = 1.0;}
864   PetscFunctionReturn(0);
865 }
866 
867 /**********************************ivec.c**************************************
868 Function rvec_set()
869 
870 Input :
871 Output:
872 Return:
873 Description:
874 ***********************************ivec.c*************************************/
875 PetscErrorCode rvec_set( PetscScalar *arg1,  PetscScalar arg2,  int n)
876 {
877   PetscFunctionBegin;
878   while (n--)  {*arg1++ = arg2;}
879   PetscFunctionReturn(0);
880 }
881 
882 
883 
884 /**********************************ivec.c**************************************
885 Function rvec_scale()
886 
887 Input :
888 Output:
889 Return:
890 Description:
891 ***********************************ivec.c*************************************/
892 PetscErrorCode rvec_scale( PetscScalar *arg1,  PetscScalar arg2,  int n)
893 {
894   PetscFunctionBegin;
895   while (n--)  {*arg1++ *= arg2;}
896   PetscFunctionReturn(0);
897 }
898 
899 
900 
901 /********************************ivec.c**************************************
902 Function rvec_add()
903 
904 Input :
905 Output:
906 Return:
907 Description:
908 *********************************ivec.c*************************************/
909 PetscErrorCode rvec_add( PetscScalar *arg1,  PetscScalar *arg2,  int n)
910 {
911   PetscFunctionBegin;
912   while (n--)  {*arg1++ += *arg2++;}
913   PetscFunctionReturn(0);
914 }
915 
916 /********************************ivec.c**************************************
917 Function rvec_mult()
918 
919 Input :
920 Output:
921 Return:
922 Description:
923 *********************************ivec.c*************************************/
924 PetscErrorCode rvec_mult( PetscScalar *arg1,  PetscScalar *arg2,  int n)
925 {
926   PetscFunctionBegin;
927   while (n--)  {*arg1++ *= *arg2++;}
928   PetscFunctionReturn(0);
929 }
930 
931 
932 
933 /********************************ivec.c**************************************
934 Function rvec_max()
935 
936 Input :
937 Output:
938 Return:
939 Description:
940 *********************************ivec.c*************************************/
941 PetscErrorCode rvec_max( PetscScalar *arg1,  PetscScalar *arg2,  int n)
942 {
943   PetscFunctionBegin;
944   while (n--)  {*arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++;}
945   PetscFunctionReturn(0);
946 }
947 
948 
949 
950 /********************************ivec.c**************************************
951 Function rvec_max_abs()
952 
953 Input :
954 Output:
955 Return:
956 Description:
957 *********************************ivec.c*************************************/
958 PetscErrorCode rvec_max_abs( PetscScalar *arg1,  PetscScalar *arg2,  int n)
959 {
960   PetscFunctionBegin;
961   while (n--)  {*arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++;}
962   PetscFunctionReturn(0);
963 }
964 
965 
966 
967 /********************************ivec.c**************************************
968 Function rvec_min()
969 
970 Input :
971 Output:
972 Return:
973 Description:
974 *********************************ivec.c*************************************/
975 PetscErrorCode rvec_min( PetscScalar *arg1,  PetscScalar *arg2,  int n)
976 {
977   PetscFunctionBegin;
978   while (n--)  {*arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++;}
979   PetscFunctionReturn(0);
980 }
981 
982 
983 
984 /********************************ivec.c**************************************
985 Function rvec_min_abs()
986 
987 Input :
988 Output:
989 Return:
990 Description:
991 *********************************ivec.c*************************************/
992 PetscErrorCode rvec_min_abs( PetscScalar *arg1,  PetscScalar *arg2,  int n)
993 {
994   PetscFunctionBegin;
995   while (n--)  {*arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++;}
996   PetscFunctionReturn(0);
997 }
998 
999 
1000 
1001 /********************************ivec.c**************************************
1002 Function rvec_exists()
1003 
1004 Input :
1005 Output:
1006 Return:
1007 Description:
1008 *********************************ivec.c*************************************/
1009 PetscErrorCode rvec_exists( PetscScalar *arg1,  PetscScalar *arg2,  int n)
1010 {
1011   PetscFunctionBegin;
1012   while (n--)  {*arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++;}
1013   PetscFunctionReturn(0);
1014 }
1015 
1016 
1017 
1018 /**********************************ivec.c**************************************
1019 Function rvec_non_uniform()
1020 
1021 Input :
1022 Output:
1023 Return:
1024 Description:
1025 ***********************************ivec.c*************************************/
1026 PetscErrorCode rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2,  int n,  int *arg3)
1027 {
1028    int i, j, type;
1029 
1030   PetscFunctionBegin;
1031   /* LATER: if we're really motivated we can sort and then unsort */
1032   for (i=0;i<n;)
1033     {
1034       /* clump 'em for now */
1035       j=i+1;
1036       type = arg3[i];
1037       while ((j<n)&&(arg3[j]==type))
1038 	{j++;}
1039 
1040       /* how many together */
1041       j -= i;
1042 
1043       /* call appropriate ivec function */
1044       if (type == GL_MAX)
1045 	{rvec_max(arg1,arg2,j);}
1046       else if (type == GL_MIN)
1047 	{rvec_min(arg1,arg2,j);}
1048       else if (type == GL_MULT)
1049 	{rvec_mult(arg1,arg2,j);}
1050       else if (type == GL_ADD)
1051 	{rvec_add(arg1,arg2,j);}
1052       else if (type == GL_MAX_ABS)
1053 	{rvec_max_abs(arg1,arg2,j);}
1054       else if (type == GL_MIN_ABS)
1055 	{rvec_min_abs(arg1,arg2,j);}
1056       else if (type == GL_EXISTS)
1057 	{rvec_exists(arg1,arg2,j);}
1058       else
1059 	{error_msg_fatal("unrecognized type passed to rvec_non_uniform()!");}
1060 
1061       arg1+=j; arg2+=j; i+=j;
1062     }
1063   PetscFunctionReturn(0);
1064 }
1065 
1066 
1067 
1068 /**********************************ivec.c**************************************
1069 Function rvec_fct_addr()
1070 
1071 Input :
1072 Output:
1073 Return:
1074 Description:
1075 ***********************************ivec.c*************************************/
1076 vfp rvec_fct_addr( int type)
1077 {
1078   if (type == NON_UNIFORM)
1079     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_non_uniform);}
1080   else if (type == GL_MAX)
1081     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_max);}
1082   else if (type == GL_MIN)
1083     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_min);}
1084   else if (type == GL_MULT)
1085     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_mult);}
1086   else if (type == GL_ADD)
1087     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_add);}
1088   else if (type == GL_MAX_ABS)
1089     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_max_abs);}
1090   else if (type == GL_MIN_ABS)
1091     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_min_abs);}
1092   else if (type == GL_EXISTS)
1093     {return((PetscErrorCode (*)(void*, void *, int, ...))&rvec_exists);}
1094 
1095   /* catch all ... not good if we get here */
1096   return(NULL);
1097 }
1098 
1099 
1100 
1101 
1102 
1103 
1104