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