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_COMM_SELF,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) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ivec_sort() :: STACK EXHAUSTED!!!"); 283 284 /* push right hand child iff length > 1 */ 285 if ((*top_s = size-((PetscInt) (pi-ar)))) 286 { 287 *(top_a++) = pi; 288 size -= *top_s+2; 289 top_s++; 290 } 291 /* set up for next loop iff there is something to do */ 292 else if (size -= *top_s+2) 293 {;} 294 /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 295 else 296 { 297 ar = *(--top_a); 298 size = *(--top_s); 299 } 300 } 301 302 /* else sort small list directly then pop another off stack */ 303 else 304 { 305 /* insertion sort for bottom */ 306 for (pj=ar+1;pj<=ar+size;pj++) 307 { 308 temp = *pj; 309 for (pi=pj-1;pi>=ar;pi--) 310 { 311 if (*pi <= temp) break; 312 *(pi+1)=*pi; 313 } 314 *(pi+1)=temp; 315 } 316 317 /* check to see if stack is exhausted ==> DONE */ 318 if (top_s==bottom_s) PetscFunctionReturn(0); 319 320 /* else pop another list from the stack */ 321 ar = *(--top_a); 322 size = *(--top_s); 323 } 324 } 325 PetscFunctionReturn(0); 326 } 327 328 /******************************************************************************/ 329 PetscErrorCode ivec_sort_companion( PetscInt *ar, PetscInt *ar2, PetscInt size) 330 { 331 PetscInt *pi, *pj, temp, temp2; 332 PetscInt **top_a = (PetscInt **)offset_stack; 333 PetscInt *top_s = size_stack, *bottom_s = size_stack; 334 PetscInt *pi2, *pj2; 335 PetscInt mid; 336 337 PetscFunctionBegin; 338 /* we're really interested in the offset of the last element */ 339 /* ==> length of the list is now size + 1 */ 340 size--; 341 342 /* do until we're done ... return when stack is exhausted */ 343 for (;;) 344 { 345 /* if list is large enough use quicksort partition exchange code */ 346 if (size > SORT_OPT) 347 { 348 /* start up pointer at element 1 and down at size */ 349 mid = size>>1; 350 pi = ar+1; 351 pj = ar+mid; 352 pi2 = ar2+1; 353 pj2 = ar2+mid; 354 355 /* find middle element in list and swap w/ element 1 */ 356 SWAP(*pi,*pj) 357 SWAP(*pi2,*pj2) 358 359 /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 360 /* note ==> pivot_value in index 0 */ 361 pj = ar+size; 362 pj2 = ar2+size; 363 if (*pi > *pj) 364 {SWAP(*pi,*pj) SWAP(*pi2,*pj2)} 365 if (*ar > *pj) 366 {SWAP(*ar,*pj) SWAP(*ar2,*pj2)} 367 else if (*pi > *ar) 368 {SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1))} 369 370 /* partition about pivot_value ... */ 371 /* note lists of length 2 are not guaranteed to be sorted */ 372 for(;;) 373 { 374 /* walk up ... and down ... swap if equal to pivot! */ 375 do {pi++; pi2++;} while (*pi<*ar); 376 do {pj--; pj2--;} while (*pj>*ar); 377 378 /* if we've crossed we're done */ 379 if (pj<pi) break; 380 381 /* else swap */ 382 SWAP(*pi,*pj) 383 SWAP(*pi2,*pj2) 384 } 385 386 /* place pivot_value in it's correct location */ 387 SWAP(*ar,*pj) 388 SWAP(*ar2,*pj2) 389 390 /* test stack_size to see if we've exhausted our stack */ 391 if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ivec_sort_companion() :: STACK EXHAUSTED!!!"); 392 393 /* push right hand child iff length > 1 */ 394 if ((*top_s = size-((PetscInt) (pi-ar)))) 395 { 396 *(top_a++) = pi; 397 *(top_a++) = pi2; 398 size -= *top_s+2; 399 top_s++; 400 } 401 /* set up for next loop iff there is something to do */ 402 else if (size -= *top_s+2) 403 {;} 404 /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 405 else 406 { 407 ar2 = *(--top_a); 408 ar = *(--top_a); 409 size = *(--top_s); 410 } 411 } 412 413 /* else sort small list directly then pop another off stack */ 414 else 415 { 416 /* insertion sort for bottom */ 417 for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) 418 { 419 temp = *pj; 420 temp2 = *pj2; 421 for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) 422 { 423 if (*pi <= temp) break; 424 *(pi+1)=*pi; 425 *(pi2+1)=*pi2; 426 } 427 *(pi+1)=temp; 428 *(pi2+1)=temp2; 429 } 430 431 /* check to see if stack is exhausted ==> DONE */ 432 if (top_s==bottom_s) PetscFunctionReturn(0); 433 434 /* else pop another list from the stack */ 435 ar2 = *(--top_a); 436 ar = *(--top_a); 437 size = *(--top_s); 438 } 439 } 440 PetscFunctionReturn(0); 441 } 442 443 /******************************************************************************/ 444 PetscErrorCode ivec_sort_companion_hack( PetscInt *ar, PetscInt **ar2, PetscInt size) 445 { 446 PetscInt *pi, *pj, temp, *ptr; 447 PetscInt **top_a = (PetscInt **)offset_stack; 448 PetscInt *top_s = size_stack, *bottom_s = size_stack; 449 PetscInt **pi2, **pj2; 450 PetscInt mid; 451 452 PetscFunctionBegin; 453 /* we're really interested in the offset of the last element */ 454 /* ==> length of the list is now size + 1 */ 455 size--; 456 457 /* do until we're done ... return when stack is exhausted */ 458 for (;;) 459 { 460 /* if list is large enough use quicksort partition exchange code */ 461 if (size > SORT_OPT) 462 { 463 /* start up pointer at element 1 and down at size */ 464 mid = size>>1; 465 pi = ar+1; 466 pj = ar+mid; 467 pi2 = ar2+1; 468 pj2 = ar2+mid; 469 470 /* find middle element in list and swap w/ element 1 */ 471 SWAP(*pi,*pj) 472 P_SWAP(*pi2,*pj2) 473 474 /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 475 /* note ==> pivot_value in index 0 */ 476 pj = ar+size; 477 pj2 = ar2+size; 478 if (*pi > *pj) 479 {SWAP(*pi,*pj) P_SWAP(*pi2,*pj2)} 480 if (*ar > *pj) 481 {SWAP(*ar,*pj) P_SWAP(*ar2,*pj2)} 482 else if (*pi > *ar) 483 {SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1))} 484 485 /* partition about pivot_value ... */ 486 /* note lists of length 2 are not guaranteed to be sorted */ 487 for(;;) 488 { 489 /* walk up ... and down ... swap if equal to pivot! */ 490 do {pi++; pi2++;} while (*pi<*ar); 491 do {pj--; pj2--;} while (*pj>*ar); 492 493 /* if we've crossed we're done */ 494 if (pj<pi) break; 495 496 /* else swap */ 497 SWAP(*pi,*pj) 498 P_SWAP(*pi2,*pj2) 499 } 500 501 /* place pivot_value in it's correct location */ 502 SWAP(*ar,*pj) 503 P_SWAP(*ar2,*pj2) 504 505 /* test stack_size to see if we've exhausted our stack */ 506 if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ivec_sort_companion_hack() :: STACK EXHAUSTED!!!"); 507 508 /* push right hand child iff length > 1 */ 509 if ((*top_s = size-((PetscInt) (pi-ar)))) 510 { 511 *(top_a++) = pi; 512 *(top_a++) = (PetscInt*) pi2; 513 size -= *top_s+2; 514 top_s++; 515 } 516 /* set up for next loop iff there is something to do */ 517 else if (size -= *top_s+2) 518 {;} 519 /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 520 else 521 { 522 ar2 = (PetscInt **) *(--top_a); 523 ar = *(--top_a); 524 size = *(--top_s); 525 } 526 } 527 528 /* else sort small list directly then pop another off stack */ 529 else 530 { 531 /* insertion sort for bottom */ 532 for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) 533 { 534 temp = *pj; 535 ptr = *pj2; 536 for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) 537 { 538 if (*pi <= temp) break; 539 *(pi+1)=*pi; 540 *(pi2+1)=*pi2; 541 } 542 *(pi+1)=temp; 543 *(pi2+1)=ptr; 544 } 545 546 /* check to see if stack is exhausted ==> DONE */ 547 if (top_s==bottom_s) PetscFunctionReturn(0); 548 549 /* else pop another list from the stack */ 550 ar2 = (PetscInt **)*(--top_a); 551 ar = *(--top_a); 552 size = *(--top_s); 553 } 554 } 555 PetscFunctionReturn(0); 556 } 557 558 /******************************************************************************/ 559 PetscErrorCode SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type) 560 { 561 PetscFunctionBegin; 562 if (type == SORT_INTEGER) { 563 if (ar2) ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size); 564 else ivec_sort((PetscInt*)ar1,size); 565 } else if (type == SORT_INT_PTR) { 566 if (ar2) ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt **)ar2,size); 567 else ivec_sort((PetscInt*)ar1,size); 568 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"SMI_sort only does SORT_INTEGER!"); 569 PetscFunctionReturn(0); 570 } 571 572 /***********************************ivec.c*************************************/ 573 PetscInt ivec_linear_search( PetscInt item, PetscInt *list, PetscInt n) 574 { 575 PetscInt tmp = n-1; 576 PetscFunctionBegin; 577 while (n--) {if (*list++ == item) {return(tmp-n);}} 578 return(-1); 579 } 580 581 /***********************************ivec.c*************************************/ 582 PetscInt ivec_binary_search( PetscInt item, PetscInt *list, PetscInt rh) 583 { 584 PetscInt mid, lh=0; 585 586 rh--; 587 while (lh<=rh) 588 { 589 mid = (lh+rh)>>1; 590 if (*(list+mid) == item) 591 {return(mid);} 592 if (*(list+mid) > item) 593 {rh = mid-1;} 594 else 595 {lh = mid+1;} 596 } 597 return(-1); 598 } 599 600 /*********************************ivec.c*************************************/ 601 PetscErrorCode rvec_copy( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 602 { 603 PetscFunctionBegin; 604 while (n--) {*arg1++ = *arg2++;} 605 PetscFunctionReturn(0); 606 } 607 608 /*********************************ivec.c*************************************/ 609 PetscErrorCode rvec_zero( PetscScalar *arg1, PetscInt n) 610 { 611 PetscFunctionBegin; 612 while (n--) {*arg1++ = 0.0;} 613 PetscFunctionReturn(0); 614 } 615 616 /***********************************ivec.c*************************************/ 617 PetscErrorCode rvec_one( PetscScalar *arg1, PetscInt n) 618 { 619 PetscFunctionBegin; 620 while (n--) {*arg1++ = 1.0;} 621 PetscFunctionReturn(0); 622 } 623 624 /***********************************ivec.c*************************************/ 625 PetscErrorCode rvec_set( PetscScalar *arg1, PetscScalar arg2, PetscInt n) 626 { 627 PetscFunctionBegin; 628 while (n--) {*arg1++ = arg2;} 629 PetscFunctionReturn(0); 630 } 631 632 /***********************************ivec.c*************************************/ 633 PetscErrorCode rvec_scale( PetscScalar *arg1, PetscScalar arg2, PetscInt n) 634 { 635 PetscFunctionBegin; 636 while (n--) {*arg1++ *= arg2;} 637 PetscFunctionReturn(0); 638 } 639 640 /*********************************ivec.c*************************************/ 641 PetscErrorCode rvec_add( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 642 { 643 PetscFunctionBegin; 644 while (n--) {*arg1++ += *arg2++;} 645 PetscFunctionReturn(0); 646 } 647 648 /*********************************ivec.c*************************************/ 649 PetscErrorCode rvec_mult( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 650 { 651 PetscFunctionBegin; 652 while (n--) {*arg1++ *= *arg2++;} 653 PetscFunctionReturn(0); 654 } 655 656 /*********************************ivec.c*************************************/ 657 PetscErrorCode rvec_max( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 658 { 659 PetscFunctionBegin; 660 while (n--) {*arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++;} 661 PetscFunctionReturn(0); 662 } 663 664 /*********************************ivec.c*************************************/ 665 PetscErrorCode rvec_max_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 666 { 667 PetscFunctionBegin; 668 while (n--) {*arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++;} 669 PetscFunctionReturn(0); 670 } 671 672 /*********************************ivec.c*************************************/ 673 PetscErrorCode rvec_min( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 674 { 675 PetscFunctionBegin; 676 while (n--) {*arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++;} 677 PetscFunctionReturn(0); 678 } 679 680 /*********************************ivec.c*************************************/ 681 PetscErrorCode rvec_min_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 682 { 683 PetscFunctionBegin; 684 while (n--) {*arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++;} 685 PetscFunctionReturn(0); 686 } 687 688 /*********************************ivec.c*************************************/ 689 PetscErrorCode rvec_exists( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 690 { 691 PetscFunctionBegin; 692 while (n--) {*arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++;} 693 PetscFunctionReturn(0); 694 } 695 696 /***********************************ivec.c*************************************/ 697 PetscErrorCode rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2, PetscInt n, PetscInt *arg3) 698 { 699 PetscInt i, j, type; 700 701 PetscFunctionBegin; 702 /* LATER: if we're really motivated we can sort and then unsort */ 703 for (i=0;i<n;) 704 { 705 /* clump 'em for now */ 706 j=i+1; 707 type = arg3[i]; 708 while ((j<n)&&(arg3[j]==type)) 709 {j++;} 710 711 /* how many together */ 712 j -= i; 713 714 /* call appropriate ivec function */ 715 if (type == GL_MAX) 716 {rvec_max(arg1,arg2,j);} 717 else if (type == GL_MIN) 718 {rvec_min(arg1,arg2,j);} 719 else if (type == GL_MULT) 720 {rvec_mult(arg1,arg2,j);} 721 else if (type == GL_ADD) 722 {rvec_add(arg1,arg2,j);} 723 else if (type == GL_MAX_ABS) 724 {rvec_max_abs(arg1,arg2,j);} 725 else if (type == GL_MIN_ABS) 726 {rvec_min_abs(arg1,arg2,j);} 727 else if (type == GL_EXISTS) 728 {rvec_exists(arg1,arg2,j);} 729 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to rvec_non_uniform()!"); 730 731 arg1+=j; arg2+=j; i+=j; 732 } 733 PetscFunctionReturn(0); 734 } 735 736 /***********************************ivec.c*************************************/ 737 vfp rvec_fct_addr( PetscInt type) 738 { 739 if (type == NON_UNIFORM) 740 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_non_uniform);} 741 else if (type == GL_MAX) 742 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_max);} 743 else if (type == GL_MIN) 744 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_min);} 745 else if (type == GL_MULT) 746 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_mult);} 747 else if (type == GL_ADD) 748 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_add);} 749 else if (type == GL_MAX_ABS) 750 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_max_abs);} 751 else if (type == GL_MIN_ABS) 752 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_min_abs);} 753 else if (type == GL_EXISTS) 754 {return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&rvec_exists);} 755 756 /* catch all ... not good if we get here */ 757 return(NULL); 758 } 759 760 761 762 763 764 765