1827bd09bSSatish Balay 2827bd09bSSatish Balay /**********************************ivec.c************************************** 3827bd09bSSatish Balay 4827bd09bSSatish Balay Author: Henry M. Tufo III 5827bd09bSSatish Balay 6827bd09bSSatish Balay e-mail: hmt@cs.brown.edu 7827bd09bSSatish Balay 8827bd09bSSatish Balay snail-mail: 9827bd09bSSatish Balay Division of Applied Mathematics 10827bd09bSSatish Balay Brown University 11827bd09bSSatish Balay Providence, RI 02912 12827bd09bSSatish Balay 13827bd09bSSatish Balay Last Modification: 14827bd09bSSatish Balay 6.21.97 15827bd09bSSatish Balay ***********************************ivec.c*************************************/ 16827bd09bSSatish Balay 17c6db04a5SJed Brown #include <../src/ksp/pc/impls/tfs/tfs.h> 18827bd09bSSatish Balay 19827bd09bSSatish Balay /* sorting args ivec.c ivec.c ... */ 20827bd09bSSatish Balay #define SORT_OPT 6 21827bd09bSSatish Balay #define SORT_STACK 50000 22827bd09bSSatish Balay 23827bd09bSSatish Balay 24827bd09bSSatish Balay /* allocate an address and size stack for sorter(s) */ 25827bd09bSSatish Balay static void *offset_stack[2*SORT_STACK]; 2652f87cdaSBarry Smith static PetscInt size_stack[SORT_STACK]; 27827bd09bSSatish Balay 287b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 29ca8e9878SJed Brown PetscInt *PCTFS_ivec_copy(PetscInt *arg1, PetscInt *arg2, PetscInt n) 30827bd09bSSatish Balay { 31827bd09bSSatish Balay while (n--) { *arg1++ = *arg2++; } 32827bd09bSSatish Balay return(arg1); 33827bd09bSSatish Balay } 34827bd09bSSatish Balay 357b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 36ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_zero(PetscInt *arg1, PetscInt n) 37827bd09bSSatish Balay { 383fdc5746SBarry Smith PetscFunctionBegin; 39827bd09bSSatish Balay while (n--) { *arg1++ = 0; } 403fdc5746SBarry Smith PetscFunctionReturn(0); 41827bd09bSSatish Balay } 42827bd09bSSatish Balay 437b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 44ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_set(PetscInt *arg1, PetscInt arg2, PetscInt n) 45827bd09bSSatish Balay { 463fdc5746SBarry Smith PetscFunctionBegin; 47827bd09bSSatish Balay while (n--) { *arg1++ = arg2; } 483fdc5746SBarry Smith PetscFunctionReturn(0); 49827bd09bSSatish Balay } 50827bd09bSSatish Balay 517b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 52ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_max(PetscInt *arg1, PetscInt *arg2, PetscInt n) 53827bd09bSSatish Balay { 543fdc5746SBarry Smith PetscFunctionBegin; 5539945688SSatish Balay while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; } 563fdc5746SBarry Smith PetscFunctionReturn(0); 57827bd09bSSatish Balay } 58827bd09bSSatish Balay 597b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 60ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_min( PetscInt *arg1, PetscInt *arg2, PetscInt n) 61827bd09bSSatish Balay { 623fdc5746SBarry Smith PetscFunctionBegin; 6339945688SSatish Balay while (n--) { *(arg1) = PetscMin(*arg1,*arg2); arg1++; arg2++; } 643fdc5746SBarry Smith PetscFunctionReturn(0); 65827bd09bSSatish Balay } 66827bd09bSSatish Balay 677b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 68ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_mult( PetscInt *arg1, PetscInt *arg2, PetscInt n) 69827bd09bSSatish Balay { 703fdc5746SBarry Smith PetscFunctionBegin; 71827bd09bSSatish Balay while (n--) { *arg1++ *= *arg2++; } 723fdc5746SBarry Smith PetscFunctionReturn(0); 73827bd09bSSatish Balay } 74827bd09bSSatish Balay 757b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 76ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_add( PetscInt *arg1, PetscInt *arg2, PetscInt n) 77827bd09bSSatish Balay { 783fdc5746SBarry Smith PetscFunctionBegin; 79827bd09bSSatish Balay while (n--) { *arg1++ += *arg2++; } 803fdc5746SBarry Smith PetscFunctionReturn(0); 81827bd09bSSatish Balay } 82827bd09bSSatish Balay 837b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 84ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_lxor( PetscInt *arg1, PetscInt *arg2, PetscInt n) 85827bd09bSSatish Balay { 863fdc5746SBarry Smith PetscFunctionBegin; 87827bd09bSSatish Balay while (n--) { *arg1=((*arg1 || *arg2) && !(*arg1 && *arg2)); arg1++; arg2++; } 883fdc5746SBarry Smith PetscFunctionReturn(0); 89827bd09bSSatish Balay } 90827bd09bSSatish Balay 917b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 92ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_xor( PetscInt *arg1, PetscInt *arg2, PetscInt n) 93827bd09bSSatish Balay { 943fdc5746SBarry Smith PetscFunctionBegin; 95827bd09bSSatish Balay while (n--) { *arg1++ ^= *arg2++; } 963fdc5746SBarry Smith PetscFunctionReturn(0); 97827bd09bSSatish Balay } 98827bd09bSSatish Balay 997b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 100ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_or( PetscInt *arg1, PetscInt *arg2, PetscInt n) 101827bd09bSSatish Balay { 1023fdc5746SBarry Smith PetscFunctionBegin; 103827bd09bSSatish Balay while (n--) { *arg1++ |= *arg2++; } 1043fdc5746SBarry Smith PetscFunctionReturn(0); 105827bd09bSSatish Balay } 106827bd09bSSatish Balay 1077b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 108ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_lor( PetscInt *arg1, PetscInt *arg2, PetscInt n) 109827bd09bSSatish Balay { 1103fdc5746SBarry Smith PetscFunctionBegin; 111827bd09bSSatish Balay while (n--) { *arg1 = (*arg1 || *arg2); arg1++; arg2++; } 1123fdc5746SBarry Smith PetscFunctionReturn(0); 113827bd09bSSatish Balay } 114827bd09bSSatish Balay 1157b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 116ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_and( PetscInt *arg1, PetscInt *arg2, PetscInt n) 117827bd09bSSatish Balay { 1183fdc5746SBarry Smith PetscFunctionBegin; 119827bd09bSSatish Balay while (n--) { *arg1++ &= *arg2++; } 1203fdc5746SBarry Smith PetscFunctionReturn(0); 121827bd09bSSatish Balay } 122827bd09bSSatish Balay 1237b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 124ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_land( PetscInt *arg1, PetscInt *arg2, PetscInt n) 125827bd09bSSatish Balay { 1263fdc5746SBarry Smith PetscFunctionBegin; 127827bd09bSSatish Balay while (n--) { *arg1 = (*arg1 && *arg2); arg1++; arg2++; } 1283fdc5746SBarry Smith PetscFunctionReturn(0); 129827bd09bSSatish Balay } 130827bd09bSSatish Balay 1317b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 132ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_and3( PetscInt *arg1, PetscInt *arg2, PetscInt *arg3, PetscInt n) 133827bd09bSSatish Balay { 1343fdc5746SBarry Smith PetscFunctionBegin; 135827bd09bSSatish Balay while (n--) { *arg1++ = (*arg2++ & *arg3++); } 1363fdc5746SBarry Smith PetscFunctionReturn(0); 137827bd09bSSatish Balay } 138827bd09bSSatish Balay 1397b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 140ca8e9878SJed Brown PetscInt PCTFS_ivec_sum( PetscInt *arg1, PetscInt n) 141827bd09bSSatish Balay { 14252f87cdaSBarry Smith PetscInt tmp = 0; 143827bd09bSSatish Balay 144827bd09bSSatish Balay while (n--) {tmp += *arg1++;} 145827bd09bSSatish Balay return(tmp); 146827bd09bSSatish Balay } 147827bd09bSSatish Balay 1487b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 149ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2, PetscInt n, PetscInt *arg3) 150827bd09bSSatish Balay { 15152f87cdaSBarry Smith PetscInt i, j, type; 152827bd09bSSatish Balay 153827bd09bSSatish Balay 1543fdc5746SBarry Smith PetscFunctionBegin; 155827bd09bSSatish Balay /* LATER: if we're really motivated we can sort and then unsort */ 156*db4deed7SKarl Rupp for (i=0;i<n;) { 157827bd09bSSatish Balay /* clump 'em for now */ 158827bd09bSSatish Balay j=i+1; 159827bd09bSSatish Balay type = arg3[i]; 160*db4deed7SKarl Rupp while ((j<n)&&(arg3[j]==type)) { j++; } 161827bd09bSSatish Balay 162827bd09bSSatish Balay /* how many together */ 163827bd09bSSatish Balay j -= i; 164827bd09bSSatish Balay 165827bd09bSSatish Balay /* call appropriate ivec function */ 166*db4deed7SKarl Rupp if (type == GL_MAX) { PCTFS_ivec_max(arg1,arg2,j); } 167*db4deed7SKarl Rupp else if (type == GL_MIN) { PCTFS_ivec_min(arg1,arg2,j); } 168*db4deed7SKarl Rupp else if (type == GL_MULT) { PCTFS_ivec_mult(arg1,arg2,j); } 169*db4deed7SKarl Rupp else if (type == GL_ADD) { PCTFS_ivec_add(arg1,arg2,j); } 170*db4deed7SKarl Rupp else if (type == GL_B_XOR) { PCTFS_ivec_xor(arg1,arg2,j); } 171*db4deed7SKarl Rupp else if (type == GL_B_OR) { PCTFS_ivec_or(arg1,arg2,j); } 172*db4deed7SKarl Rupp else if (type == GL_B_AND) { PCTFS_ivec_and(arg1,arg2,j); } 173*db4deed7SKarl Rupp else if (type == GL_L_XOR) { PCTFS_ivec_lxor(arg1,arg2,j); } 174*db4deed7SKarl Rupp else if (type == GL_L_OR) { PCTFS_ivec_lor(arg1,arg2,j); } 175*db4deed7SKarl Rupp else if (type == GL_L_AND) { PCTFS_ivec_land(arg1,arg2,j); } 176*db4deed7SKarl Rupp else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_ivec_non_uniform()!"); 177827bd09bSSatish Balay 178827bd09bSSatish Balay arg1+=j; arg2+=j; i+=j; 179827bd09bSSatish Balay } 1803fdc5746SBarry Smith PetscFunctionReturn(0); 181827bd09bSSatish Balay } 182827bd09bSSatish Balay 1837b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 184ca8e9878SJed Brown vfp PCTFS_ivec_fct_addr( PetscInt type) 185827bd09bSSatish Balay { 1863fdc5746SBarry Smith PetscFunctionBegin; 187*db4deed7SKarl Rupp if (type == NON_UNIFORM) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_non_uniform); } 188*db4deed7SKarl Rupp else if (type == GL_MAX) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_max); } 189*db4deed7SKarl Rupp else if (type == GL_MIN) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_min); } 190*db4deed7SKarl Rupp else if (type == GL_MULT) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_mult); } 191*db4deed7SKarl Rupp else if (type == GL_ADD) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_add); } 192*db4deed7SKarl Rupp else if (type == GL_B_XOR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_xor); } 193*db4deed7SKarl Rupp else if (type == GL_B_OR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_or); } 194*db4deed7SKarl Rupp else if (type == GL_B_AND) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_and); } 195*db4deed7SKarl Rupp else if (type == GL_L_XOR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lxor); } 196*db4deed7SKarl Rupp else if (type == GL_L_OR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lor); } 197*db4deed7SKarl Rupp else if (type == GL_L_AND) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_land); } 198827bd09bSSatish Balay 199827bd09bSSatish Balay /* catch all ... not good if we get here */ 200827bd09bSSatish Balay return(NULL); 201827bd09bSSatish Balay } 202827bd09bSSatish Balay 2037b1ae94cSBarry Smith /******************************************************************************/ 204ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort( PetscInt *ar, PetscInt size) 205827bd09bSSatish Balay { 20652f87cdaSBarry Smith PetscInt *pi, *pj, temp; 20752f87cdaSBarry Smith PetscInt **top_a = (PetscInt **) offset_stack; 20852f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 209827bd09bSSatish Balay 210827bd09bSSatish Balay 211827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 212827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 213827bd09bSSatish Balay size--; 214827bd09bSSatish Balay 215827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 216*db4deed7SKarl Rupp for (;;) { 217827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 218*db4deed7SKarl Rupp if (size > SORT_OPT) { 219827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 220827bd09bSSatish Balay pi = ar+1; 221827bd09bSSatish Balay pj = ar+size; 222827bd09bSSatish Balay 223827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 224827bd09bSSatish Balay SWAP(*(ar+(size>>1)),*pi) 225827bd09bSSatish Balay 226827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 227827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 228*db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) } 229*db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) } 230*db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) } 231827bd09bSSatish Balay 232827bd09bSSatish Balay /* partition about pivot_value ... */ 233827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 234*db4deed7SKarl Rupp for (;;) { 235827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 236827bd09bSSatish Balay do pi++; while (*pi<*ar); 237827bd09bSSatish Balay do pj--; while (*pj>*ar); 238827bd09bSSatish Balay 239827bd09bSSatish Balay /* if we've crossed we're done */ 240827bd09bSSatish Balay if (pj<pi) break; 241827bd09bSSatish Balay 242827bd09bSSatish Balay /* else swap */ 243827bd09bSSatish Balay SWAP(*pi,*pj) 244827bd09bSSatish Balay } 245827bd09bSSatish Balay 246827bd09bSSatish Balay /* place pivot_value in it's correct location */ 247827bd09bSSatish Balay SWAP(*ar,*pj) 248827bd09bSSatish Balay 249827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 250ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort() :: STACK EXHAUSTED!!!"); 251827bd09bSSatish Balay 252827bd09bSSatish Balay /* push right hand child iff length > 1 */ 253*db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 254827bd09bSSatish Balay *(top_a++) = pi; 255827bd09bSSatish Balay size -= *top_s+2; 256827bd09bSSatish Balay top_s++; 257*db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 258*db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 259827bd09bSSatish Balay ar = *(--top_a); 260827bd09bSSatish Balay size = *(--top_s); 261827bd09bSSatish Balay } 262*db4deed7SKarl Rupp } else { /* else sort small list directly then pop another off stack */ 263827bd09bSSatish Balay 264827bd09bSSatish Balay /* insertion sort for bottom */ 265*db4deed7SKarl Rupp for (pj=ar+1;pj<=ar+size;pj++) { 266827bd09bSSatish Balay temp = *pj; 267*db4deed7SKarl Rupp for (pi=pj-1;pi>=ar;pi--) { 268827bd09bSSatish Balay if (*pi <= temp) break; 269827bd09bSSatish Balay *(pi+1)=*pi; 270827bd09bSSatish Balay } 271827bd09bSSatish Balay *(pi+1)=temp; 272827bd09bSSatish Balay } 273827bd09bSSatish Balay 274827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 2753fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 276827bd09bSSatish Balay 277827bd09bSSatish Balay /* else pop another list from the stack */ 278827bd09bSSatish Balay ar = *(--top_a); 279827bd09bSSatish Balay size = *(--top_s); 280827bd09bSSatish Balay } 281827bd09bSSatish Balay } 2823fdc5746SBarry Smith PetscFunctionReturn(0); 283827bd09bSSatish Balay } 284827bd09bSSatish Balay 2857b1ae94cSBarry Smith /******************************************************************************/ 286ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort_companion( PetscInt *ar, PetscInt *ar2, PetscInt size) 287827bd09bSSatish Balay { 28852f87cdaSBarry Smith PetscInt *pi, *pj, temp, temp2; 28952f87cdaSBarry Smith PetscInt **top_a = (PetscInt **)offset_stack; 29052f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 29152f87cdaSBarry Smith PetscInt *pi2, *pj2; 29252f87cdaSBarry Smith PetscInt mid; 293827bd09bSSatish Balay 2943fdc5746SBarry Smith PetscFunctionBegin; 295827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 296827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 297827bd09bSSatish Balay size--; 298827bd09bSSatish Balay 299827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 300*db4deed7SKarl Rupp for (;;) { 301*db4deed7SKarl Rupp 302827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 303*db4deed7SKarl Rupp if (size > SORT_OPT) { 304*db4deed7SKarl Rupp 305827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 306827bd09bSSatish Balay mid = size>>1; 307827bd09bSSatish Balay pi = ar+1; 308827bd09bSSatish Balay pj = ar+mid; 309827bd09bSSatish Balay pi2 = ar2+1; 310827bd09bSSatish Balay pj2 = ar2+mid; 311827bd09bSSatish Balay 312827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 313827bd09bSSatish Balay SWAP(*pi,*pj) 314827bd09bSSatish Balay SWAP(*pi2,*pj2) 315827bd09bSSatish Balay 316827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 317827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 318827bd09bSSatish Balay pj = ar+size; 319827bd09bSSatish Balay pj2 = ar2+size; 320*db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) SWAP(*pi2,*pj2) } 321*db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) SWAP(*ar2,*pj2) } 322*db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1)) } 323827bd09bSSatish Balay 324827bd09bSSatish Balay /* partition about pivot_value ... */ 325827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 326*db4deed7SKarl Rupp for (;;) { 327827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 328827bd09bSSatish Balay do { pi++; pi2++; } while (*pi<*ar); 329827bd09bSSatish Balay do { pj--; pj2--; } while (*pj>*ar); 330827bd09bSSatish Balay 331827bd09bSSatish Balay /* if we've crossed we're done */ 332827bd09bSSatish Balay if (pj<pi) break; 333827bd09bSSatish Balay 334827bd09bSSatish Balay /* else swap */ 335827bd09bSSatish Balay SWAP(*pi,*pj) 336827bd09bSSatish Balay SWAP(*pi2,*pj2) 337827bd09bSSatish Balay } 338827bd09bSSatish Balay 339827bd09bSSatish Balay /* place pivot_value in it's correct location */ 340827bd09bSSatish Balay SWAP(*ar,*pj) 341827bd09bSSatish Balay SWAP(*ar2,*pj2) 342827bd09bSSatish Balay 343827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 344ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion() :: STACK EXHAUSTED!!!"); 345827bd09bSSatish Balay 346827bd09bSSatish Balay /* push right hand child iff length > 1 */ 347*db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 348827bd09bSSatish Balay *(top_a++) = pi; 349827bd09bSSatish Balay *(top_a++) = pi2; 350827bd09bSSatish Balay size -= *top_s+2; 351827bd09bSSatish Balay top_s++; 352*db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 353*db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 354827bd09bSSatish Balay ar2 = *(--top_a); 355827bd09bSSatish Balay ar = *(--top_a); 356827bd09bSSatish Balay size = *(--top_s); 357827bd09bSSatish Balay } 358*db4deed7SKarl Rupp } else { /* else sort small list directly then pop another off stack */ 359827bd09bSSatish Balay 360827bd09bSSatish Balay /* insertion sort for bottom */ 361*db4deed7SKarl Rupp for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { 362827bd09bSSatish Balay temp = *pj; 363827bd09bSSatish Balay temp2 = *pj2; 364*db4deed7SKarl Rupp for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { 365827bd09bSSatish Balay if (*pi <= temp) break; 366827bd09bSSatish Balay *(pi+1)=*pi; 367827bd09bSSatish Balay *(pi2+1)=*pi2; 368827bd09bSSatish Balay } 369827bd09bSSatish Balay *(pi+1)=temp; 370827bd09bSSatish Balay *(pi2+1)=temp2; 371827bd09bSSatish Balay } 372827bd09bSSatish Balay 373827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 3743fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 375827bd09bSSatish Balay 376827bd09bSSatish Balay /* else pop another list from the stack */ 377827bd09bSSatish Balay ar2 = *(--top_a); 378827bd09bSSatish Balay ar = *(--top_a); 379827bd09bSSatish Balay size = *(--top_s); 380827bd09bSSatish Balay } 381827bd09bSSatish Balay } 3823fdc5746SBarry Smith PetscFunctionReturn(0); 383827bd09bSSatish Balay } 384827bd09bSSatish Balay 3857b1ae94cSBarry Smith /******************************************************************************/ 386ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort_companion_hack( PetscInt *ar, PetscInt **ar2, PetscInt size) 387827bd09bSSatish Balay { 38852f87cdaSBarry Smith PetscInt *pi, *pj, temp, *ptr; 38952f87cdaSBarry Smith PetscInt **top_a = (PetscInt **)offset_stack; 39052f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 39152f87cdaSBarry Smith PetscInt **pi2, **pj2; 39252f87cdaSBarry Smith PetscInt mid; 393827bd09bSSatish Balay 3943fdc5746SBarry Smith PetscFunctionBegin; 395827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 396827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 397827bd09bSSatish Balay size--; 398827bd09bSSatish Balay 399827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 400*db4deed7SKarl Rupp for (;;) { 401*db4deed7SKarl Rupp 402827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 403*db4deed7SKarl Rupp if (size > SORT_OPT) { 404*db4deed7SKarl Rupp 405827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 406827bd09bSSatish Balay mid = size>>1; 407827bd09bSSatish Balay pi = ar+1; 408827bd09bSSatish Balay pj = ar+mid; 409827bd09bSSatish Balay pi2 = ar2+1; 410827bd09bSSatish Balay pj2 = ar2+mid; 411827bd09bSSatish Balay 412827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 413827bd09bSSatish Balay SWAP(*pi,*pj) 414827bd09bSSatish Balay P_SWAP(*pi2,*pj2) 415827bd09bSSatish Balay 416827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 417827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 418827bd09bSSatish Balay pj = ar+size; 419827bd09bSSatish Balay pj2 = ar2+size; 420*db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) } 421*db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) P_SWAP(*ar2,*pj2) } 422*db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1)) } 423827bd09bSSatish Balay 424827bd09bSSatish Balay /* partition about pivot_value ... */ 425827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 426*db4deed7SKarl Rupp for (;;) { 427*db4deed7SKarl Rupp 428827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 429827bd09bSSatish Balay do {pi++; pi2++;} while (*pi<*ar); 430827bd09bSSatish Balay do {pj--; pj2--;} while (*pj>*ar); 431827bd09bSSatish Balay 432827bd09bSSatish Balay /* if we've crossed we're done */ 433827bd09bSSatish Balay if (pj<pi) break; 434827bd09bSSatish Balay 435827bd09bSSatish Balay /* else swap */ 436827bd09bSSatish Balay SWAP(*pi,*pj) 437827bd09bSSatish Balay P_SWAP(*pi2,*pj2) 438827bd09bSSatish Balay } 439827bd09bSSatish Balay 440827bd09bSSatish Balay /* place pivot_value in it's correct location */ 441827bd09bSSatish Balay SWAP(*ar,*pj) 442827bd09bSSatish Balay P_SWAP(*ar2,*pj2) 443827bd09bSSatish Balay 444827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 445ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion_hack() :: STACK EXHAUSTED!!!"); 446827bd09bSSatish Balay 447827bd09bSSatish Balay /* push right hand child iff length > 1 */ 448*db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 449827bd09bSSatish Balay *(top_a++) = pi; 45052f87cdaSBarry Smith *(top_a++) = (PetscInt*) pi2; 451827bd09bSSatish Balay size -= *top_s+2; 452827bd09bSSatish Balay top_s++; 453*db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 454*db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 45552f87cdaSBarry Smith ar2 = (PetscInt **) *(--top_a); 456827bd09bSSatish Balay ar = *(--top_a); 457827bd09bSSatish Balay size = *(--top_s); 458827bd09bSSatish Balay } 459827bd09bSSatish Balay } 460827bd09bSSatish Balay 461*db4deed7SKarl Rupp 462*db4deed7SKarl Rupp else { /* else sort small list directly then pop another off stack */ 463827bd09bSSatish Balay /* insertion sort for bottom */ 464*db4deed7SKarl Rupp for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { 465827bd09bSSatish Balay temp = *pj; 466827bd09bSSatish Balay ptr = *pj2; 467*db4deed7SKarl Rupp for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { 468827bd09bSSatish Balay if (*pi <= temp) break; 469827bd09bSSatish Balay *(pi+1)=*pi; 470827bd09bSSatish Balay *(pi2+1)=*pi2; 471827bd09bSSatish Balay } 472827bd09bSSatish Balay *(pi+1)=temp; 473827bd09bSSatish Balay *(pi2+1)=ptr; 474827bd09bSSatish Balay } 475827bd09bSSatish Balay 476827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 4773fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 478827bd09bSSatish Balay 479827bd09bSSatish Balay /* else pop another list from the stack */ 48052f87cdaSBarry Smith ar2 = (PetscInt **)*(--top_a); 481827bd09bSSatish Balay ar = *(--top_a); 482827bd09bSSatish Balay size = *(--top_s); 483827bd09bSSatish Balay } 484827bd09bSSatish Balay } 4853fdc5746SBarry Smith PetscFunctionReturn(0); 486827bd09bSSatish Balay } 487827bd09bSSatish Balay 4887b1ae94cSBarry Smith /******************************************************************************/ 489ca8e9878SJed Brown PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type) 490827bd09bSSatish Balay { 4913fdc5746SBarry Smith PetscFunctionBegin; 492e7e72b3dSBarry Smith if (type == SORT_INTEGER) { 493ca8e9878SJed Brown if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size); 494ca8e9878SJed Brown else PCTFS_ivec_sort((PetscInt*)ar1,size); 495e7e72b3dSBarry Smith } else if (type == SORT_INT_PTR) { 496ca8e9878SJed Brown if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt **)ar2,size); 497ca8e9878SJed Brown else PCTFS_ivec_sort((PetscInt*)ar1,size); 498ca8e9878SJed Brown } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!"); 4993fdc5746SBarry Smith PetscFunctionReturn(0); 500827bd09bSSatish Balay } 501827bd09bSSatish Balay 5027b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 503ca8e9878SJed Brown PetscInt PCTFS_ivec_linear_search( PetscInt item, PetscInt *list, PetscInt n) 504827bd09bSSatish Balay { 50552f87cdaSBarry Smith PetscInt tmp = n-1; 5063fdc5746SBarry Smith PetscFunctionBegin; 507827bd09bSSatish Balay while (n--) { if (*list++ == item) { return(tmp-n); } } 508827bd09bSSatish Balay return(-1); 509827bd09bSSatish Balay } 510827bd09bSSatish Balay 5117b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 512ca8e9878SJed Brown PetscInt PCTFS_ivec_binary_search( PetscInt item, PetscInt *list, PetscInt rh) 513827bd09bSSatish Balay { 51452f87cdaSBarry Smith PetscInt mid, lh=0; 515827bd09bSSatish Balay 516827bd09bSSatish Balay rh--; 517*db4deed7SKarl Rupp while (lh<=rh) { 518827bd09bSSatish Balay mid = (lh+rh)>>1; 519*db4deed7SKarl Rupp if (*(list+mid) == item) { return(mid); } 520*db4deed7SKarl Rupp if (*(list+mid) > item) { rh = mid-1; } 521*db4deed7SKarl Rupp else { lh = mid+1; } 522827bd09bSSatish Balay } 523827bd09bSSatish Balay return(-1); 524827bd09bSSatish Balay } 525827bd09bSSatish Balay 5267b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 527ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_copy( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 528827bd09bSSatish Balay { 5293fdc5746SBarry Smith PetscFunctionBegin; 530827bd09bSSatish Balay while (n--) { *arg1++ = *arg2++; } 5313fdc5746SBarry Smith PetscFunctionReturn(0); 532827bd09bSSatish Balay } 533827bd09bSSatish Balay 5347b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 535ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_zero( PetscScalar *arg1, PetscInt n) 536827bd09bSSatish Balay { 5373fdc5746SBarry Smith PetscFunctionBegin; 538827bd09bSSatish Balay while (n--) { *arg1++ = 0.0; } 5393fdc5746SBarry Smith PetscFunctionReturn(0); 540827bd09bSSatish Balay } 541827bd09bSSatish Balay 5427b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 543ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_one( PetscScalar *arg1, PetscInt n) 544827bd09bSSatish Balay { 5453fdc5746SBarry Smith PetscFunctionBegin; 546827bd09bSSatish Balay while (n--) { *arg1++ = 1.0; } 5473fdc5746SBarry Smith PetscFunctionReturn(0); 548827bd09bSSatish Balay } 549827bd09bSSatish Balay 5507b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 551ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_set( PetscScalar *arg1, PetscScalar arg2, PetscInt n) 552827bd09bSSatish Balay { 5533fdc5746SBarry Smith PetscFunctionBegin; 554827bd09bSSatish Balay while (n--) { *arg1++ = arg2; } 5553fdc5746SBarry Smith PetscFunctionReturn(0); 556827bd09bSSatish Balay } 557827bd09bSSatish Balay 5587b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 559ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_scale( PetscScalar *arg1, PetscScalar arg2, PetscInt n) 560827bd09bSSatish Balay { 5613fdc5746SBarry Smith PetscFunctionBegin; 562827bd09bSSatish Balay while (n--) { *arg1++ *= arg2; } 5633fdc5746SBarry Smith PetscFunctionReturn(0); 564827bd09bSSatish Balay } 565827bd09bSSatish Balay 5667b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 567ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_add( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 568827bd09bSSatish Balay { 5693fdc5746SBarry Smith PetscFunctionBegin; 570827bd09bSSatish Balay while (n--) { *arg1++ += *arg2++; } 5713fdc5746SBarry Smith PetscFunctionReturn(0); 572827bd09bSSatish Balay } 573827bd09bSSatish Balay 5747b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 575ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_mult( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 576827bd09bSSatish Balay { 5773fdc5746SBarry Smith PetscFunctionBegin; 578827bd09bSSatish Balay while (n--) { *arg1++ *= *arg2++; } 5793fdc5746SBarry Smith PetscFunctionReturn(0); 580827bd09bSSatish Balay } 581827bd09bSSatish Balay 5827b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 583ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_max( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 584827bd09bSSatish Balay { 5853fdc5746SBarry Smith PetscFunctionBegin; 58639945688SSatish Balay while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; } 5873fdc5746SBarry Smith PetscFunctionReturn(0); 588827bd09bSSatish Balay } 589827bd09bSSatish Balay 5907b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 591ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_max_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 592827bd09bSSatish Balay { 5933fdc5746SBarry Smith PetscFunctionBegin; 594827bd09bSSatish Balay while (n--) { *arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++; } 5953fdc5746SBarry Smith PetscFunctionReturn(0); 596827bd09bSSatish Balay } 597827bd09bSSatish Balay 5987b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 599ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_min( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 600827bd09bSSatish Balay { 6013fdc5746SBarry Smith PetscFunctionBegin; 60239945688SSatish Balay while (n--) { *arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++; } 6033fdc5746SBarry Smith PetscFunctionReturn(0); 604827bd09bSSatish Balay } 605827bd09bSSatish Balay 6067b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 607ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_min_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 608827bd09bSSatish Balay { 6093fdc5746SBarry Smith PetscFunctionBegin; 610827bd09bSSatish Balay while (n--) { *arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++; } 6113fdc5746SBarry Smith PetscFunctionReturn(0); 612827bd09bSSatish Balay } 613827bd09bSSatish Balay 6147b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 615ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_exists( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 616827bd09bSSatish Balay { 6173fdc5746SBarry Smith PetscFunctionBegin; 618827bd09bSSatish Balay while (n--) { *arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++; } 6193fdc5746SBarry Smith PetscFunctionReturn(0); 620827bd09bSSatish Balay } 621827bd09bSSatish Balay 6227b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 623ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2, PetscInt n, PetscInt *arg3) 624827bd09bSSatish Balay { 62552f87cdaSBarry Smith PetscInt i, j, type; 626827bd09bSSatish Balay 6273fdc5746SBarry Smith PetscFunctionBegin; 628827bd09bSSatish Balay /* LATER: if we're really motivated we can sort and then unsort */ 629*db4deed7SKarl Rupp for (i=0;i<n;) { 630*db4deed7SKarl Rupp 631827bd09bSSatish Balay /* clump 'em for now */ 632827bd09bSSatish Balay j=i+1; 633827bd09bSSatish Balay type = arg3[i]; 634*db4deed7SKarl Rupp while ((j<n)&&(arg3[j]==type)) { j++; } 635827bd09bSSatish Balay 636827bd09bSSatish Balay /* how many together */ 637827bd09bSSatish Balay j -= i; 638827bd09bSSatish Balay 639827bd09bSSatish Balay /* call appropriate ivec function */ 640*db4deed7SKarl Rupp if (type == GL_MAX) { PCTFS_rvec_max(arg1,arg2,j); } 641*db4deed7SKarl Rupp else if (type == GL_MIN) { PCTFS_rvec_min(arg1,arg2,j); } 642*db4deed7SKarl Rupp else if (type == GL_MULT) { PCTFS_rvec_mult(arg1,arg2,j); } 643*db4deed7SKarl Rupp else if (type == GL_ADD) { PCTFS_rvec_add(arg1,arg2,j); } 644*db4deed7SKarl Rupp else if (type == GL_MAX_ABS) { PCTFS_rvec_max_abs(arg1,arg2,j); } 645*db4deed7SKarl Rupp else if (type == GL_MIN_ABS) { PCTFS_rvec_min_abs(arg1,arg2,j); } 646*db4deed7SKarl Rupp else if (type == GL_EXISTS) { PCTFS_rvec_exists(arg1,arg2,j); } 647ca8e9878SJed Brown else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_rvec_non_uniform()!"); 648827bd09bSSatish Balay 649827bd09bSSatish Balay arg1+=j; arg2+=j; i+=j; 650827bd09bSSatish Balay } 6513fdc5746SBarry Smith PetscFunctionReturn(0); 652827bd09bSSatish Balay } 653827bd09bSSatish Balay 6547b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 655ca8e9878SJed Brown vfp PCTFS_rvec_fct_addr( PetscInt type) 656827bd09bSSatish Balay { 657*db4deed7SKarl Rupp if (type == NON_UNIFORM) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_non_uniform); } 658*db4deed7SKarl Rupp else if (type == GL_MAX) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max); } 659*db4deed7SKarl Rupp else if (type == GL_MIN) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min); } 660*db4deed7SKarl Rupp else if (type == GL_MULT) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_mult); } 661*db4deed7SKarl Rupp else if (type == GL_ADD) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_add); } 662*db4deed7SKarl Rupp else if (type == GL_MAX_ABS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max_abs); } 663*db4deed7SKarl Rupp else if (type == GL_MIN_ABS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min_abs); } 664*db4deed7SKarl Rupp else if (type == GL_EXISTS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_exists); } 665827bd09bSSatish Balay 666827bd09bSSatish Balay /* catch all ... not good if we get here */ 667827bd09bSSatish Balay return(NULL); 668827bd09bSSatish Balay } 669827bd09bSSatish Balay 670827bd09bSSatish Balay 671827bd09bSSatish Balay 672827bd09bSSatish Balay 673827bd09bSSatish Balay 674827bd09bSSatish Balay 675