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