1 2 /* 3 Routines to compute overlapping regions of a parallel MPI matrix. 4 Used for finding submatrices that were shared across processors. 5 */ 6 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h> 7 #include <petscbt.h> 8 9 static PetscErrorCode MatIncreaseOverlap_MPISBAIJ_Once(Mat,PetscInt,IS*); 10 static PetscErrorCode MatIncreaseOverlap_MPISBAIJ_Local(Mat,PetscInt*,PetscInt,PetscInt*,PetscBT*); 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ" 14 PetscErrorCode MatIncreaseOverlap_MPISBAIJ(Mat C,PetscInt is_max,IS is[],PetscInt ov) 15 { 16 PetscErrorCode ierr; 17 PetscInt i,N=C->cmap->N, bs=C->rmap->bs,M=C->rmap->N,Mbs=M/bs,*nidx,isz,iov; 18 IS *is_new,*is_row; 19 Mat *submats; 20 Mat_MPISBAIJ *c=(Mat_MPISBAIJ*)C->data; 21 Mat_SeqSBAIJ *asub_i; 22 PetscBT table; 23 PetscInt *ai,brow,nz,nis,l,nmax,nstages_local,nstages,max_no,pos; 24 const PetscInt *idx; 25 26 PetscFunctionBegin; 27 ierr = PetscMalloc(is_max*sizeof(IS),&is_new);CHKERRQ(ierr); 28 /* Convert the indices into block format */ 29 ierr = ISCompressIndicesGeneral(N,bs,is_max,is,is_new);CHKERRQ(ierr); 30 if (ov < 0){ SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 31 32 /* ----- previous non-scalable implementation ----- */ 33 PetscBool flg=PETSC_FALSE; 34 ierr = PetscOptionsHasName(PETSC_NULL, "-IncreaseOverlap_old", &flg);CHKERRQ(ierr); 35 if (flg){ /* previous non-scalable implementation */ 36 printf("use previous non-scalable implementation...\n"); 37 for (i=0; i<ov; ++i) { 38 ierr = MatIncreaseOverlap_MPISBAIJ_Once(C,is_max,is_new);CHKERRQ(ierr); 39 } 40 } else { /* scalable implementation using modified BAIJ routines */ 41 42 ierr = PetscMalloc((Mbs+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 43 ierr = PetscBTCreate(Mbs,table);CHKERRQ(ierr); /* for column search */ 44 45 /* Create is_row */ 46 ierr = PetscMalloc(is_max*sizeof(IS **),&is_row);CHKERRQ(ierr); 47 ierr = ISCreateStride(PETSC_COMM_SELF,Mbs,0,1,&is_row[0]);CHKERRQ(ierr); 48 for (i=1; i<is_max; i++) is_row[i] = is_row[0]; /* reuse is_row[0] */ 49 50 /* Allocate memory to hold all the submatrices - Modified from MatGetSubMatrices_MPIBAIJ() */ 51 ierr = PetscMalloc((is_max+1)*sizeof(Mat),&submats);CHKERRQ(ierr); 52 53 /* Determine the number of stages through which submatrices are done */ 54 nmax = 20*1000000 / (c->Nbs * sizeof(PetscInt)); 55 if (!nmax) nmax = 1; 56 nstages_local = is_max/nmax + ((is_max % nmax)?1:0); 57 58 /* Make sure every processor loops through the nstages */ 59 ierr = MPI_Allreduce(&nstages_local,&nstages,1,MPIU_INT,MPI_MAX,((PetscObject)C)->comm);CHKERRQ(ierr); 60 61 for (iov=0; iov<ov; ++iov) { 62 /* 1) Get submats for column search */ 63 for (i=0,pos=0; i<nstages; i++) { 64 if (pos+nmax <= is_max) max_no = nmax; 65 else if (pos == is_max) max_no = 0; 66 else max_no = is_max-pos; 67 68 c->ijonly = PETSC_TRUE; 69 ierr = MatGetSubMatrices_MPIBAIJ_local(C,max_no,is_row+pos,is_new+pos,MAT_INITIAL_MATRIX,submats+pos);CHKERRQ(ierr); 70 pos += max_no; 71 } 72 73 /* 2) Row search */ 74 ierr = MatIncreaseOverlap_MPIBAIJ_Once(C,is_max,is_new);CHKERRQ(ierr); 75 76 /* 3) Column search */ 77 for (i=0; i<is_max; i++){ 78 asub_i = (Mat_SeqSBAIJ*)submats[i]->data; 79 ai=asub_i->i;; 80 81 /* put is_new obtained from MatIncreaseOverlap_MPIBAIJ() to table */ 82 ierr = PetscBTMemzero(Mbs,table);CHKERRQ(ierr); 83 84 ierr = ISGetIndices(is_new[i],&idx);CHKERRQ(ierr); 85 ierr = ISGetLocalSize(is_new[i],&nis);CHKERRQ(ierr); 86 for (l=0; l<nis; l++) { 87 ierr = PetscBTSet(table,idx[l]);CHKERRQ(ierr); 88 nidx[l] = idx[l]; 89 } 90 isz = nis; 91 92 /* add column entries to table */ 93 for (brow=0; brow<Mbs; brow++){ 94 nz = ai[brow+1] - ai[brow]; 95 if (nz) { 96 if (!PetscBTLookupSet(table,brow)) nidx[isz++] = brow; 97 } 98 } 99 ierr = ISRestoreIndices(is_new[i],&idx);CHKERRQ(ierr); 100 ierr = ISDestroy(is_new[i]);CHKERRQ(ierr); 101 102 /* create updated is_new */ 103 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,is_new+i);CHKERRQ(ierr); 104 } 105 106 for (i=0; i<is_max; i++){ierr = MatDestroy(submats[i]);CHKERRQ(ierr);} 107 } 108 109 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 110 ierr = PetscFree(submats);CHKERRQ(ierr); 111 ierr = ISDestroy(is_row[0]);CHKERRQ(ierr); 112 ierr = PetscFree(is_row);CHKERRQ(ierr); 113 ierr = PetscFree(nidx);CHKERRQ(ierr); 114 115 } 116 //--------end of new---------- 117 118 for (i=0; i<is_max; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 119 ierr = ISExpandIndicesGeneral(N,bs,is_max,is_new,is);CHKERRQ(ierr); 120 121 for (i=0; i<is_max; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 122 ierr = PetscFree(is_new);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 typedef enum {MINE,OTHER} WhoseOwner; 127 /* data1, odata1 and odata2 are packed in the format (for communication): 128 data[0] = is_max, no of is 129 data[1] = size of is[0] 130 ... 131 data[is_max] = size of is[is_max-1] 132 data[is_max + 1] = data(is[0]) 133 ... 134 data[is_max+1+sum(size of is[k]), k=0,...,i-1] = data(is[i]) 135 ... 136 data2 is packed in the format (for creating output is[]): 137 data[0] = is_max, no of is 138 data[1] = size of is[0] 139 ... 140 data[is_max] = size of is[is_max-1] 141 data[is_max + 1] = data(is[0]) 142 ... 143 data[is_max + 1 + Mbs*i) = data(is[i]) 144 ... 145 */ 146 #undef __FUNCT__ 147 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Once" 148 static PetscErrorCode MatIncreaseOverlap_MPISBAIJ_Once(Mat C,PetscInt is_max,IS is[]) 149 { 150 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 151 PetscErrorCode ierr; 152 PetscMPIInt size,rank,tag1,tag2,*len_s,nrqr,nrqs,*id_r1,*len_r1,flag,len; 153 const PetscInt *idx_i; 154 PetscInt idx,isz,col,*n,*data1,**data1_start,*data2,*data2_i,*data,*data_i, 155 Mbs,i,j,k,*odata1,*odata2, 156 proc_id,**odata2_ptr,*ctable=0,*btable,len_max,len_est; 157 PetscInt proc_end=0,*iwork,len_unused,nodata2; 158 PetscInt ois_max; /* max no of is[] in each of processor */ 159 char *t_p; 160 MPI_Comm comm; 161 MPI_Request *s_waits1,*s_waits2,r_req; 162 MPI_Status *s_status,r_status; 163 PetscBT *table; /* mark indices of this processor's is[] */ 164 PetscBT table_i; 165 PetscBT otable; /* mark indices of other processors' is[] */ 166 PetscInt bs=C->rmap->bs,Bn = c->B->cmap->n,Bnbs = Bn/bs,*Bowners; 167 IS garray_local,garray_gl; 168 169 PetscFunctionBegin; 170 comm = ((PetscObject)C)->comm; 171 size = c->size; 172 rank = c->rank; 173 Mbs = c->Mbs; 174 175 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 176 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 177 178 /* create tables used in 179 step 1: table[i] - mark c->garray of proc [i] 180 step 3: table[i] - mark indices of is[i] when whose=MINE 181 table[0] - mark incideces of is[] when whose=OTHER */ 182 len = PetscMax(is_max, size);CHKERRQ(ierr); 183 ierr = PetscMalloc2(len,PetscBT,&table,(Mbs/PETSC_BITS_PER_BYTE+1)*len,char,&t_p);CHKERRQ(ierr); 184 for (i=0; i<len; i++) { 185 table[i] = t_p + (Mbs/PETSC_BITS_PER_BYTE+1)*i; 186 } 187 188 ierr = MPI_Allreduce(&is_max,&ois_max,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 189 190 /* 1. Send this processor's is[] to other processors */ 191 /*---------------------------------------------------*/ 192 /* allocate spaces */ 193 ierr = PetscMalloc(is_max*sizeof(PetscInt),&n);CHKERRQ(ierr); 194 len = 0; 195 for (i=0; i<is_max; i++) { 196 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 197 len += n[i]; 198 } 199 if (!len) { 200 is_max = 0; 201 } else { 202 len += 1 + is_max; /* max length of data1 for one processor */ 203 } 204 205 206 ierr = PetscMalloc((size*len+1)*sizeof(PetscInt),&data1);CHKERRQ(ierr); 207 ierr = PetscMalloc(size*sizeof(PetscInt*),&data1_start);CHKERRQ(ierr); 208 for (i=0; i<size; i++) data1_start[i] = data1 + i*len; 209 210 ierr = PetscMalloc4(size,PetscInt,&len_s,size,PetscInt,&btable,size,PetscInt,&iwork,size+1,PetscInt,&Bowners);CHKERRQ(ierr); 211 212 /* gather c->garray from all processors */ 213 ierr = ISCreateGeneral(comm,Bnbs,c->garray,PETSC_COPY_VALUES,&garray_local);CHKERRQ(ierr); 214 ierr = ISAllGather(garray_local, &garray_gl);CHKERRQ(ierr); 215 ierr = ISDestroy(garray_local);CHKERRQ(ierr); 216 ierr = MPI_Allgather(&Bnbs,1,MPIU_INT,Bowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 217 Bowners[0] = 0; 218 for (i=0; i<size; i++) Bowners[i+1] += Bowners[i]; 219 220 if (is_max){ 221 /* hash table ctable which maps c->row to proc_id) */ 222 ierr = PetscMalloc(Mbs*sizeof(PetscInt),&ctable);CHKERRQ(ierr); 223 for (proc_id=0,j=0; proc_id<size; proc_id++) { 224 for (; j<C->rmap->range[proc_id+1]/bs; j++) { 225 ctable[j] = proc_id; 226 } 227 } 228 229 /* hash tables marking c->garray */ 230 ierr = ISGetIndices(garray_gl,&idx_i); 231 for (i=0; i<size; i++){ 232 table_i = table[i]; 233 ierr = PetscBTMemzero(Mbs,table_i);CHKERRQ(ierr); 234 for (j = Bowners[i]; j<Bowners[i+1]; j++){ /* go through B cols of proc[i]*/ 235 ierr = PetscBTSet(table_i,idx_i[j]);CHKERRQ(ierr); 236 } 237 } 238 ierr = ISRestoreIndices(garray_gl,&idx_i);CHKERRQ(ierr); 239 } /* if (is_max) */ 240 ierr = ISDestroy(garray_gl);CHKERRQ(ierr); 241 242 /* evaluate communication - mesg to who, length, and buffer space */ 243 for (i=0; i<size; i++) len_s[i] = 0; 244 245 /* header of data1 */ 246 for (proc_id=0; proc_id<size; proc_id++){ 247 iwork[proc_id] = 0; 248 *data1_start[proc_id] = is_max; 249 data1_start[proc_id]++; 250 for (j=0; j<is_max; j++) { 251 if (proc_id == rank){ 252 *data1_start[proc_id] = n[j]; 253 } else { 254 *data1_start[proc_id] = 0; 255 } 256 data1_start[proc_id]++; 257 } 258 } 259 260 for (i=0; i<is_max; i++) { 261 ierr = ISGetIndices(is[i],&idx_i);CHKERRQ(ierr); 262 for (j=0; j<n[i]; j++){ 263 idx = idx_i[j]; 264 *data1_start[rank] = idx; data1_start[rank]++; /* for local proccessing */ 265 proc_end = ctable[idx]; 266 for (proc_id=0; proc_id<=proc_end; proc_id++){ /* for others to process */ 267 if (proc_id == rank ) continue; /* done before this loop */ 268 if (proc_id < proc_end && !PetscBTLookup(table[proc_id],idx)) 269 continue; /* no need for sending idx to [proc_id] */ 270 *data1_start[proc_id] = idx; data1_start[proc_id]++; 271 len_s[proc_id]++; 272 } 273 } 274 /* update header data */ 275 for (proc_id=0; proc_id<size; proc_id++){ 276 if (proc_id== rank) continue; 277 *(data1 + proc_id*len + 1 + i) = len_s[proc_id] - iwork[proc_id]; 278 iwork[proc_id] = len_s[proc_id] ; 279 } 280 ierr = ISRestoreIndices(is[i],&idx_i);CHKERRQ(ierr); 281 } 282 283 nrqs = 0; nrqr = 0; 284 for (i=0; i<size; i++){ 285 data1_start[i] = data1 + i*len; 286 if (len_s[i]){ 287 nrqs++; 288 len_s[i] += 1 + is_max; /* add no. of header msg */ 289 } 290 } 291 292 for (i=0; i<is_max; i++) { 293 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 294 } 295 ierr = PetscFree(n);CHKERRQ(ierr); 296 ierr = PetscFree(ctable);CHKERRQ(ierr); 297 298 /* Determine the number of messages to expect, their lengths, from from-ids */ 299 ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&nrqr);CHKERRQ(ierr); 300 ierr = PetscGatherMessageLengths(comm,nrqs,nrqr,len_s,&id_r1,&len_r1);CHKERRQ(ierr); 301 302 /* Now post the sends */ 303 ierr = PetscMalloc2(size,MPI_Request,&s_waits1,size,MPI_Request,&s_waits2);CHKERRQ(ierr); 304 k = 0; 305 for (proc_id=0; proc_id<size; proc_id++){ /* send data1 to processor [proc_id] */ 306 if (len_s[proc_id]){ 307 ierr = MPI_Isend(data1_start[proc_id],len_s[proc_id],MPIU_INT,proc_id,tag1,comm,s_waits1+k);CHKERRQ(ierr); 308 k++; 309 } 310 } 311 312 /* 2. Receive other's is[] and process. Then send back */ 313 /*-----------------------------------------------------*/ 314 len = 0; 315 for (i=0; i<nrqr; i++){ 316 if (len_r1[i] > len)len = len_r1[i]; 317 } 318 ierr = PetscFree(len_r1);CHKERRQ(ierr); 319 ierr = PetscFree(id_r1);CHKERRQ(ierr); 320 321 for (proc_id=0; proc_id<size; proc_id++) 322 len_s[proc_id] = iwork[proc_id] = 0; 323 324 ierr = PetscMalloc((len+1)*sizeof(PetscInt),&odata1);CHKERRQ(ierr); 325 ierr = PetscMalloc(size*sizeof(PetscInt**),&odata2_ptr);CHKERRQ(ierr); 326 ierr = PetscBTCreate(Mbs,otable);CHKERRQ(ierr); 327 328 len_max = ois_max*(Mbs+1); /* max space storing all is[] for each receive */ 329 len_est = 2*len_max; /* estimated space of storing is[] for all receiving messages */ 330 ierr = PetscMalloc((len_est+1)*sizeof(PetscInt),&odata2);CHKERRQ(ierr); 331 nodata2 = 0; /* nodata2+1: num of PetscMalloc(,&odata2_ptr[]) called */ 332 odata2_ptr[nodata2] = odata2; 333 len_unused = len_est; /* unused space in the array odata2_ptr[nodata2]-- needs to be >= len_max */ 334 335 k = 0; 336 while (k < nrqr){ 337 /* Receive messages */ 338 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag1,comm,&flag,&r_status);CHKERRQ(ierr); 339 if (flag){ 340 ierr = MPI_Get_count(&r_status,MPIU_INT,&len);CHKERRQ(ierr); 341 proc_id = r_status.MPI_SOURCE; 342 ierr = MPI_Irecv(odata1,len,MPIU_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 343 ierr = MPI_Wait(&r_req,&r_status);CHKERRQ(ierr); 344 345 /* Process messages */ 346 /* make sure there is enough unused space in odata2 array */ 347 if (len_unused < len_max){ /* allocate more space for odata2 */ 348 ierr = PetscMalloc((len_est+1)*sizeof(PetscInt),&odata2);CHKERRQ(ierr); 349 odata2_ptr[++nodata2] = odata2; 350 len_unused = len_est; 351 } 352 353 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,odata1,OTHER,odata2,&otable);CHKERRQ(ierr); 354 len = 1 + odata2[0]; 355 for (i=0; i<odata2[0]; i++){ 356 len += odata2[1 + i]; 357 } 358 359 /* Send messages back */ 360 ierr = MPI_Isend(odata2,len,MPIU_INT,proc_id,tag2,comm,s_waits2+k);CHKERRQ(ierr); 361 k++; 362 odata2 += len; 363 len_unused -= len; 364 len_s[proc_id] = len; /* num of messages sending back to [proc_id] by this proc */ 365 } 366 } 367 ierr = PetscFree(odata1);CHKERRQ(ierr); 368 ierr = PetscBTDestroy(otable);CHKERRQ(ierr); 369 370 /* 3. Do local work on this processor's is[] */ 371 /*-------------------------------------------*/ 372 /* make sure there is enough unused space in odata2(=data) array */ 373 len_max = is_max*(Mbs+1); /* max space storing all is[] for this processor */ 374 if (len_unused < len_max){ /* allocate more space for odata2 */ 375 ierr = PetscMalloc((len_est+1)*sizeof(PetscInt),&odata2);CHKERRQ(ierr); 376 odata2_ptr[++nodata2] = odata2; 377 len_unused = len_est; 378 } 379 380 data = odata2; 381 ierr = MatIncreaseOverlap_MPISBAIJ_Local(C,data1_start[rank],MINE,data,table);CHKERRQ(ierr); 382 ierr = PetscFree(data1_start);CHKERRQ(ierr); 383 384 /* 4. Receive work done on other processors, then merge */ 385 /*------------------------------------------------------*/ 386 /* get max number of messages that this processor expects to recv */ 387 ierr = MPI_Allreduce(len_s,iwork,size,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 388 ierr = PetscMalloc((iwork[rank]+1)*sizeof(PetscInt),&data2);CHKERRQ(ierr); 389 ierr = PetscFree4(len_s,btable,iwork,Bowners);CHKERRQ(ierr); 390 391 k = 0; 392 while (k < nrqs){ 393 /* Receive messages */ 394 ierr = MPI_Iprobe(MPI_ANY_SOURCE,tag2,comm,&flag,&r_status); 395 if (flag){ 396 ierr = MPI_Get_count(&r_status,MPIU_INT,&len);CHKERRQ(ierr); 397 proc_id = r_status.MPI_SOURCE; 398 ierr = MPI_Irecv(data2,len,MPIU_INT,proc_id,r_status.MPI_TAG,comm,&r_req);CHKERRQ(ierr); 399 ierr = MPI_Wait(&r_req,&r_status);CHKERRQ(ierr); 400 if (len > 1+is_max){ /* Add data2 into data */ 401 data2_i = data2 + 1 + is_max; 402 for (i=0; i<is_max; i++){ 403 table_i = table[i]; 404 data_i = data + 1 + is_max + Mbs*i; 405 isz = data[1+i]; 406 for (j=0; j<data2[1+i]; j++){ 407 col = data2_i[j]; 408 if (!PetscBTLookupSet(table_i,col)) {data_i[isz++] = col;} 409 } 410 data[1+i] = isz; 411 if (i < is_max - 1) data2_i += data2[1+i]; 412 } 413 } 414 k++; 415 } 416 } 417 ierr = PetscFree(data2);CHKERRQ(ierr); 418 ierr = PetscFree2(table,t_p);CHKERRQ(ierr); 419 420 /* phase 1 sends are complete */ 421 ierr = PetscMalloc(size*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 422 if (nrqs) {ierr = MPI_Waitall(nrqs,s_waits1,s_status);CHKERRQ(ierr);} 423 ierr = PetscFree(data1);CHKERRQ(ierr); 424 425 /* phase 2 sends are complete */ 426 if (nrqr){ierr = MPI_Waitall(nrqr,s_waits2,s_status);CHKERRQ(ierr);} 427 ierr = PetscFree2(s_waits1,s_waits2);CHKERRQ(ierr); 428 ierr = PetscFree(s_status);CHKERRQ(ierr); 429 430 /* 5. Create new is[] */ 431 /*--------------------*/ 432 for (i=0; i<is_max; i++) { 433 data_i = data + 1 + is_max + Mbs*i; 434 ierr = ISCreateGeneral(PETSC_COMM_SELF,data[1+i],data_i,PETSC_COPY_VALUES,is+i);CHKERRQ(ierr); 435 } 436 for (k=0; k<=nodata2; k++){ 437 ierr = PetscFree(odata2_ptr[k]);CHKERRQ(ierr); 438 } 439 ierr = PetscFree(odata2_ptr);CHKERRQ(ierr); 440 441 PetscFunctionReturn(0); 442 } 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Local" 446 /* 447 MatIncreaseOverlap_MPISBAIJ_Local - Called by MatIncreaseOverlap, to do 448 the work on the local processor. 449 450 Inputs: 451 C - MAT_MPISBAIJ; 452 data - holds is[]. See MatIncreaseOverlap_MPISBAIJ_Once() for the format. 453 whose - whose is[] to be processed, 454 MINE: this processor's is[] 455 OTHER: other processor's is[] 456 Output: 457 nidx - whose = MINE: 458 holds input and newly found indices in the same format as data 459 whose = OTHER: 460 only holds the newly found indices 461 table - table[i]: mark the indices of is[i], i=0,...,is_max. Used only in the case 'whose=MINE'. 462 */ 463 /* Would computation be reduced by swapping the loop 'for each is' and 'for each row'? */ 464 static PetscErrorCode MatIncreaseOverlap_MPISBAIJ_Local(Mat C,PetscInt *data,PetscInt whose,PetscInt *nidx,PetscBT *table) 465 { 466 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 467 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)(c->A)->data; 468 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(c->B)->data; 469 PetscErrorCode ierr; 470 PetscInt row,mbs,Mbs,*nidx_i,col,col_max,isz,isz0,*ai,*aj,*bi,*bj,*garray,rstart,l; 471 PetscInt a_start,a_end,b_start,b_end,i,j,k,is_max,*idx_i,n; 472 PetscBT table0; /* mark the indices of input is[] for look up */ 473 PetscBT table_i; /* poits to i-th table. When whose=OTHER, a single table is used for all is[] */ 474 475 PetscFunctionBegin; 476 Mbs = c->Mbs; mbs = a->mbs; 477 ai = a->i; aj = a->j; 478 bi = b->i; bj = b->j; 479 garray = c->garray; 480 rstart = c->rstartbs; 481 is_max = data[0]; 482 483 ierr = PetscBTCreate(Mbs,table0);CHKERRQ(ierr); 484 485 nidx[0] = is_max; 486 idx_i = data + is_max + 1; /* ptr to input is[0] array */ 487 nidx_i = nidx + is_max + 1; /* ptr to output is[0] array */ 488 for (i=0; i<is_max; i++) { /* for each is */ 489 isz = 0; 490 n = data[1+i]; /* size of input is[i] */ 491 492 /* initialize and set table_i(mark idx and nidx) and table0(only mark idx) */ 493 if (whose == MINE){ /* process this processor's is[] */ 494 table_i = table[i]; 495 nidx_i = nidx + 1+ is_max + Mbs*i; 496 } else { /* process other processor's is[] - only use one temp table */ 497 table_i = table[0]; 498 } 499 ierr = PetscBTMemzero(Mbs,table_i);CHKERRQ(ierr); 500 ierr = PetscBTMemzero(Mbs,table0);CHKERRQ(ierr); 501 if (n==0) { 502 nidx[1+i] = 0; /* size of new is[i] */ 503 continue; 504 } 505 506 isz0 = 0; col_max = 0; 507 for (j=0; j<n; j++){ 508 col = idx_i[j]; 509 if (col >= Mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"index col %D >= Mbs %D",col,Mbs); 510 if(!PetscBTLookupSet(table_i,col)) { 511 ierr = PetscBTSet(table0,col);CHKERRQ(ierr); 512 if (whose == MINE) {nidx_i[isz0] = col;} 513 if (col_max < col) col_max = col; 514 isz0++; 515 } 516 } 517 518 if (whose == MINE) {isz = isz0;} 519 k = 0; /* no. of indices from input is[i] that have been examined */ 520 for (row=0; row<mbs; row++){ 521 a_start = ai[row]; a_end = ai[row+1]; 522 b_start = bi[row]; b_end = bi[row+1]; 523 if (PetscBTLookup(table0,row+rstart)){ /* row is on input is[i]: 524 do row search: collect all col in this row */ 525 for (l = a_start; l<a_end ; l++){ /* Amat */ 526 col = aj[l] + rstart; 527 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 528 } 529 for (l = b_start; l<b_end ; l++){ /* Bmat */ 530 col = garray[bj[l]]; 531 if (!PetscBTLookupSet(table_i,col)) {nidx_i[isz++] = col;} 532 } 533 k++; 534 if (k >= isz0) break; /* for (row=0; row<mbs; row++) */ 535 } else { /* row is not on input is[i]: 536 do col serach: add row onto nidx_i if there is a col in nidx_i */ 537 for (l = a_start; l<a_end ; l++){ /* Amat */ 538 col = aj[l] + rstart; 539 if (col > col_max) break; 540 if (PetscBTLookup(table0,col)){ 541 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 542 break; /* for l = start; l<end ; l++) */ 543 } 544 } 545 for (l = b_start; l<b_end ; l++){ /* Bmat */ 546 col = garray[bj[l]]; 547 if (col > col_max) break; 548 if (PetscBTLookup(table0,col)){ 549 if (!PetscBTLookupSet(table_i,row+rstart)) {nidx_i[isz++] = row+rstart;} 550 break; /* for l = start; l<end ; l++) */ 551 } 552 } 553 } 554 } 555 556 if (i < is_max - 1){ 557 idx_i += n; /* ptr to input is[i+1] array */ 558 nidx_i += isz; /* ptr to output is[i+1] array */ 559 } 560 nidx[1+i] = isz; /* size of new is[i] */ 561 } /* for each is */ 562 ierr = PetscBTDestroy(table0);CHKERRQ(ierr); 563 564 PetscFunctionReturn(0); 565 } 566 567 568