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