/**********************************ivec.c************************************** Author: Henry M. Tufo III e-mail: hmt@cs.brown.edu snail-mail: Division of Applied Mathematics Brown University Providence, RI 02912 Last Modification: 6.21.97 ***********************************ivec.c*************************************/ #include <../src/ksp/pc/impls/tfs/tfs.h> /* sorting args ivec.c ivec.c ... */ #define SORT_OPT 6 #define SORT_STACK 50000 /* allocate an address and size stack for sorter(s) */ static void *offset_stack[2*SORT_STACK]; static PetscInt size_stack[SORT_STACK]; /***********************************ivec.c*************************************/ PetscInt *PCTFS_ivec_copy(PetscInt *arg1, PetscInt *arg2, PetscInt n) { while (n--) { *arg1++ = *arg2++; } return(arg1); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_zero(PetscInt *arg1, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = 0; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_set(PetscInt *arg1, PetscInt arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = arg2; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_max(PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_min( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *(arg1) = PetscMin(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_mult( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ *= *arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_add( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ += *arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_lxor( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1=((*arg1 || *arg2) && !(*arg1 && *arg2)); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_xor( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ ^= *arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_or( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ |= *arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_lor( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = (*arg1 || *arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_and( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ &= *arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_land( PetscInt *arg1, PetscInt *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = (*arg1 && *arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_and3( PetscInt *arg1, PetscInt *arg2, PetscInt *arg3, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = (*arg2++ & *arg3++); } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscInt PCTFS_ivec_sum( PetscInt *arg1, PetscInt n) { PetscInt tmp = 0; while (n--) {tmp += *arg1++;} return(tmp); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_ivec_non_uniform(PetscInt *arg1, PetscInt *arg2, PetscInt n, PetscInt *arg3) { PetscInt i, j, type; PetscFunctionBegin; /* LATER: if we're really motivated we can sort and then unsort */ for (i=0;i length of the list is now size + 1 */ size--; /* do until we're done ... return when stack is exhausted */ for (;;) { /* if list is large enough use quicksort partition exchange code */ if (size > SORT_OPT) { /* start up pointer at element 1 and down at size */ pi = ar+1; pj = ar+size; /* find middle element in list and swap w/ element 1 */ SWAP(*(ar+(size>>1)),*pi) /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ /* note ==> pivot_value in index 0 */ if (*pi > *pj) { SWAP(*pi,*pj) } if (*ar > *pj) { SWAP(*ar,*pj) } else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) } /* partition about pivot_value ... */ /* note lists of length 2 are not guaranteed to be sorted */ for (;;) { /* walk up ... and down ... swap if equal to pivot! */ do pi++; while (*pi<*ar); do pj--; while (*pj>*ar); /* if we've crossed we're done */ if (pj= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort() :: STACK EXHAUSTED!!!"); /* push right hand child iff length > 1 */ if ((*top_s = size-((PetscInt) (pi-ar)))) { *(top_a++) = pi; size -= *top_s+2; top_s++; } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ ar = *(--top_a); size = *(--top_s); } } else { /* else sort small list directly then pop another off stack */ /* insertion sort for bottom */ for (pj=ar+1;pj<=ar+size;pj++) { temp = *pj; for (pi=pj-1;pi>=ar;pi--) { if (*pi <= temp) break; *(pi+1)=*pi; } *(pi+1)=temp; } /* check to see if stack is exhausted ==> DONE */ if (top_s==bottom_s) PetscFunctionReturn(0); /* else pop another list from the stack */ ar = *(--top_a); size = *(--top_s); } } PetscFunctionReturn(0); } /******************************************************************************/ PetscErrorCode PCTFS_ivec_sort_companion( PetscInt *ar, PetscInt *ar2, PetscInt size) { PetscInt *pi, *pj, temp, temp2; PetscInt **top_a = (PetscInt **)offset_stack; PetscInt *top_s = size_stack, *bottom_s = size_stack; PetscInt *pi2, *pj2; PetscInt mid; PetscFunctionBegin; /* we're really interested in the offset of the last element */ /* ==> length of the list is now size + 1 */ size--; /* do until we're done ... return when stack is exhausted */ for (;;) { /* if list is large enough use quicksort partition exchange code */ if (size > SORT_OPT) { /* start up pointer at element 1 and down at size */ mid = size>>1; pi = ar+1; pj = ar+mid; pi2 = ar2+1; pj2 = ar2+mid; /* find middle element in list and swap w/ element 1 */ SWAP(*pi,*pj) SWAP(*pi2,*pj2) /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ /* note ==> pivot_value in index 0 */ pj = ar+size; pj2 = ar2+size; if (*pi > *pj) { SWAP(*pi,*pj) SWAP(*pi2,*pj2) } if (*ar > *pj) { SWAP(*ar,*pj) SWAP(*ar2,*pj2) } else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) SWAP(*(ar2),*(ar2+1)) } /* partition about pivot_value ... */ /* note lists of length 2 are not guaranteed to be sorted */ for (;;) { /* walk up ... and down ... swap if equal to pivot! */ do { pi++; pi2++; } while (*pi<*ar); do { pj--; pj2--; } while (*pj>*ar); /* if we've crossed we're done */ if (pj= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion() :: STACK EXHAUSTED!!!"); /* push right hand child iff length > 1 */ if ((*top_s = size-((PetscInt) (pi-ar)))) { *(top_a++) = pi; *(top_a++) = pi2; size -= *top_s+2; top_s++; } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ ar2 = *(--top_a); ar = *(--top_a); size = *(--top_s); } } else { /* else sort small list directly then pop another off stack */ /* insertion sort for bottom */ for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { temp = *pj; temp2 = *pj2; for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { if (*pi <= temp) break; *(pi+1)=*pi; *(pi2+1)=*pi2; } *(pi+1)=temp; *(pi2+1)=temp2; } /* check to see if stack is exhausted ==> DONE */ if (top_s==bottom_s) PetscFunctionReturn(0); /* else pop another list from the stack */ ar2 = *(--top_a); ar = *(--top_a); size = *(--top_s); } } PetscFunctionReturn(0); } /******************************************************************************/ PetscErrorCode PCTFS_ivec_sort_companion_hack( PetscInt *ar, PetscInt **ar2, PetscInt size) { PetscInt *pi, *pj, temp, *ptr; PetscInt **top_a = (PetscInt **)offset_stack; PetscInt *top_s = size_stack, *bottom_s = size_stack; PetscInt **pi2, **pj2; PetscInt mid; PetscFunctionBegin; /* we're really interested in the offset of the last element */ /* ==> length of the list is now size + 1 */ size--; /* do until we're done ... return when stack is exhausted */ for (;;) { /* if list is large enough use quicksort partition exchange code */ if (size > SORT_OPT) { /* start up pointer at element 1 and down at size */ mid = size>>1; pi = ar+1; pj = ar+mid; pi2 = ar2+1; pj2 = ar2+mid; /* find middle element in list and swap w/ element 1 */ SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) /* order element 0,1,size-1 st {M,L,...,U} w/L<=M<=U */ /* note ==> pivot_value in index 0 */ pj = ar+size; pj2 = ar2+size; if (*pi > *pj) { SWAP(*pi,*pj) P_SWAP(*pi2,*pj2) } if (*ar > *pj) { SWAP(*ar,*pj) P_SWAP(*ar2,*pj2) } else if (*pi > *ar) { SWAP(*(ar),*(ar+1)) P_SWAP(*(ar2),*(ar2+1)) } /* partition about pivot_value ... */ /* note lists of length 2 are not guaranteed to be sorted */ for (;;) { /* walk up ... and down ... swap if equal to pivot! */ do {pi++; pi2++;} while (*pi<*ar); do {pj--; pj2--;} while (*pj>*ar); /* if we've crossed we're done */ if (pj= SORT_STACK) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_ivec_sort_companion_hack() :: STACK EXHAUSTED!!!"); /* push right hand child iff length > 1 */ if ((*top_s = size-((PetscInt) (pi-ar)))) { *(top_a++) = pi; *(top_a++) = (PetscInt*) pi2; size -= *top_s+2; top_s++; } else if (size -= *top_s+2) {;} /* set up for next loop iff there is something to do */ else { /* might as well pop - note NR_OPT >=2 ==> we're ok! */ ar2 = (PetscInt **) *(--top_a); ar = *(--top_a); size = *(--top_s); } } else { /* else sort small list directly then pop another off stack */ /* insertion sort for bottom */ for (pj=ar+1, pj2=ar2+1;pj<=ar+size;pj++,pj2++) { temp = *pj; ptr = *pj2; for (pi=pj-1,pi2=pj2-1;pi>=ar;pi--,pi2--) { if (*pi <= temp) break; *(pi+1)=*pi; *(pi2+1)=*pi2; } *(pi+1)=temp; *(pi2+1)=ptr; } /* check to see if stack is exhausted ==> DONE */ if (top_s==bottom_s) PetscFunctionReturn(0); /* else pop another list from the stack */ ar2 = (PetscInt **)*(--top_a); ar = *(--top_a); size = *(--top_s); } } PetscFunctionReturn(0); } /******************************************************************************/ PetscErrorCode PCTFS_SMI_sort(void *ar1, void *ar2, PetscInt size, PetscInt type) { PetscFunctionBegin; if (type == SORT_INTEGER) { if (ar2) PCTFS_ivec_sort_companion((PetscInt*)ar1,(PetscInt*)ar2,size); else PCTFS_ivec_sort((PetscInt*)ar1,size); } else if (type == SORT_INT_PTR) { if (ar2) PCTFS_ivec_sort_companion_hack((PetscInt*)ar1,(PetscInt **)ar2,size); else PCTFS_ivec_sort((PetscInt*)ar1,size); } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCTFS_SMI_sort only does SORT_INTEGER!"); PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscInt PCTFS_ivec_linear_search( PetscInt item, PetscInt *list, PetscInt n) { PetscInt tmp = n-1; PetscFunctionBegin; while (n--) { if (*list++ == item) { return(tmp-n); } } return(-1); } /***********************************ivec.c*************************************/ PetscInt PCTFS_ivec_binary_search( PetscInt item, PetscInt *list, PetscInt rh) { PetscInt mid, lh=0; rh--; while (lh<=rh) { mid = (lh+rh)>>1; if (*(list+mid) == item) { return(mid); } if (*(list+mid) > item) { rh = mid-1; } else { lh = mid+1; } } return(-1); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_copy( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = *arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_zero( PetscScalar *arg1, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = 0.0; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_one( PetscScalar *arg1, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = 1.0; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_set( PetscScalar *arg1, PetscScalar arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ = arg2; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_scale( PetscScalar *arg1, PetscScalar arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ *= arg2; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_add( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ += *arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_mult( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1++ *= *arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_max( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = PetscMax(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_max_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = MAX_FABS(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_min( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = PetscMin(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_min_abs( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = MIN_FABS(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /*********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_exists( PetscScalar *arg1, PetscScalar *arg2, PetscInt n) { PetscFunctionBegin; while (n--) { *arg1 = EXISTS(*arg1,*arg2); arg1++; arg2++; } PetscFunctionReturn(0); } /***********************************ivec.c*************************************/ PetscErrorCode PCTFS_rvec_non_uniform(PetscScalar *arg1, PetscScalar *arg2, PetscInt n, PetscInt *arg3) { PetscInt i, j, type; PetscFunctionBegin; /* LATER: if we're really motivated we can sort and then unsort */ for (i=0;i