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 1533fdc5746SBarry Smith PetscFunctionBegin; 154827bd09bSSatish Balay /* LATER: if we're really motivated we can sort and then unsort */ 155db4deed7SKarl Rupp for (i=0;i<n;) { 156827bd09bSSatish Balay /* clump 'em for now */ 157827bd09bSSatish Balay j=i+1; 158827bd09bSSatish Balay type = arg3[i]; 159db4deed7SKarl Rupp while ((j<n)&&(arg3[j]==type)) { j++; } 160827bd09bSSatish Balay 161827bd09bSSatish Balay /* how many together */ 162827bd09bSSatish Balay j -= i; 163827bd09bSSatish Balay 164827bd09bSSatish Balay /* call appropriate ivec function */ 165db4deed7SKarl Rupp if (type == GL_MAX) { PCTFS_ivec_max(arg1,arg2,j); } 166db4deed7SKarl Rupp else if (type == GL_MIN) { PCTFS_ivec_min(arg1,arg2,j); } 167db4deed7SKarl Rupp else if (type == GL_MULT) { PCTFS_ivec_mult(arg1,arg2,j); } 168db4deed7SKarl Rupp else if (type == GL_ADD) { PCTFS_ivec_add(arg1,arg2,j); } 169db4deed7SKarl Rupp else if (type == GL_B_XOR) { PCTFS_ivec_xor(arg1,arg2,j); } 170db4deed7SKarl Rupp else if (type == GL_B_OR) { PCTFS_ivec_or(arg1,arg2,j); } 171db4deed7SKarl Rupp else if (type == GL_B_AND) { PCTFS_ivec_and(arg1,arg2,j); } 172db4deed7SKarl Rupp else if (type == GL_L_XOR) { PCTFS_ivec_lxor(arg1,arg2,j); } 173db4deed7SKarl Rupp else if (type == GL_L_OR) { PCTFS_ivec_lor(arg1,arg2,j); } 174db4deed7SKarl Rupp else if (type == GL_L_AND) { PCTFS_ivec_land(arg1,arg2,j); } 175db4deed7SKarl Rupp else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_ivec_non_uniform()!"); 176827bd09bSSatish Balay 177827bd09bSSatish Balay arg1+=j; arg2+=j; i+=j; 178827bd09bSSatish Balay } 1793fdc5746SBarry Smith PetscFunctionReturn(0); 180827bd09bSSatish Balay } 181827bd09bSSatish Balay 1827b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 183ca8e9878SJed Brown vfp PCTFS_ivec_fct_addr(PetscInt type) 184827bd09bSSatish Balay { 185db4deed7SKarl Rupp if (type == NON_UNIFORM) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_non_uniform); } 186db4deed7SKarl Rupp else if (type == GL_MAX) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_max); } 187db4deed7SKarl Rupp else if (type == GL_MIN) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_min); } 188db4deed7SKarl Rupp else if (type == GL_MULT) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_mult); } 189db4deed7SKarl Rupp else if (type == GL_ADD) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_add); } 190db4deed7SKarl Rupp else if (type == GL_B_XOR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_xor); } 191db4deed7SKarl Rupp else if (type == GL_B_OR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_or); } 192db4deed7SKarl Rupp else if (type == GL_B_AND) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_and); } 193db4deed7SKarl Rupp else if (type == GL_L_XOR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lxor); } 194db4deed7SKarl Rupp else if (type == GL_L_OR) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_lor); } 195db4deed7SKarl Rupp else if (type == GL_L_AND) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_ivec_land); } 196827bd09bSSatish Balay 197827bd09bSSatish Balay /* catch all ... not good if we get here */ 198827bd09bSSatish Balay return(NULL); 199827bd09bSSatish Balay } 200827bd09bSSatish Balay 2017b1ae94cSBarry Smith /******************************************************************************/ 202ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort(PetscInt *ar, PetscInt size) 203827bd09bSSatish Balay { 20452f87cdaSBarry Smith PetscInt *pi, *pj, temp; 20552f87cdaSBarry Smith PetscInt **top_a = (PetscInt **) offset_stack; 20652f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 207827bd09bSSatish Balay 208827bd09bSSatish Balay 209827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 210827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 211827bd09bSSatish Balay size--; 212827bd09bSSatish Balay 213827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 214db4deed7SKarl Rupp for (;;) { 215827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 216db4deed7SKarl Rupp if (size > SORT_OPT) { 217827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 218827bd09bSSatish Balay pi = ar+1; 219827bd09bSSatish Balay pj = ar+size; 220827bd09bSSatish Balay 221827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 222827bd09bSSatish Balay SWAP(*(ar+(size>>1)),*pi) 223827bd09bSSatish Balay 224827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 225827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 226db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) } 227db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) } 228db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) } 229827bd09bSSatish Balay 230827bd09bSSatish Balay /* partition about pivot_value ... */ 231827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 232db4deed7SKarl Rupp for (;;) { 233827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 234827bd09bSSatish Balay do pi++; while (*pi<*ar); 235827bd09bSSatish Balay do pj--; while (*pj>*ar); 236827bd09bSSatish Balay 237827bd09bSSatish Balay /* if we've crossed we're done */ 238827bd09bSSatish Balay if (pj<pi) break; 239827bd09bSSatish Balay 240827bd09bSSatish Balay /* else swap */ 241827bd09bSSatish Balay SWAP(*pi,*pj) 242827bd09bSSatish Balay } 243827bd09bSSatish Balay 244827bd09bSSatish Balay /* place pivot_value in it's correct location */ 245827bd09bSSatish Balay SWAP(*ar,*pj) 246827bd09bSSatish Balay 247827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 248ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort() :: STACK EXHAUSTED!!!"); 249827bd09bSSatish Balay 250827bd09bSSatish Balay /* push right hand child iff length > 1 */ 251db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 252827bd09bSSatish Balay *(top_a++) = pi; 253827bd09bSSatish Balay size -= *top_s+2; 254827bd09bSSatish Balay top_s++; 255db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 256db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 257827bd09bSSatish Balay ar = *(--top_a); 258827bd09bSSatish Balay size = *(--top_s); 259827bd09bSSatish Balay } 260db4deed7SKarl Rupp } else { /* else sort small list directly then pop another off stack */ 261827bd09bSSatish Balay 262827bd09bSSatish Balay /* insertion sort for bottom */ 263db4deed7SKarl Rupp for (pj=ar+1;pj<=ar+size;pj++) { 264827bd09bSSatish Balay temp = *pj; 265db4deed7SKarl Rupp for (pi=pj-1;pi>=ar;pi--) { 266827bd09bSSatish Balay if (*pi <= temp) break; 267827bd09bSSatish Balay *(pi+1)=*pi; 268827bd09bSSatish Balay } 269827bd09bSSatish Balay *(pi+1)=temp; 270827bd09bSSatish Balay } 271827bd09bSSatish Balay 272827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 2733fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 274827bd09bSSatish Balay 275827bd09bSSatish Balay /* else pop another list from the stack */ 276827bd09bSSatish Balay ar = *(--top_a); 277827bd09bSSatish Balay size = *(--top_s); 278827bd09bSSatish Balay } 279827bd09bSSatish Balay } 2803fdc5746SBarry Smith PetscFunctionReturn(0); 281827bd09bSSatish Balay } 282827bd09bSSatish Balay 2837b1ae94cSBarry Smith /******************************************************************************/ 284ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort_companion(PetscInt *ar, PetscInt *ar2, PetscInt size) 285827bd09bSSatish Balay { 28652f87cdaSBarry Smith PetscInt *pi, *pj, temp, temp2; 28752f87cdaSBarry Smith PetscInt **top_a = (PetscInt **)offset_stack; 28852f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 28952f87cdaSBarry Smith PetscInt *pi2, *pj2; 29052f87cdaSBarry Smith PetscInt mid; 291827bd09bSSatish Balay 2923fdc5746SBarry Smith PetscFunctionBegin; 293827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 294827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 295827bd09bSSatish Balay size--; 296827bd09bSSatish Balay 297827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 298db4deed7SKarl Rupp for (;;) { 299db4deed7SKarl Rupp 300827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 301db4deed7SKarl Rupp if (size > SORT_OPT) { 302db4deed7SKarl Rupp 303827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 304827bd09bSSatish Balay mid = size>>1; 305827bd09bSSatish Balay pi = ar+1; 306827bd09bSSatish Balay pj = ar+mid; 307827bd09bSSatish Balay pi2 = ar2+1; 308827bd09bSSatish Balay pj2 = ar2+mid; 309827bd09bSSatish Balay 310827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 311827bd09bSSatish Balay SWAP(*pi,*pj) 312827bd09bSSatish Balay SWAP(*pi2,*pj2) 313827bd09bSSatish Balay 314827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 315827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 316827bd09bSSatish Balay pj = ar+size; 317827bd09bSSatish Balay pj2 = ar2+size; 318db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) SWAP(*pi2,*pj2) } 319db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) SWAP(*ar2,*pj2) } 320db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1)) } 321827bd09bSSatish Balay 322827bd09bSSatish Balay /* partition about pivot_value ... */ 323827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 324db4deed7SKarl Rupp for (;;) { 325827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 326827bd09bSSatish Balay do { pi++; pi2++; } while (*pi<*ar); 327827bd09bSSatish Balay do { pj--; pj2--; } while (*pj>*ar); 328827bd09bSSatish Balay 329827bd09bSSatish Balay /* if we've crossed we're done */ 330827bd09bSSatish Balay if (pj<pi) break; 331827bd09bSSatish Balay 332827bd09bSSatish Balay /* else swap */ 333827bd09bSSatish Balay SWAP(*pi,*pj) 334827bd09bSSatish Balay SWAP(*pi2,*pj2) 335827bd09bSSatish Balay } 336827bd09bSSatish Balay 337827bd09bSSatish Balay /* place pivot_value in it's correct location */ 338827bd09bSSatish Balay SWAP(*ar,*pj) 339827bd09bSSatish Balay SWAP(*ar2,*pj2) 340827bd09bSSatish Balay 341827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 342ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion() :: STACK EXHAUSTED!!!"); 343827bd09bSSatish Balay 344827bd09bSSatish Balay /* push right hand child iff length > 1 */ 345db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 346827bd09bSSatish Balay *(top_a++) = pi; 347827bd09bSSatish Balay *(top_a++) = pi2; 348827bd09bSSatish Balay size -= *top_s+2; 349827bd09bSSatish Balay top_s++; 350db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 351db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 352827bd09bSSatish Balay ar2 = *(--top_a); 353827bd09bSSatish Balay ar = *(--top_a); 354827bd09bSSatish Balay size = *(--top_s); 355827bd09bSSatish Balay } 356db4deed7SKarl Rupp } else { /* else sort small list directly then pop another off stack */ 357827bd09bSSatish Balay 358827bd09bSSatish Balay /* insertion sort for bottom */ 359db4deed7SKarl Rupp for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { 360827bd09bSSatish Balay temp = *pj; 361827bd09bSSatish Balay temp2 = *pj2; 362db4deed7SKarl Rupp for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { 363827bd09bSSatish Balay if (*pi <= temp) break; 364827bd09bSSatish Balay *(pi+1)=*pi; 365827bd09bSSatish Balay *(pi2+1)=*pi2; 366827bd09bSSatish Balay } 367827bd09bSSatish Balay *(pi+1)=temp; 368827bd09bSSatish Balay *(pi2+1)=temp2; 369827bd09bSSatish Balay } 370827bd09bSSatish Balay 371827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 3723fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 373827bd09bSSatish Balay 374827bd09bSSatish Balay /* else pop another list from the stack */ 375827bd09bSSatish Balay ar2 = *(--top_a); 376827bd09bSSatish Balay ar = *(--top_a); 377827bd09bSSatish Balay size = *(--top_s); 378827bd09bSSatish Balay } 379827bd09bSSatish Balay } 3803fdc5746SBarry Smith PetscFunctionReturn(0); 381827bd09bSSatish Balay } 382827bd09bSSatish Balay 3837b1ae94cSBarry Smith /******************************************************************************/ 384ca8e9878SJed Brown PetscErrorCode PCTFS_ivec_sort_companion_hack(PetscInt *ar, PetscInt **ar2, PetscInt size) 385827bd09bSSatish Balay { 38652f87cdaSBarry Smith PetscInt *pi, *pj, temp, *ptr; 38752f87cdaSBarry Smith PetscInt **top_a = (PetscInt **)offset_stack; 38852f87cdaSBarry Smith PetscInt *top_s = size_stack, *bottom_s = size_stack; 38952f87cdaSBarry Smith PetscInt **pi2, **pj2; 39052f87cdaSBarry Smith PetscInt mid; 391827bd09bSSatish Balay 3923fdc5746SBarry Smith PetscFunctionBegin; 393827bd09bSSatish Balay /* we're really interested in the offset of the last element */ 394827bd09bSSatish Balay /* ==> length of the list is now size + 1 */ 395827bd09bSSatish Balay size--; 396827bd09bSSatish Balay 397827bd09bSSatish Balay /* do until we're done ... return when stack is exhausted */ 398db4deed7SKarl Rupp for (;;) { 399db4deed7SKarl Rupp 400827bd09bSSatish Balay /* if list is large enough use quicksort partition exchange code */ 401db4deed7SKarl Rupp if (size > SORT_OPT) { 402db4deed7SKarl Rupp 403827bd09bSSatish Balay /* start up pointer at element 1 and down at size */ 404827bd09bSSatish Balay mid = size>>1; 405827bd09bSSatish Balay pi = ar+1; 406827bd09bSSatish Balay pj = ar+mid; 407827bd09bSSatish Balay pi2 = ar2+1; 408827bd09bSSatish Balay pj2 = ar2+mid; 409827bd09bSSatish Balay 410827bd09bSSatish Balay /* find middle element in list and swap w/ element 1 */ 411827bd09bSSatish Balay SWAP(*pi,*pj) 412827bd09bSSatish Balay P_SWAP(*pi2,*pj2) 413827bd09bSSatish Balay 414827bd09bSSatish Balay /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ 415827bd09bSSatish Balay /* note ==> pivot_value in index 0 */ 416827bd09bSSatish Balay pj = ar+size; 417827bd09bSSatish Balay pj2 = ar2+size; 418db4deed7SKarl Rupp if (*pi > *pj) { SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) } 419db4deed7SKarl Rupp if (*ar > *pj) { SWAP(*ar,*pj) P_SWAP(*ar2,*pj2) } 420db4deed7SKarl Rupp else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1)) } 421827bd09bSSatish Balay 422827bd09bSSatish Balay /* partition about pivot_value ... */ 423827bd09bSSatish Balay /* note lists of length 2 are not guaranteed to be sorted */ 424db4deed7SKarl Rupp for (;;) { 425db4deed7SKarl Rupp 426827bd09bSSatish Balay /* walk up ... and down ... swap if equal to pivot! */ 427827bd09bSSatish Balay do {pi++; pi2++;} while (*pi<*ar); 428827bd09bSSatish Balay do {pj--; pj2--;} while (*pj>*ar); 429827bd09bSSatish Balay 430827bd09bSSatish Balay /* if we've crossed we're done */ 431827bd09bSSatish Balay if (pj<pi) break; 432827bd09bSSatish Balay 433827bd09bSSatish Balay /* else swap */ 434827bd09bSSatish Balay SWAP(*pi,*pj) 435827bd09bSSatish Balay P_SWAP(*pi2,*pj2) 436827bd09bSSatish Balay } 437827bd09bSSatish Balay 438827bd09bSSatish Balay /* place pivot_value in it's correct location */ 439827bd09bSSatish Balay SWAP(*ar,*pj) 440827bd09bSSatish Balay P_SWAP(*ar2,*pj2) 441827bd09bSSatish Balay 442827bd09bSSatish Balay /* test stack_size to see if we've exhausted our stack */ 443ca8e9878SJed Brown if (top_s-bottom_s >= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion_hack() :: STACK EXHAUSTED!!!"); 444827bd09bSSatish Balay 445827bd09bSSatish Balay /* push right hand child iff length > 1 */ 446db4deed7SKarl Rupp if ((*top_s = size-((PetscInt) (pi-ar)))) { 447827bd09bSSatish Balay *(top_a++) = pi; 44852f87cdaSBarry Smith *(top_a++) = (PetscInt*) pi2; 449827bd09bSSatish Balay size -= *top_s+2; 450827bd09bSSatish Balay top_s++; 451db4deed7SKarl Rupp } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ 452db4deed7SKarl Rupp else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ 45352f87cdaSBarry Smith ar2 = (PetscInt **) *(--top_a); 454827bd09bSSatish Balay ar = *(--top_a); 455827bd09bSSatish Balay size = *(--top_s); 456827bd09bSSatish Balay } 457827bd09bSSatish Balay } 458827bd09bSSatish Balay 459db4deed7SKarl Rupp 460db4deed7SKarl Rupp else { /* else sort small list directly then pop another off stack */ 461827bd09bSSatish Balay /* insertion sort for bottom */ 462db4deed7SKarl Rupp for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { 463827bd09bSSatish Balay temp = *pj; 464827bd09bSSatish Balay ptr = *pj2; 465db4deed7SKarl Rupp for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { 466827bd09bSSatish Balay if (*pi <= temp) break; 467827bd09bSSatish Balay *(pi+1)=*pi; 468827bd09bSSatish Balay *(pi2+1)=*pi2; 469827bd09bSSatish Balay } 470827bd09bSSatish Balay *(pi+1)=temp; 471827bd09bSSatish Balay *(pi2+1)=ptr; 472827bd09bSSatish Balay } 473827bd09bSSatish Balay 474827bd09bSSatish Balay /* check to see if stack is exhausted ==> DONE */ 4753fdc5746SBarry Smith if (top_s==bottom_s) PetscFunctionReturn(0); 476827bd09bSSatish Balay 477827bd09bSSatish Balay /* else pop another list from the stack */ 47852f87cdaSBarry Smith ar2 = (PetscInt **)*(--top_a); 479827bd09bSSatish Balay ar = *(--top_a); 480827bd09bSSatish Balay size = *(--top_s); 481827bd09bSSatish Balay } 482827bd09bSSatish Balay } 4833fdc5746SBarry Smith PetscFunctionReturn(0); 484827bd09bSSatish Balay } 485827bd09bSSatish Balay 4867b1ae94cSBarry Smith /******************************************************************************/ 487ca8e9878SJed Brown PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type) 488827bd09bSSatish Balay { 4893fdc5746SBarry Smith PetscFunctionBegin; 490e7e72b3dSBarry Smith if (type == SORT_INTEGER) { 491ca8e9878SJed Brown if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size); 492ca8e9878SJed Brown else PCTFS_ivec_sort((PetscInt*)ar1,size); 493e7e72b3dSBarry Smith } else if (type == SORT_INT_PTR) { 494ca8e9878SJed Brown if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt **)ar2,size); 495ca8e9878SJed Brown else PCTFS_ivec_sort((PetscInt*)ar1,size); 496ca8e9878SJed Brown } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!"); 4973fdc5746SBarry Smith PetscFunctionReturn(0); 498827bd09bSSatish Balay } 499827bd09bSSatish Balay 5007b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 501ca8e9878SJed Brown PetscInt PCTFS_ivec_linear_search(PetscInt item, PetscInt *list, PetscInt n) 502827bd09bSSatish Balay { 50352f87cdaSBarry Smith PetscInt tmp = n-1; 504*5fd66863SKarl Rupp 5053fdc5746SBarry Smith PetscFunctionBegin; 506827bd09bSSatish Balay while (n--) { if (*list++ == item) { return(tmp-n); } } 507827bd09bSSatish Balay return(-1); 508827bd09bSSatish Balay } 509827bd09bSSatish Balay 5107b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 511ca8e9878SJed Brown PetscInt PCTFS_ivec_binary_search(PetscInt item, PetscInt *list, PetscInt rh) 512827bd09bSSatish Balay { 51352f87cdaSBarry Smith PetscInt mid, lh=0; 514827bd09bSSatish Balay 515827bd09bSSatish Balay rh--; 516db4deed7SKarl Rupp while (lh<=rh) { 517827bd09bSSatish Balay mid = (lh+rh)>>1; 518db4deed7SKarl Rupp if (*(list+mid) == item) { return(mid); } 519db4deed7SKarl Rupp if (*(list+mid) > item) { rh = mid-1; } 520db4deed7SKarl Rupp else { lh = mid+1; } 521827bd09bSSatish Balay } 522827bd09bSSatish Balay return(-1); 523827bd09bSSatish Balay } 524827bd09bSSatish Balay 5257b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 526ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_copy(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 527827bd09bSSatish Balay { 5283fdc5746SBarry Smith PetscFunctionBegin; 529827bd09bSSatish Balay while (n--) { *arg1++ = *arg2++; } 5303fdc5746SBarry Smith PetscFunctionReturn(0); 531827bd09bSSatish Balay } 532827bd09bSSatish Balay 5337b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 534ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_zero(PetscScalar *arg1, PetscInt n) 535827bd09bSSatish Balay { 5363fdc5746SBarry Smith PetscFunctionBegin; 537827bd09bSSatish Balay while (n--) { *arg1++ = 0.0; } 5383fdc5746SBarry Smith PetscFunctionReturn(0); 539827bd09bSSatish Balay } 540827bd09bSSatish Balay 5417b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 542ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_one(PetscScalar *arg1, PetscInt n) 543827bd09bSSatish Balay { 5443fdc5746SBarry Smith PetscFunctionBegin; 545827bd09bSSatish Balay while (n--) { *arg1++ = 1.0; } 5463fdc5746SBarry Smith PetscFunctionReturn(0); 547827bd09bSSatish Balay } 548827bd09bSSatish Balay 5497b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 550ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_set(PetscScalar *arg1, PetscScalar arg2, PetscInt n) 551827bd09bSSatish Balay { 5523fdc5746SBarry Smith PetscFunctionBegin; 553827bd09bSSatish Balay while (n--) { *arg1++ = arg2; } 5543fdc5746SBarry Smith PetscFunctionReturn(0); 555827bd09bSSatish Balay } 556827bd09bSSatish Balay 5577b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 558ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_scale(PetscScalar *arg1, PetscScalar arg2, PetscInt n) 559827bd09bSSatish Balay { 5603fdc5746SBarry Smith PetscFunctionBegin; 561827bd09bSSatish Balay while (n--) { *arg1++ *= arg2; } 5623fdc5746SBarry Smith PetscFunctionReturn(0); 563827bd09bSSatish Balay } 564827bd09bSSatish Balay 5657b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 566ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_add(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 567827bd09bSSatish Balay { 5683fdc5746SBarry Smith PetscFunctionBegin; 569827bd09bSSatish Balay while (n--) { *arg1++ += *arg2++; } 5703fdc5746SBarry Smith PetscFunctionReturn(0); 571827bd09bSSatish Balay } 572827bd09bSSatish Balay 5737b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 574ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_mult(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 575827bd09bSSatish Balay { 5763fdc5746SBarry Smith PetscFunctionBegin; 577827bd09bSSatish Balay while (n--) { *arg1++ *= *arg2++; } 5783fdc5746SBarry Smith PetscFunctionReturn(0); 579827bd09bSSatish Balay } 580827bd09bSSatish Balay 5817b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 582ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_max(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 583827bd09bSSatish Balay { 5843fdc5746SBarry Smith PetscFunctionBegin; 58539945688SSatish Balay while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; } 5863fdc5746SBarry Smith PetscFunctionReturn(0); 587827bd09bSSatish Balay } 588827bd09bSSatish Balay 5897b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 590ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_max_abs(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 591827bd09bSSatish Balay { 5923fdc5746SBarry Smith PetscFunctionBegin; 593827bd09bSSatish Balay while (n--) { *arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++; } 5943fdc5746SBarry Smith PetscFunctionReturn(0); 595827bd09bSSatish Balay } 596827bd09bSSatish Balay 5977b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 598ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_min(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 599827bd09bSSatish Balay { 6003fdc5746SBarry Smith PetscFunctionBegin; 60139945688SSatish Balay while (n--) { *arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++; } 6023fdc5746SBarry Smith PetscFunctionReturn(0); 603827bd09bSSatish Balay } 604827bd09bSSatish Balay 6057b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 606ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_min_abs(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 607827bd09bSSatish Balay { 6083fdc5746SBarry Smith PetscFunctionBegin; 609827bd09bSSatish Balay while (n--) { *arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++; } 6103fdc5746SBarry Smith PetscFunctionReturn(0); 611827bd09bSSatish Balay } 612827bd09bSSatish Balay 6137b1ae94cSBarry Smith /*********************************ivec.c*************************************/ 614ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_exists(PetscScalar *arg1, PetscScalar *arg2, PetscInt n) 615827bd09bSSatish Balay { 6163fdc5746SBarry Smith PetscFunctionBegin; 617827bd09bSSatish Balay while (n--) { *arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++; } 6183fdc5746SBarry Smith PetscFunctionReturn(0); 619827bd09bSSatish Balay } 620827bd09bSSatish Balay 6217b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 622ca8e9878SJed Brown PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2, PetscInt n, PetscInt *arg3) 623827bd09bSSatish Balay { 62452f87cdaSBarry Smith PetscInt i, j, type; 625827bd09bSSatish Balay 6263fdc5746SBarry Smith PetscFunctionBegin; 627827bd09bSSatish Balay /* LATER: if we're really motivated we can sort and then unsort */ 628db4deed7SKarl Rupp for (i=0;i<n;) { 629db4deed7SKarl Rupp 630827bd09bSSatish Balay /* clump 'em for now */ 631827bd09bSSatish Balay j=i+1; 632827bd09bSSatish Balay type = arg3[i]; 633db4deed7SKarl Rupp while ((j<n)&&(arg3[j]==type)) { j++; } 634827bd09bSSatish Balay 635827bd09bSSatish Balay /* how many together */ 636827bd09bSSatish Balay j -= i; 637827bd09bSSatish Balay 638827bd09bSSatish Balay /* call appropriate ivec function */ 639db4deed7SKarl Rupp if (type == GL_MAX) { PCTFS_rvec_max(arg1,arg2,j); } 640db4deed7SKarl Rupp else if (type == GL_MIN) { PCTFS_rvec_min(arg1,arg2,j); } 641db4deed7SKarl Rupp else if (type == GL_MULT) { PCTFS_rvec_mult(arg1,arg2,j); } 642db4deed7SKarl Rupp else if (type == GL_ADD) { PCTFS_rvec_add(arg1,arg2,j); } 643db4deed7SKarl Rupp else if (type == GL_MAX_ABS) { PCTFS_rvec_max_abs(arg1,arg2,j); } 644db4deed7SKarl Rupp else if (type == GL_MIN_ABS) { PCTFS_rvec_min_abs(arg1,arg2,j); } 645db4deed7SKarl Rupp else if (type == GL_EXISTS) { PCTFS_rvec_exists(arg1,arg2,j); } 646ca8e9878SJed Brown else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"unrecognized type passed to PCTFS_rvec_non_uniform()!"); 647827bd09bSSatish Balay 648827bd09bSSatish Balay arg1+=j; arg2+=j; i+=j; 649827bd09bSSatish Balay } 6503fdc5746SBarry Smith PetscFunctionReturn(0); 651827bd09bSSatish Balay } 652827bd09bSSatish Balay 6537b1ae94cSBarry Smith /***********************************ivec.c*************************************/ 654ca8e9878SJed Brown vfp PCTFS_rvec_fct_addr(PetscInt type) 655827bd09bSSatish Balay { 656db4deed7SKarl Rupp if (type == NON_UNIFORM) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_non_uniform); } 657db4deed7SKarl Rupp else if (type == GL_MAX) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max); } 658db4deed7SKarl Rupp else if (type == GL_MIN) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min); } 659db4deed7SKarl Rupp else if (type == GL_MULT) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_mult); } 660db4deed7SKarl Rupp else if (type == GL_ADD) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_add); } 661db4deed7SKarl Rupp else if (type == GL_MAX_ABS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_max_abs); } 662db4deed7SKarl Rupp else if (type == GL_MIN_ABS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_min_abs); } 663db4deed7SKarl Rupp else if (type == GL_EXISTS) { return((PetscErrorCode (*)(void*, void *, PetscInt, ...))&PCTFS_rvec_exists); } 664827bd09bSSatish Balay 665827bd09bSSatish Balay /* catch all ... not good if we get here */ 666827bd09bSSatish Balay return(NULL); 667827bd09bSSatish Balay } 668827bd09bSSatish Balay 669827bd09bSSatish Balay 670827bd09bSSatish Balay 671827bd09bSSatish Balay 672827bd09bSSatish Balay 673827bd09bSSatish Balay 674