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