1 #ifndef lint 2 static char vcid[] = "$Id: baij2.c,v 1.1 1996/04/30 18:51:06 balay Exp balay $"; 3 #endif 4 5 #include "baij.h" 6 #include "petsc.h" 7 8 9 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 10 { 11 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 12 int row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 13 int start, end, *ai, *aj; 14 char *table; 15 16 m = a->mbs; 17 ai = a->i; 18 aj = a->j; 19 20 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqBAIJ: illegal overlap value used"); 21 22 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 23 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 24 25 for ( i=0; i<is_max; i++ ) { 26 /* Initialise the two local arrays */ 27 isz = 0; 28 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 29 30 /* Extract the indices, assume there can be duplicate entries */ 31 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 32 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 33 34 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 35 for ( j=0; j<n ; ++j){ 36 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 37 } 38 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 39 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 40 41 k = 0; 42 for ( j=0; j<ov; j++){ /* for each overlap*/ 43 n = isz; 44 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 45 row = nidx[k]; 46 start = ai[row]; 47 end = ai[row+1]; 48 for ( l = start; l<end ; l++){ 49 val = aj[l]; 50 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 51 } 52 } 53 } 54 ierr = ISCreateSeq(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 55 } 56 PetscFree(table); 57 PetscFree(nidx); 58 return 0; 59 } 60