1 /*$Id: sbaijov.c,v 1.65 2001/08/06 21:15:42 bsmith Exp $*/ 2 3 /* 4 Routines to compute overlapping regions of a parallel MPI matrix. 5 Used for finding submatrices that were shared across processors. 6 */ 7 #include "src/mat/impls/sbaij/mpi/mpisbaij.h" 8 #include "petscbt.h" 9 10 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat,int,IS *); 11 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat,int,char **,int*,int**); 12 static int MatIncreaseOverlap_MPISBAIJ_Receive(Mat,int,int **,int**,int*); 13 14 /* this function is sasme as MatCompressIndicesGeneral_MPIBAIJ -- should be removed! */ 15 #undef __FUNCT__ 16 #define __FUNCT__ "MatCompressIndicesGeneral_MPISBAIJ" 17 static int MatCompressIndicesGeneral_MPISBAIJ(Mat C,int imax,const IS is_in[],IS is_out[]) 18 { 19 Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)C->data; 20 int ierr,isz,bs = baij->bs,n,i,j,*idx,ival; 21 #if defined (PETSC_USE_CTABLE) 22 PetscTable gid1_lid1; 23 int tt, gid1, *nidx; 24 PetscTablePosition tpos; 25 #else 26 int Nbs,*nidx; 27 PetscBT table; 28 #endif 29 30 PetscFunctionBegin; 31 /* printf(" ...MatCompressIndicesGeneral_MPISBAIJ is called ...\n"); */ 32 #if defined (PETSC_USE_CTABLE) 33 ierr = PetscTableCreate(baij->mbs,&gid1_lid1);CHKERRQ(ierr); 34 #else 35 Nbs = baij->Nbs; 36 ierr = PetscMalloc((Nbs+1)*sizeof(int),&nidx);CHKERRQ(ierr); 37 ierr = PetscBTCreate(Nbs,table);CHKERRQ(ierr); 38 #endif 39 for (i=0; i<imax; i++) { 40 isz = 0; 41 #if defined (PETSC_USE_CTABLE) 42 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 43 #else 44 ierr = PetscBTMemzero(Nbs,table);CHKERRQ(ierr); 45 #endif 46 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 47 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 48 for (j=0; j<n ; j++) { 49 ival = idx[j]/bs; /* convert the indices into block indices */ 50 #if defined (PETSC_USE_CTABLE) 51 ierr = PetscTableFind(gid1_lid1,ival+1,&tt);CHKERRQ(ierr); 52 if (!tt) { 53 ierr = PetscTableAdd(gid1_lid1,ival+1,isz+1);CHKERRQ(ierr); 54 isz++; 55 } 56 #else 57 if (ival>Nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than mat-dim"); 58 if(!PetscBTLookupSet(table,ival)) { nidx[isz++] = ival;} 59 #endif 60 } 61 ierr = ISRestoreIndices(is_in[i],&idx);CHKERRQ(ierr); 62 #if defined (PETSC_USE_CTABLE) 63 ierr = PetscMalloc((isz+1)*sizeof(int),&nidx);CHKERRQ(ierr); 64 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 65 j = 0; 66 while (tpos) { 67 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid1,&tt);CHKERRQ(ierr); 68 if (tt-- > isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"index greater than array-dim"); } 69 nidx[tt] = gid1 - 1; 70 j++; 71 } 72 if (j != isz) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"table error: jj != isz"); } 73 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));CHKERRQ(ierr); 74 ierr = PetscFree(nidx);CHKERRQ(ierr); 75 #else 76 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is_out+i));CHKERRQ(ierr); 77 #endif 78 } 79 #if defined (PETSC_USE_CTABLE) 80 ierr = PetscTableDelete(gid1_lid1);CHKERRQ(ierr); 81 #else 82 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 83 ierr = PetscFree(nidx);CHKERRQ(ierr); 84 #endif 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "MatCompressIndicesSorted_MPISBAIJ" 90 static int MatCompressIndicesSorted_MPISBAIJ(Mat C,int imax,const IS is_in[],IS is_out[]) 91 { 92 Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)C->data; 93 int ierr,bs=baij->bs,i,j,k,val,n,*idx,*nidx,*idx_local; 94 PetscTruth flg; 95 #if defined (PETSC_USE_CTABLE) 96 int maxsz; 97 #else 98 int Nbs=baij->Nbs; 99 #endif 100 PetscFunctionBegin; 101 printf(" ... MatCompressIndicesSorted_MPISBAIJ is called ...\n"); 102 for (i=0; i<imax; i++) { 103 ierr = ISSorted(is_in[i],&flg);CHKERRQ(ierr); 104 if (!flg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Indices are not sorted"); 105 } 106 #if defined (PETSC_USE_CTABLE) 107 /* Now check max size */ 108 for (i=0,maxsz=0; i<imax; i++) { 109 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 110 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 111 if (n%bs !=0) SETERRQ(1,"Indices are not block ordered"); 112 n = n/bs; /* The reduced index size */ 113 if (n > maxsz) maxsz = n; 114 } 115 ierr = PetscMalloc((maxsz+1)*sizeof(int),&nidx);CHKERRQ(ierr); 116 #else 117 ierr = PetscMalloc((Nbs+1)*sizeof(int),&nidx);CHKERRQ(ierr); 118 #endif 119 /* Now check if the indices are in block order */ 120 for (i=0; i<imax; i++) { 121 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 122 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 123 if (n%bs !=0) SETERRQ(1,"Indices are not block ordered"); 124 125 n = n/bs; /* The reduced index size */ 126 idx_local = idx; 127 for (j=0; j<n ; j++) { 128 val = idx_local[0]; 129 if (val%bs != 0) SETERRQ(1,"Indices are not block ordered"); 130 for (k=0; k<bs; k++) { 131 if (val+k != idx_local[k]) SETERRQ(1,"Indices are not block ordered"); 132 } 133 nidx[j] = val/bs; 134 idx_local +=bs; 135 } 136 ierr = ISRestoreIndices(is_in[i],&idx);CHKERRQ(ierr); 137 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,nidx,(is_out+i));CHKERRQ(ierr); 138 } 139 ierr = PetscFree(nidx);CHKERRQ(ierr); 140 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "MatExpandIndices_MPISBAIJ" 146 static int MatExpandIndices_MPISBAIJ(Mat C,int imax,const IS is_in[],IS is_out[]) 147 { 148 Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)C->data; 149 int ierr,bs = baij->bs,n,i,j,k,*idx,*nidx; 150 #if defined (PETSC_USE_CTABLE) 151 int maxsz; 152 #else 153 int Nbs = baij->Nbs; 154 #endif 155 156 PetscFunctionBegin; 157 /* printf(" ... MatExpandIndices_MPISBAIJ is called ...\n"); */ 158 #if defined (PETSC_USE_CTABLE) 159 /* Now check max size */ 160 for (i=0,maxsz=0; i<imax; i++) { 161 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 162 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 163 if (n*bs > maxsz) maxsz = n*bs; 164 } 165 ierr = PetscMalloc((maxsz+1)*sizeof(int),&nidx);CHKERRQ(ierr); 166 #else 167 ierr = PetscMalloc((Nbs*bs+1)*sizeof(int),&nidx);CHKERRQ(ierr); 168 #endif 169 170 for (i=0; i<imax; i++) { 171 ierr = ISGetIndices(is_in[i],&idx);CHKERRQ(ierr); 172 ierr = ISGetLocalSize(is_in[i],&n);CHKERRQ(ierr); 173 for (j=0; j<n ; ++j){ 174 for (k=0; k<bs; k++) 175 nidx[j*bs+k] = idx[j]*bs+k; 176 } 177 ierr = ISRestoreIndices(is_in[i],&idx);CHKERRQ(ierr); 178 ierr = ISCreateGeneral(PETSC_COMM_SELF,n*bs,nidx,is_out+i);CHKERRQ(ierr); 179 } 180 ierr = PetscFree(nidx);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 185 #undef __FUNCT__ 186 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ" 187 int MatIncreaseOverlap_MPISBAIJ(Mat C,int imax,IS is[],int ov) 188 { 189 int i,ierr; 190 IS *is_new; 191 192 PetscFunctionBegin; 193 ierr = PetscMalloc(imax*sizeof(IS),&is_new);CHKERRQ(ierr); 194 /* Convert the indices into block format */ 195 ierr = MatCompressIndicesGeneral_MPISBAIJ(C,imax,is,is_new);CHKERRQ(ierr); 196 if (ov < 0){ SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 197 for (i=0; i<ov; ++i) { 198 ierr = MatIncreaseOverlap_MPISBAIJ_Once(C,imax,is_new);CHKERRQ(ierr); 199 } 200 for (i=0; i<imax; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 201 ierr = MatExpandIndices_MPISBAIJ(C,imax,is_new,is);CHKERRQ(ierr); 202 for (i=0; i<imax; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 203 ierr = PetscFree(is_new);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Once" 209 static int MatIncreaseOverlap_MPISBAIJ_Once(Mat C,int imax,IS is[]) 210 { 211 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 212 int **idx,*n,len,*idx_i; 213 int size,rank,Mbs,i,j,k,ierr,**rbuf,row,proc,nrqs,msz,*outdat,*indat; 214 int *onodes1,*olengths1,tag1,tag2,*onodes2,*olengths2,flag,proc_id; 215 PetscBT *table; 216 MPI_Comm comm; 217 MPI_Request *s_waits1,*r_waits1,*s_waits2,*r_waits2; 218 MPI_Status *s_status,r_status; 219 220 PetscFunctionBegin; 221 222 comm = C->comm; 223 size = c->size; 224 rank = c->rank; 225 Mbs = c->Mbs; 226 227 /* 1. Send is[] to all other processors */ 228 /*--------------------------------------*/ 229 /* This processor sends its is[] to all other processors in the format: 230 outdat[0] = is_max, no of is in this processor 231 outdat[1] = n[0], size of is[0] 232 ... 233 outdat[is_max] = n[is_max-1], size of is[is_max-1] 234 outdat[is_max + 1] = data(is[0]) 235 ... 236 outdat[is_max + i] = data(is[i]) 237 ... 238 */ 239 len = (imax+1)*sizeof(int*)+ (imax)*sizeof(int); 240 ierr = PetscMalloc(len,&idx);CHKERRQ(ierr); 241 n = (int*)(idx + imax); 242 243 /* Allocate Memory for outgoing messages */ 244 len = 1 + imax; 245 for (i=0; i<imax; i++) { 246 ierr = ISGetIndices(is[i],&idx[i]);CHKERRQ(ierr); 247 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 248 len += n[i]; 249 } 250 ierr = PetscMalloc(len*sizeof(int),&outdat);CHKERRQ(ierr); 251 252 /* Form the outgoing messages */ 253 outdat[0] = imax; 254 for (i=0; i<imax; i++) { 255 outdat[i+1] = n[i]; 256 } 257 k = imax + 1; 258 for (i=0; i<imax; i++) { /* for is[i] */ 259 idx_i = idx[i]; 260 for (j=0; j<n[i]; j++){ 261 outdat[k] = *(idx_i); 262 /* if (!rank) printf(" outdat[%d] = %d\n",k,outdat[k] ); */ 263 k++; idx_i++; 264 } 265 /* printf(" [%d] n[%d]=%d, k: %d, \n",rank,i,n[i],k); */ 266 } 267 if (k != len) SETERRQ3(1,"[%d] Error on forming the outgoing messages: k %d != len %d",rank,k,len); 268 269 /* Now post the sends */ 270 ierr = PetscMalloc(size*sizeof(MPI_Request),&s_waits1);CHKERRQ(ierr); 271 272 k = 0; 273 for (i=0; i<size; ++i) { /* send outdat to processor [i] */ 274 if (i != rank){ 275 ierr = MPI_Isend(outdat,len,MPI_INT,i,rank,comm,&s_waits1[k]);CHKERRQ(ierr); 276 printf(" [%d] send %d msg to [%d] \n",rank,len,i); 277 k++; 278 } 279 } 280 281 /* 2. Do local work */ 282 /*------------------*/ 283 284 /* No longer need the original indices*/ 285 for (i=0; i<imax; ++i) { 286 ierr = ISRestoreIndices(is[i],idx+i);CHKERRQ(ierr); 287 } 288 ierr = PetscFree(idx);CHKERRQ(ierr); 289 290 /* 3. Receive other's is[] and process. Then send back */ 291 /*----------------------------------------------------*/ 292 /* Send is done */ 293 nrqs = size-1; 294 ierr = PetscMalloc(size*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 295 ierr = MPI_Waitall(nrqs,s_waits1,s_status);CHKERRQ(ierr); 296 ierr = PetscFree(outdat);CHKERRQ(ierr); 297 298 ierr = PetscMalloc(size*sizeof(MPI_Request),&r_waits1);CHKERRQ(ierr); 299 ierr = PetscMalloc(size*sizeof(MPI_Request),&s_waits2);CHKERRQ(ierr); 300 k = 0; 301 do { 302 /* Receive messages */ 303 ierr = MPI_Iprobe(MPI_ANY_SOURCE,MPI_ANY_TAG,comm,&flag,&r_status); 304 if (flag){ 305 ierr = MPI_Get_count(&r_status,MPI_INT,&len); 306 proc_id = r_status.MPI_SOURCE; 307 ierr = ierr = PetscMalloc(len*sizeof(int),&indat);CHKERRQ(ierr); 308 ierr = MPI_Irecv(indat,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_waits1[k]); 309 printf(" [%d] recv %d msg from [%d] \n",rank,len,proc_id); 310 311 /* Process messages -- not done yet */ 312 len = indat[0]; 313 ierr = PetscMalloc(len*sizeof(int),&outdat);CHKERRQ(ierr); 314 for (i=0; i<len; i++){outdat[i] = indat[i+1];} 315 316 /* Send messages back */ 317 printf(" [%d] send %d msg back to [%d] \n",rank,len,proc_id); 318 ierr = MPI_Isend(outdat,len,MPI_INT,proc_id,rank,comm,&s_waits2[k]);CHKERRQ(ierr); 319 320 k++; 321 ierr = PetscFree(outdat);CHKERRQ(ierr); 322 ierr = PetscFree(indat);CHKERRQ(ierr); 323 } 324 } while (k < nrqs); 325 326 /* 4. Receive work done on other processors, then process */ 327 /*--------------------------------------------------------*/ 328 ierr = MPI_Waitall(nrqs,s_waits2,s_status);CHKERRQ(ierr); 329 ierr = PetscMalloc(size*sizeof(MPI_Request),&r_waits2);CHKERRQ(ierr); 330 k = 0; 331 do { 332 /* Receive messages */ 333 ierr = MPI_Iprobe(MPI_ANY_SOURCE,MPI_ANY_TAG,comm,&flag,&r_status); 334 if (flag){ 335 ierr = MPI_Get_count(&r_status,MPI_INT,&len); 336 proc_id = r_status.MPI_SOURCE; 337 ierr = ierr = PetscMalloc(len*sizeof(int),&indat);CHKERRQ(ierr); 338 ierr = MPI_Irecv(indat,len,MPI_INT,proc_id,r_status.MPI_TAG,comm,&r_waits2[k]); 339 printf(" [%d] recv %d msg from [%d] \n",rank,len,proc_id); 340 341 /* Process messages -- not done yet */ 342 343 344 k++; 345 ierr = PetscFree(indat);CHKERRQ(ierr); 346 } 347 } while (k < nrqs); 348 349 #ifdef OLD 350 for (i=0; i<imax; ++i) { 351 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz[i],data[i],is+i);CHKERRQ(ierr); 352 } 353 354 355 ierr = PetscFree(onodes2);CHKERRQ(ierr); 356 ierr = PetscFree(olengths2);CHKERRQ(ierr); 357 358 ierr = PetscFree(pa);CHKERRQ(ierr); 359 ierr = PetscFree(rbuf2);CHKERRQ(ierr); 360 ierr = PetscFree(s_waits1);CHKERRQ(ierr); 361 ierr = PetscFree(r_waits1);CHKERRQ(ierr); 362 ierr = PetscFree(s_waits2);CHKERRQ(ierr); 363 ierr = PetscFree(r_waits2);CHKERRQ(ierr); 364 ierr = PetscFree(table);CHKERRQ(ierr); 365 ierr = PetscFree(s_status);CHKERRQ(ierr); 366 ierr = PetscFree(recv_status);CHKERRQ(ierr); 367 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 368 ierr = PetscFree(xdata);CHKERRQ(ierr); 369 ierr = PetscFree(isz1);CHKERRQ(ierr); 370 #endif /* OLD */ 371 ierr = PetscFree(s_waits1);CHKERRQ(ierr); 372 ierr = PetscFree(r_waits1);CHKERRQ(ierr); 373 ierr = PetscFree(s_waits2);CHKERRQ(ierr); 374 ierr = PetscFree(r_waits2);CHKERRQ(ierr); 375 ierr = PetscFree(s_status);CHKERRQ(ierr); 376 PetscFunctionReturn(0); 377 } 378 379 #undef __FUNCT__ 380 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Local" 381 /* 382 MatIncreaseOverlap_MPISBAIJ_Local - Called by MatincreaseOverlap, to do 383 the work on the local processor. 384 385 Inputs: 386 C - MAT_MPISBAIJ; 387 imax - total no of index sets processed at a time; 388 table - an array of char - size = Mbs bits. 389 390 Output: 391 isz - array containing the count of the solution elements correspondign 392 to each index set; 393 data - pointer to the solutions 394 */ 395 static int MatIncreaseOverlap_MPISBAIJ_Local(Mat C,int imax,PetscBT *table,int *isz,int **data) 396 { 397 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 398 Mat A = c->A,B = c->B; 399 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 400 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)B->data; 401 int start,end,val,max,rstart,cstart,*ai,*aj; 402 int *bi,*bj,*garray,i,j,k,row,*data_i,isz_i; 403 PetscBT table_i; 404 405 PetscFunctionBegin; 406 rstart = c->rstart; 407 cstart = c->cstart; 408 ai = a->i; 409 aj = a->j; 410 bi = b->i; 411 bj = b->j; 412 garray = c->garray; 413 414 415 for (i=0; i<imax; i++) { 416 data_i = data[i]; 417 table_i = table[i]; 418 isz_i = isz[i]; 419 for (j=0,max=isz[i]; j<max; j++) { 420 row = data_i[j] - rstart; 421 start = ai[row]; 422 end = ai[row+1]; 423 for (k=start; k<end; k++) { /* Amat */ 424 val = aj[k] + cstart; 425 if (!PetscBTLookupSet(table_i,val)) { data_i[isz_i++] = val;} 426 } 427 start = bi[row]; 428 end = bi[row+1]; 429 for (k=start; k<end; k++) { /* Bmat */ 430 val = garray[bj[k]]; 431 if (!PetscBTLookupSet(table_i,val)) { data_i[isz_i++] = val;} 432 } 433 } 434 isz[i] = isz_i; 435 } 436 PetscFunctionReturn(0); 437 } 438 #undef __FUNCT__ 439 #define __FUNCT__ "MatIncreaseOverlap_MPISBAIJ_Receive" 440 /* 441 MatIncreaseOverlap_MPISBAIJ_Receive - Process the recieved messages, 442 and return the output 443 444 Input: 445 C - the matrix 446 nrqr - no of messages being processed. 447 rbuf - an array of pointers to the recieved requests 448 449 Output: 450 xdata - array of messages to be sent back 451 isz1 - size of each message 452 453 For better efficiency perhaps we should malloc seperately each xdata[i], 454 then if a remalloc is required we need only copy the data for that one row 455 rather then all previous rows as it is now where a single large chunck of 456 memory is used. 457 458 */ 459 static int MatIncreaseOverlap_MPISBAIJ_Receive(Mat C,int nrqr,int **rbuf,int **xdata,int * isz1) 460 { 461 Mat_MPISBAIJ *c = (Mat_MPISBAIJ*)C->data; 462 Mat A = c->A,B = c->B; 463 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 464 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)B->data; 465 int rstart,cstart,*ai,*aj,*bi,*bj,*garray,i,j,k; 466 int row,total_sz,ct,ct1,ct2,ct3,mem_estimate,oct2,l,start,end; 467 int val,max1,max2,rank,Mbs,no_malloc =0,*tmp,new_estimate,ctr; 468 int *rbuf_i,kmax,rbuf_0,ierr; 469 PetscBT xtable; 470 471 PetscFunctionBegin; 472 rank = c->rank; 473 Mbs = c->Mbs; 474 rstart = c->rstart; 475 cstart = c->cstart; 476 ai = a->i; 477 aj = a->j; 478 bi = b->i; 479 bj = b->j; 480 garray = c->garray; 481 482 483 for (i=0,ct=0,total_sz=0; i<nrqr; ++i) { 484 rbuf_i = rbuf[i]; 485 rbuf_0 = rbuf_i[0]; 486 ct += rbuf_0; 487 for (j=1; j<=rbuf_0; j++) { total_sz += rbuf_i[2*j]; } 488 } 489 490 if (c->Mbs) max1 = ct*(a->nz +b->nz)/c->Mbs; 491 else max1 = 1; 492 mem_estimate = 3*((total_sz > max1 ? total_sz : max1)+1); 493 ierr = PetscMalloc(mem_estimate*sizeof(int),&xdata[0]);CHKERRQ(ierr); 494 ++no_malloc; 495 ierr = PetscBTCreate(Mbs,xtable);CHKERRQ(ierr); 496 ierr = PetscMemzero(isz1,nrqr*sizeof(int));CHKERRQ(ierr); 497 498 ct3 = 0; 499 for (i=0; i<nrqr; i++) { /* for easch mesg from proc i */ 500 rbuf_i = rbuf[i]; 501 rbuf_0 = rbuf_i[0]; 502 ct1 = 2*rbuf_0+1; 503 ct2 = ct1; 504 ct3 += ct1; 505 for (j=1; j<=rbuf_0; j++) { /* for each IS from proc i*/ 506 ierr = PetscBTMemzero(Mbs,xtable);CHKERRQ(ierr); 507 oct2 = ct2; 508 kmax = rbuf_i[2*j]; 509 for (k=0; k<kmax; k++,ct1++) { 510 row = rbuf_i[ct1]; 511 if (!PetscBTLookupSet(xtable,row)) { 512 if (!(ct3 < mem_estimate)) { 513 new_estimate = (int)(1.5*mem_estimate)+1; 514 ierr = PetscMalloc(new_estimate * sizeof(int),&tmp);CHKERRQ(ierr); 515 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(int));CHKERRQ(ierr); 516 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 517 xdata[0] = tmp; 518 mem_estimate = new_estimate; ++no_malloc; 519 for (ctr=1; ctr<=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 520 } 521 xdata[i][ct2++] = row; 522 ct3++; 523 } 524 } 525 for (k=oct2,max2=ct2; k<max2; k++) { 526 row = xdata[i][k] - rstart; 527 start = ai[row]; 528 end = ai[row+1]; 529 for (l=start; l<end; l++) { 530 val = aj[l] + cstart; 531 if (!PetscBTLookupSet(xtable,val)) { 532 if (!(ct3 < mem_estimate)) { 533 new_estimate = (int)(1.5*mem_estimate)+1; 534 ierr = PetscMalloc(new_estimate * sizeof(int),&tmp);CHKERRQ(ierr); 535 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(int));CHKERRQ(ierr); 536 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 537 xdata[0] = tmp; 538 mem_estimate = new_estimate; ++no_malloc; 539 for (ctr=1; ctr<=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 540 } 541 xdata[i][ct2++] = val; 542 ct3++; 543 } 544 } 545 start = bi[row]; 546 end = bi[row+1]; 547 for (l=start; l<end; l++) { 548 val = garray[bj[l]]; 549 if (!PetscBTLookupSet(xtable,val)) { 550 if (!(ct3 < mem_estimate)) { 551 new_estimate = (int)(1.5*mem_estimate)+1; 552 ierr = PetscMalloc(new_estimate * sizeof(int),&tmp);CHKERRQ(ierr); 553 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(int));CHKERRQ(ierr); 554 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 555 xdata[0] = tmp; 556 mem_estimate = new_estimate; ++no_malloc; 557 for (ctr =1; ctr <=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 558 } 559 xdata[i][ct2++] = val; 560 ct3++; 561 } 562 } 563 } 564 /* Update the header*/ 565 xdata[i][2*j] = ct2 - oct2; /* Undo the vector isz1 and use only a var*/ 566 xdata[i][2*j-1] = rbuf_i[2*j-1]; 567 } 568 xdata[i][0] = rbuf_0; 569 xdata[i+1] = xdata[i] + ct2; 570 isz1[i] = ct2; /* size of each message */ 571 } 572 ierr = PetscBTDestroy(xtable);CHKERRQ(ierr); 573 PetscLogInfo(0,"MatIncreaseOverlap_MPISBAIJ:[%d] Allocated %d bytes, required %d, no of mallocs = %d\n",rank,mem_estimate,ct3,no_malloc); 574 PetscFunctionReturn(0); 575 } 576 577 578