1 /* 2 Routines to compute overlapping regions of a parallel MPI matrix 3 and to find submatrices that were shared across processors. 4 */ 5 #include "src/mat/impls/baij/mpi/mpibaij.h" 6 #include "petscbt.h" 7 8 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Once(Mat,PetscInt,IS *); 9 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Local(Mat,PetscInt,char **,PetscInt*,PetscInt**); 10 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Receive(Mat,PetscInt,PetscInt **,PetscInt**,PetscInt*); 11 EXTERN PetscErrorCode MatGetRow_MPIBAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 12 EXTERN PetscErrorCode MatRestoreRow_MPIBAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatIncreaseOverlap_MPIBAIJ" 16 PetscErrorCode MatIncreaseOverlap_MPIBAIJ(Mat C,PetscInt imax,IS is[],PetscInt ov) 17 { 18 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 19 PetscErrorCode ierr; 20 PetscInt i,N=C->N, bs=c->bs; 21 IS *is_new; 22 23 PetscFunctionBegin; 24 ierr = PetscMalloc(imax*sizeof(IS),&is_new);CHKERRQ(ierr); 25 /* Convert the indices into block format */ 26 ierr = ISCompressIndicesGeneral(N,bs,imax,is,is_new);CHKERRQ(ierr); 27 if (ov < 0){ SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified\n");} 28 for (i=0; i<ov; ++i) { 29 ierr = MatIncreaseOverlap_MPIBAIJ_Once(C,imax,is_new);CHKERRQ(ierr); 30 } 31 for (i=0; i<imax; i++) {ierr = ISDestroy(is[i]);CHKERRQ(ierr);} 32 ierr = ISExpandIndicesGeneral(N,bs,imax,is_new,is);CHKERRQ(ierr); 33 for (i=0; i<imax; i++) {ierr = ISDestroy(is_new[i]);CHKERRQ(ierr);} 34 ierr = PetscFree(is_new);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 /* 39 Sample message format: 40 If a processor A wants processor B to process some elements corresponding 41 to index sets is[1], is[5] 42 mesg [0] = 2 (no of index sets in the mesg) 43 ----------- 44 mesg [1] = 1 => is[1] 45 mesg [2] = sizeof(is[1]); 46 ----------- 47 mesg [5] = 5 => is[5] 48 mesg [6] = sizeof(is[5]); 49 ----------- 50 mesg [7] 51 mesg [n] data(is[1]) 52 ----------- 53 mesg[n+1] 54 mesg[m] data(is[5]) 55 ----------- 56 57 Notes: 58 nrqs - no of requests sent (or to be sent out) 59 nrqr - no of requests recieved (which have to be or which have been processed 60 */ 61 #undef __FUNCT__ 62 #define __FUNCT__ "MatIncreaseOverlap_MPIBAIJ_Once" 63 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Once(Mat C,PetscInt imax,IS is[]) 64 { 65 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 66 PetscInt **idx,*n,*w3,*w4,*rtable,**data,len,*idx_i; 67 PetscErrorCode ierr; 68 PetscMPIInt size,rank,tag1,tag2,*w2,*w1,nrqr; 69 PetscInt Mbs,i,j,k,**rbuf,row,proc,nrqs,msz,**outdat,**ptr; 70 PetscInt *ctr,*pa,*tmp,*isz,*isz1,**xdata,**rbuf2; 71 PetscMPIInt *onodes1,*olengths1,*onodes2,*olengths2; 72 PetscBT *table; 73 MPI_Comm comm; 74 MPI_Request *s_waits1,*r_waits1,*s_waits2,*r_waits2; 75 MPI_Status *s_status,*recv_status; 76 77 PetscFunctionBegin; 78 comm = C->comm; 79 size = c->size; 80 rank = c->rank; 81 Mbs = c->Mbs; 82 83 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 84 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 85 86 len = (imax+1)*sizeof(PetscInt*)+ (imax + Mbs)*sizeof(PetscInt); 87 ierr = PetscMalloc(len,&idx);CHKERRQ(ierr); 88 n = (PetscInt*)(idx + imax); 89 rtable = n + imax; 90 91 for (i=0; i<imax; i++) { 92 ierr = ISGetIndices(is[i],&idx[i]);CHKERRQ(ierr); 93 ierr = ISGetLocalSize(is[i],&n[i]);CHKERRQ(ierr); 94 } 95 96 /* Create hash table for the mapping :row -> proc*/ 97 for (i=0,j=0; i<size; i++) { 98 len = c->rowners[i+1]; 99 for (; j<len; j++) { 100 rtable[j] = i; 101 } 102 } 103 104 /* evaluate communication - mesg to who,length of mesg, and buffer space 105 required. Based on this, buffers are allocated, and data copied into them*/ 106 ierr = PetscMalloc(size*2*sizeof(PetscInt)+size*2*sizeof(PetscMPIInt),&w1);CHKERRQ(ierr);/* mesg size */ 107 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 108 w3 = (PetscInt*) (w2 + size); /* no of IS that needs to be sent to proc i */ 109 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 110 ierr = PetscMemzero(w1,size*sizeof(PetscInt)+2*size*sizeof(PetscMPIInt));CHKERRQ(ierr); /* initialise work vector*/ 111 for (i=0; i<imax; i++) { 112 ierr = PetscMemzero(w4,size*sizeof(PetscInt));CHKERRQ(ierr); /* initialise work vector*/ 113 idx_i = idx[i]; 114 len = n[i]; 115 for (j=0; j<len; j++) { 116 row = idx_i[j]; 117 if (row < 0) { 118 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Index set cannot have negative entries"); 119 } 120 proc = rtable[row]; 121 w4[proc]++; 122 } 123 for (j=0; j<size; j++){ 124 if (w4[j]) { w1[j] += w4[j]; w3[j]++;} 125 } 126 } 127 128 nrqs = 0; /* no of outgoing messages */ 129 msz = 0; /* total mesg length (for all proc */ 130 w1[rank] = 0; /* no mesg sent to itself */ 131 w3[rank] = 0; 132 for (i=0; i<size; i++) { 133 if (w1[i]) {w2[i] = 1; nrqs++;} /* there exists a message to proc i */ 134 } 135 /* pa - is list of processors to communicate with */ 136 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt),&pa);CHKERRQ(ierr); 137 for (i=0,j=0; i<size; i++) { 138 if (w1[i]) {pa[j] = i; j++;} 139 } 140 141 /* Each message would have a header = 1 + 2*(no of IS) + data */ 142 for (i=0; i<nrqs; i++) { 143 j = pa[i]; 144 w1[j] += w2[j] + 2*w3[j]; 145 msz += w1[j]; 146 } 147 148 /* Determine the number of messages to expect, their lengths, from from-ids */ 149 ierr = PetscGatherNumberOfMessages(comm,w2,w1,&nrqr);CHKERRQ(ierr); 150 ierr = PetscGatherMessageLengths(comm,nrqs,nrqr,w1,&onodes1,&olengths1);CHKERRQ(ierr); 151 152 /* Now post the Irecvs corresponding to these messages */ 153 ierr = PetscPostIrecvInt(comm,tag1,nrqr,onodes1,olengths1,&rbuf,&r_waits1);CHKERRQ(ierr); 154 155 /* Allocate Memory for outgoing messages */ 156 len = 2*size*sizeof(PetscInt*) + (size+msz)*sizeof(PetscInt); 157 ierr = PetscMalloc(len,&outdat);CHKERRQ(ierr); 158 ptr = outdat + size; /* Pointers to the data in outgoing buffers */ 159 ierr = PetscMemzero(outdat,2*size*sizeof(PetscInt*));CHKERRQ(ierr); 160 tmp = (PetscInt*)(outdat + 2*size); 161 ctr = tmp + msz; 162 163 { 164 PetscInt *iptr = tmp,ict = 0; 165 for (i=0; i<nrqs; i++) { 166 j = pa[i]; 167 iptr += ict; 168 outdat[j] = iptr; 169 ict = w1[j]; 170 } 171 } 172 173 /* Form the outgoing messages */ 174 /*plug in the headers*/ 175 for (i=0; i<nrqs; i++) { 176 j = pa[i]; 177 outdat[j][0] = 0; 178 ierr = PetscMemzero(outdat[j]+1,2*w3[j]*sizeof(PetscInt));CHKERRQ(ierr); 179 ptr[j] = outdat[j] + 2*w3[j] + 1; 180 } 181 182 /* Memory for doing local proc's work*/ 183 { 184 PetscInt *d_p; 185 char *t_p; 186 187 len = (imax)*(sizeof(PetscBT) + sizeof(PetscInt*)+ sizeof(PetscInt)) + 188 (Mbs)*imax*sizeof(PetscInt) + (Mbs/PETSC_BITS_PER_BYTE+1)*imax*sizeof(char) + 1; 189 ierr = PetscMalloc(len,&table);CHKERRQ(ierr); 190 ierr = PetscMemzero(table,len);CHKERRQ(ierr); 191 data = (PetscInt **)(table + imax); 192 isz = (PetscInt *)(data + imax); 193 d_p = (PetscInt *)(isz + imax); 194 t_p = (char *)(d_p + Mbs*imax); 195 for (i=0; i<imax; i++) { 196 table[i] = t_p + (Mbs/PETSC_BITS_PER_BYTE+1)*i; 197 data[i] = d_p + (Mbs)*i; 198 } 199 } 200 201 /* Parse the IS and update local tables and the outgoing buf with the data*/ 202 { 203 PetscInt n_i,*data_i,isz_i,*outdat_j,ctr_j; 204 PetscBT table_i; 205 206 for (i=0; i<imax; i++) { 207 ierr = PetscMemzero(ctr,size*sizeof(PetscInt));CHKERRQ(ierr); 208 n_i = n[i]; 209 table_i = table[i]; 210 idx_i = idx[i]; 211 data_i = data[i]; 212 isz_i = isz[i]; 213 for (j=0; j<n_i; j++) { /* parse the indices of each IS */ 214 row = idx_i[j]; 215 proc = rtable[row]; 216 if (proc != rank) { /* copy to the outgoing buffer */ 217 ctr[proc]++; 218 *ptr[proc] = row; 219 ptr[proc]++; 220 } 221 else { /* Update the local table */ 222 if (!PetscBTLookupSet(table_i,row)) { data_i[isz_i++] = row;} 223 } 224 } 225 /* Update the headers for the current IS */ 226 for (j=0; j<size; j++) { /* Can Optimise this loop by using pa[] */ 227 if ((ctr_j = ctr[j])) { 228 outdat_j = outdat[j]; 229 k = ++outdat_j[0]; 230 outdat_j[2*k] = ctr_j; 231 outdat_j[2*k-1] = i; 232 } 233 } 234 isz[i] = isz_i; 235 } 236 } 237 238 /* Now post the sends */ 239 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&s_waits1);CHKERRQ(ierr); 240 for (i=0; i<nrqs; ++i) { 241 j = pa[i]; 242 ierr = MPI_Isend(outdat[j],w1[j],MPIU_INT,j,tag1,comm,s_waits1+i);CHKERRQ(ierr); 243 } 244 245 /* No longer need the original indices*/ 246 for (i=0; i<imax; ++i) { 247 ierr = ISRestoreIndices(is[i],idx+i);CHKERRQ(ierr); 248 } 249 ierr = PetscFree(idx);CHKERRQ(ierr); 250 251 for (i=0; i<imax; ++i) { 252 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 253 } 254 255 /* Do Local work*/ 256 ierr = MatIncreaseOverlap_MPIBAIJ_Local(C,imax,table,isz,data);CHKERRQ(ierr); 257 258 /* Receive messages*/ 259 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&recv_status);CHKERRQ(ierr); 260 ierr = MPI_Waitall(nrqr,r_waits1,recv_status);CHKERRQ(ierr); 261 262 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 263 ierr = MPI_Waitall(nrqs,s_waits1,s_status);CHKERRQ(ierr); 264 265 /* Phase 1 sends are complete - deallocate buffers */ 266 ierr = PetscFree(outdat);CHKERRQ(ierr); 267 ierr = PetscFree(w1);CHKERRQ(ierr); 268 269 ierr = PetscMalloc((nrqr+1)*sizeof(PetscInt*),&xdata);CHKERRQ(ierr); 270 ierr = PetscMalloc((nrqr+1)*sizeof(PetscInt),&isz1);CHKERRQ(ierr); 271 ierr = MatIncreaseOverlap_MPIBAIJ_Receive(C,nrqr,rbuf,xdata,isz1);CHKERRQ(ierr); 272 ierr = PetscFree(rbuf);CHKERRQ(ierr); 273 274 /* Send the data back*/ 275 /* Do a global reduction to know the buffer space req for incoming messages*/ 276 { 277 PetscMPIInt *rw1; 278 279 ierr = PetscMalloc(size*sizeof(PetscInt),&rw1);CHKERRQ(ierr); 280 ierr = PetscMemzero(rw1,size*sizeof(PetscInt));CHKERRQ(ierr); 281 282 for (i=0; i<nrqr; ++i) { 283 proc = recv_status[i].MPI_SOURCE; 284 if (proc != onodes1[i]) SETERRQ(PETSC_ERR_PLIB,"MPI_SOURCE mismatch"); 285 rw1[proc] = isz1[i]; 286 } 287 288 ierr = PetscFree(onodes1);CHKERRQ(ierr); 289 ierr = PetscFree(olengths1);CHKERRQ(ierr); 290 291 /* Determine the number of messages to expect, their lengths, from from-ids */ 292 ierr = PetscGatherMessageLengths(comm,nrqr,nrqs,rw1,&onodes2,&olengths2);CHKERRQ(ierr); 293 PetscFree(rw1); 294 } 295 /* Now post the Irecvs corresponding to these messages */ 296 ierr = PetscPostIrecvInt(comm,tag2,nrqs,onodes2,olengths2,&rbuf2,&r_waits2);CHKERRQ(ierr); 297 298 /* Now post the sends */ 299 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits2);CHKERRQ(ierr); 300 for (i=0; i<nrqr; ++i) { 301 j = recv_status[i].MPI_SOURCE; 302 ierr = MPI_Isend(xdata[i],isz1[i],MPIU_INT,j,tag2,comm,s_waits2+i);CHKERRQ(ierr); 303 } 304 305 /* receive work done on other processors*/ 306 { 307 PetscMPIInt idex; 308 PetscInt is_no,ct1,max,*rbuf2_i,isz_i,*data_i,jmax; 309 PetscBT table_i; 310 MPI_Status *status2; 311 312 ierr = PetscMalloc((PetscMax(nrqr,nrqs)+1)*sizeof(MPI_Status),&status2);CHKERRQ(ierr); 313 314 for (i=0; i<nrqs; ++i) { 315 ierr = MPI_Waitany(nrqs,r_waits2,&idex,status2+i);CHKERRQ(ierr); 316 /* Process the message*/ 317 rbuf2_i = rbuf2[idex]; 318 ct1 = 2*rbuf2_i[0]+1; 319 jmax = rbuf2[idex][0]; 320 for (j=1; j<=jmax; j++) { 321 max = rbuf2_i[2*j]; 322 is_no = rbuf2_i[2*j-1]; 323 isz_i = isz[is_no]; 324 data_i = data[is_no]; 325 table_i = table[is_no]; 326 for (k=0; k<max; k++,ct1++) { 327 row = rbuf2_i[ct1]; 328 if (!PetscBTLookupSet(table_i,row)) { data_i[isz_i++] = row;} 329 } 330 isz[is_no] = isz_i; 331 } 332 } 333 ierr = MPI_Waitall(nrqr,s_waits2,status2);CHKERRQ(ierr); 334 ierr = PetscFree(status2);CHKERRQ(ierr); 335 } 336 337 for (i=0; i<imax; ++i) { 338 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz[i],data[i],is+i);CHKERRQ(ierr); 339 } 340 341 342 ierr = PetscFree(onodes2);CHKERRQ(ierr); 343 ierr = PetscFree(olengths2);CHKERRQ(ierr); 344 345 ierr = PetscFree(pa);CHKERRQ(ierr); 346 ierr = PetscFree(rbuf2);CHKERRQ(ierr); 347 ierr = PetscFree(s_waits1);CHKERRQ(ierr); 348 ierr = PetscFree(r_waits1);CHKERRQ(ierr); 349 ierr = PetscFree(s_waits2);CHKERRQ(ierr); 350 ierr = PetscFree(r_waits2);CHKERRQ(ierr); 351 ierr = PetscFree(table);CHKERRQ(ierr); 352 ierr = PetscFree(s_status);CHKERRQ(ierr); 353 ierr = PetscFree(recv_status);CHKERRQ(ierr); 354 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 355 ierr = PetscFree(xdata);CHKERRQ(ierr); 356 ierr = PetscFree(isz1);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "MatIncreaseOverlap_MPIBAIJ_Local" 362 /* 363 MatIncreaseOverlap_MPIBAIJ_Local - Called by MatincreaseOverlap, to do 364 the work on the local processor. 365 366 Inputs: 367 C - MAT_MPIBAIJ; 368 imax - total no of index sets processed at a time; 369 table - an array of char - size = Mbs bits. 370 371 Output: 372 isz - array containing the count of the solution elements corresponding 373 to each index set; 374 data - pointer to the solutions 375 */ 376 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Local(Mat C,PetscInt imax,PetscBT *table,PetscInt *isz,PetscInt **data) 377 { 378 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 379 Mat A = c->A,B = c->B; 380 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)B->data; 381 PetscInt start,end,val,max,rstart,cstart,*ai,*aj; 382 PetscInt *bi,*bj,*garray,i,j,k,row,*data_i,isz_i; 383 PetscBT table_i; 384 385 PetscFunctionBegin; 386 rstart = c->rstart; 387 cstart = c->cstart; 388 ai = a->i; 389 aj = a->j; 390 bi = b->i; 391 bj = b->j; 392 garray = c->garray; 393 394 395 for (i=0; i<imax; i++) { 396 data_i = data[i]; 397 table_i = table[i]; 398 isz_i = isz[i]; 399 for (j=0,max=isz[i]; j<max; j++) { 400 row = data_i[j] - rstart; 401 start = ai[row]; 402 end = ai[row+1]; 403 for (k=start; k<end; k++) { /* Amat */ 404 val = aj[k] + cstart; 405 if (!PetscBTLookupSet(table_i,val)) { data_i[isz_i++] = val;} 406 } 407 start = bi[row]; 408 end = bi[row+1]; 409 for (k=start; k<end; k++) { /* Bmat */ 410 val = garray[bj[k]]; 411 if (!PetscBTLookupSet(table_i,val)) { data_i[isz_i++] = val;} 412 } 413 } 414 isz[i] = isz_i; 415 } 416 PetscFunctionReturn(0); 417 } 418 #undef __FUNCT__ 419 #define __FUNCT__ "MatIncreaseOverlap_MPIBAIJ_Receive" 420 /* 421 MatIncreaseOverlap_MPIBAIJ_Receive - Process the recieved messages, 422 and return the output 423 424 Input: 425 C - the matrix 426 nrqr - no of messages being processed. 427 rbuf - an array of pointers to the recieved requests 428 429 Output: 430 xdata - array of messages to be sent back 431 isz1 - size of each message 432 433 For better efficiency perhaps we should malloc seperately each xdata[i], 434 then if a remalloc is required we need only copy the data for that one row 435 rather than all previous rows as it is now where a single large chunck of 436 memory is used. 437 438 */ 439 static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Receive(Mat C,PetscInt nrqr,PetscInt **rbuf,PetscInt **xdata,PetscInt * isz1) 440 { 441 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 442 Mat A = c->A,B = c->B; 443 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)B->data; 444 PetscErrorCode ierr; 445 PetscInt rstart,cstart,*ai,*aj,*bi,*bj,*garray,i,j,k; 446 PetscInt row,total_sz,ct,ct1,ct2,ct3,mem_estimate,oct2,l,start,end; 447 PetscInt val,max1,max2,rank,Mbs,no_malloc =0,*tmp,new_estimate,ctr; 448 PetscInt *rbuf_i,kmax,rbuf_0; 449 PetscBT xtable; 450 451 PetscFunctionBegin; 452 rank = c->rank; 453 Mbs = c->Mbs; 454 rstart = c->rstart; 455 cstart = c->cstart; 456 ai = a->i; 457 aj = a->j; 458 bi = b->i; 459 bj = b->j; 460 garray = c->garray; 461 462 463 for (i=0,ct=0,total_sz=0; i<nrqr; ++i) { 464 rbuf_i = rbuf[i]; 465 rbuf_0 = rbuf_i[0]; 466 ct += rbuf_0; 467 for (j=1; j<=rbuf_0; j++) { total_sz += rbuf_i[2*j]; } 468 } 469 470 if (c->Mbs) max1 = ct*(a->nz +b->nz)/c->Mbs; 471 else max1 = 1; 472 mem_estimate = 3*((total_sz > max1 ? total_sz : max1)+1); 473 ierr = PetscMalloc(mem_estimate*sizeof(PetscInt),&xdata[0]);CHKERRQ(ierr); 474 ++no_malloc; 475 ierr = PetscBTCreate(Mbs,xtable);CHKERRQ(ierr); 476 ierr = PetscMemzero(isz1,nrqr*sizeof(PetscInt));CHKERRQ(ierr); 477 478 ct3 = 0; 479 for (i=0; i<nrqr; i++) { /* for easch mesg from proc i */ 480 rbuf_i = rbuf[i]; 481 rbuf_0 = rbuf_i[0]; 482 ct1 = 2*rbuf_0+1; 483 ct2 = ct1; 484 ct3 += ct1; 485 for (j=1; j<=rbuf_0; j++) { /* for each IS from proc i*/ 486 ierr = PetscBTMemzero(Mbs,xtable);CHKERRQ(ierr); 487 oct2 = ct2; 488 kmax = rbuf_i[2*j]; 489 for (k=0; k<kmax; k++,ct1++) { 490 row = rbuf_i[ct1]; 491 if (!PetscBTLookupSet(xtable,row)) { 492 if (!(ct3 < mem_estimate)) { 493 new_estimate = (PetscInt)(1.5*mem_estimate)+1; 494 ierr = PetscMalloc(new_estimate * sizeof(PetscInt),&tmp);CHKERRQ(ierr); 495 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(PetscInt));CHKERRQ(ierr); 496 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 497 xdata[0] = tmp; 498 mem_estimate = new_estimate; ++no_malloc; 499 for (ctr=1; ctr<=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 500 } 501 xdata[i][ct2++] = row; 502 ct3++; 503 } 504 } 505 for (k=oct2,max2=ct2; k<max2; k++) { 506 row = xdata[i][k] - rstart; 507 start = ai[row]; 508 end = ai[row+1]; 509 for (l=start; l<end; l++) { 510 val = aj[l] + cstart; 511 if (!PetscBTLookupSet(xtable,val)) { 512 if (!(ct3 < mem_estimate)) { 513 new_estimate = (PetscInt)(1.5*mem_estimate)+1; 514 ierr = PetscMalloc(new_estimate * sizeof(PetscInt),&tmp);CHKERRQ(ierr); 515 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(PetscInt));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++] = val; 522 ct3++; 523 } 524 } 525 start = bi[row]; 526 end = bi[row+1]; 527 for (l=start; l<end; l++) { 528 val = garray[bj[l]]; 529 if (!PetscBTLookupSet(xtable,val)) { 530 if (!(ct3 < mem_estimate)) { 531 new_estimate = (PetscInt)(1.5*mem_estimate)+1; 532 ierr = PetscMalloc(new_estimate * sizeof(PetscInt),&tmp);CHKERRQ(ierr); 533 ierr = PetscMemcpy(tmp,xdata[0],mem_estimate*sizeof(PetscInt));CHKERRQ(ierr); 534 ierr = PetscFree(xdata[0]);CHKERRQ(ierr); 535 xdata[0] = tmp; 536 mem_estimate = new_estimate; ++no_malloc; 537 for (ctr =1; ctr <=i; ctr++) { xdata[ctr] = xdata[ctr-1] + isz1[ctr-1];} 538 } 539 xdata[i][ct2++] = val; 540 ct3++; 541 } 542 } 543 } 544 /* Update the header*/ 545 xdata[i][2*j] = ct2 - oct2; /* Undo the vector isz1 and use only a var*/ 546 xdata[i][2*j-1] = rbuf_i[2*j-1]; 547 } 548 xdata[i][0] = rbuf_0; 549 xdata[i+1] = xdata[i] + ct2; 550 isz1[i] = ct2; /* size of each message */ 551 } 552 ierr = PetscBTDestroy(xtable);CHKERRQ(ierr); 553 PetscLogInfo(0,"MatIncreaseOverlap_MPIBAIJ:[%d] Allocated %D bytes, required %D, no of mallocs = %D\n",rank,mem_estimate,ct3,no_malloc); 554 PetscFunctionReturn(0); 555 } 556 557 static PetscErrorCode MatGetSubMatrices_MPIBAIJ_local(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *); 558 559 #undef __FUNCT__ 560 #define __FUNCT__ "MatGetSubMatrices_MPIBAIJ" 561 PetscErrorCode MatGetSubMatrices_MPIBAIJ(Mat C,PetscInt ismax,const IS isrow[],const IS iscol[],MatReuse scall,Mat *submat[]) 562 { 563 IS *isrow_new,*iscol_new; 564 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 565 PetscErrorCode ierr; 566 PetscInt nmax,nstages_local,nstages,i,pos,max_no,N=C->N,bs=c->bs; 567 568 PetscFunctionBegin; 569 /* The compression and expansion should be avoided. Does'nt point 570 out errors might change the indices hence buggey */ 571 572 ierr = PetscMalloc(2*(ismax+1)*sizeof(IS),&isrow_new);CHKERRQ(ierr); 573 iscol_new = isrow_new + ismax; 574 ierr = ISCompressIndicesSorted(N,bs,ismax,isrow,isrow_new);CHKERRQ(ierr); 575 ierr = ISCompressIndicesSorted(N,bs,ismax,iscol,iscol_new);CHKERRQ(ierr); 576 577 /* Allocate memory to hold all the submatrices */ 578 if (scall != MAT_REUSE_MATRIX) { 579 ierr = PetscMalloc((ismax+1)*sizeof(Mat),submat);CHKERRQ(ierr); 580 } 581 /* Determine the number of stages through which submatrices are done */ 582 nmax = 20*1000000 / (c->Nbs * sizeof(PetscInt)); 583 if (!nmax) nmax = 1; 584 nstages_local = ismax/nmax + ((ismax % nmax)?1:0); 585 586 /* Make sure every processor loops through the nstages */ 587 ierr = MPI_Allreduce(&nstages_local,&nstages,1,MPIU_INT,MPI_MAX,C->comm);CHKERRQ(ierr); 588 for (i=0,pos=0; i<nstages; i++) { 589 if (pos+nmax <= ismax) max_no = nmax; 590 else if (pos == ismax) max_no = 0; 591 else max_no = ismax-pos; 592 ierr = MatGetSubMatrices_MPIBAIJ_local(C,max_no,isrow_new+pos,iscol_new+pos,scall,*submat+pos);CHKERRQ(ierr); 593 pos += max_no; 594 } 595 596 for (i=0; i<ismax; i++) { 597 ierr = ISDestroy(isrow_new[i]);CHKERRQ(ierr); 598 ierr = ISDestroy(iscol_new[i]);CHKERRQ(ierr); 599 } 600 ierr = PetscFree(isrow_new);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603 604 #if defined (PETSC_USE_CTABLE) 605 #undef __FUNCT__ 606 #define __FUNCT__ "PetscGetProc" 607 PetscErrorCode PetscGetProc(const PetscInt row, const PetscMPIInt size, const PetscInt proc_gnode[], PetscMPIInt *rank) 608 { 609 PetscInt nGlobalNd = proc_gnode[size]; 610 PetscMPIInt fproc = (PetscMPIInt) ((float)row * (float)size / (float)nGlobalNd + 0.5); 611 612 PetscFunctionBegin; 613 if (fproc > size) fproc = size; 614 while (row < proc_gnode[fproc] || row >= proc_gnode[fproc+1]) { 615 if (row < proc_gnode[fproc]) fproc--; 616 else fproc++; 617 } 618 *rank = fproc; 619 PetscFunctionReturn(0); 620 } 621 #endif 622 623 /* -------------------------------------------------------------------------*/ 624 /* This code is used for BAIJ and SBAIJ matrices (unfortunate dependency) */ 625 #undef __FUNCT__ 626 #define __FUNCT__ "MatGetSubMatrices_MPIBAIJ_local" 627 static PetscErrorCode MatGetSubMatrices_MPIBAIJ_local(Mat C,PetscInt ismax,const IS isrow[],const IS iscol[],MatReuse scall,Mat *submats) 628 { 629 Mat_MPIBAIJ *c = (Mat_MPIBAIJ*)C->data; 630 Mat A = c->A; 631 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)c->B->data,*mat; 632 PetscInt **irow,**icol,*nrow,*ncol,*w3,*w4,start; 633 PetscErrorCode ierr; 634 PetscMPIInt size,tag0,tag1,tag2,tag3,*w1,*w2,nrqr,idex,end; 635 PetscInt **sbuf1,**sbuf2,rank,i,j,k,l,ct1,ct2,**rbuf1,row,proc; 636 PetscInt nrqs,msz,**ptr,*req_size,*ctr,*pa,*tmp,tcol; 637 PetscInt **rbuf3,*req_source,**sbuf_aj,**rbuf2,max1,max2; 638 PetscInt **lens,is_no,ncols,*cols,mat_i,*mat_j,tmp2,jmax,*irow_i; 639 PetscInt len,ctr_j,*sbuf1_j,*sbuf_aj_i,*rbuf1_i,kmax,*lens_i; 640 PetscInt bs=c->bs,bs2=c->bs2,*a_j=a->j,*b_j=b->j,*cworkA,*cworkB; 641 PetscInt cstart = c->cstart,nzA,nzB,*a_i=a->i,*b_i=b->i,imark; 642 PetscInt *bmap = c->garray,ctmp,rstart=c->rstart; 643 MPI_Request *s_waits1,*r_waits1,*s_waits2,*r_waits2,*r_waits3; 644 MPI_Request *r_waits4,*s_waits3,*s_waits4; 645 MPI_Status *r_status1,*r_status2,*s_status1,*s_status3,*s_status2; 646 MPI_Status *r_status3,*r_status4,*s_status4; 647 MPI_Comm comm; 648 MatScalar **rbuf4,**sbuf_aa,*vals,*mat_a,*sbuf_aa_i,*vworkA,*vworkB; 649 MatScalar *a_a=a->a,*b_a=b->a; 650 PetscTruth flag; 651 PetscMPIInt *onodes1,*olengths1; 652 653 #if defined (PETSC_USE_CTABLE) 654 PetscInt tt; 655 PetscTable *rowmaps,*colmaps,lrow1_grow1,lcol1_gcol1; 656 #else 657 PetscInt **cmap,*cmap_i,*rtable,*rmap_i,**rmap, Mbs = c->Mbs; 658 #endif 659 660 PetscFunctionBegin; 661 comm = C->comm; 662 tag0 = C->tag; 663 size = c->size; 664 rank = c->rank; 665 666 /* Get some new tags to keep the communication clean */ 667 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 668 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 669 ierr = PetscObjectGetNewTag((PetscObject)C,&tag3);CHKERRQ(ierr); 670 671 /* Check if the col indices are sorted */ 672 for (i=0; i<ismax; i++) { 673 ierr = ISSorted(iscol[i],(PetscTruth*)&j);CHKERRQ(ierr); 674 if (!j) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IS is not sorted"); 675 } 676 677 len = (2*ismax+1)*(sizeof(PetscInt*)+ sizeof(PetscInt)); 678 #if !defined (PETSC_USE_CTABLE) 679 len += (Mbs+1)*sizeof(PetscInt); 680 #endif 681 ierr = PetscMalloc(len,&irow);CHKERRQ(ierr); 682 icol = irow + ismax; 683 nrow = (PetscInt*)(icol + ismax); 684 ncol = nrow + ismax; 685 #if !defined (PETSC_USE_CTABLE) 686 rtable = ncol + ismax; 687 /* Create hash table for the mapping :row -> proc*/ 688 for (i=0,j=0; i<size; i++) { 689 jmax = c->rowners[i+1]; 690 for (; j<jmax; j++) { 691 rtable[j] = i; 692 } 693 } 694 #endif 695 696 for (i=0; i<ismax; i++) { 697 ierr = ISGetIndices(isrow[i],&irow[i]);CHKERRQ(ierr); 698 ierr = ISGetIndices(iscol[i],&icol[i]);CHKERRQ(ierr); 699 ierr = ISGetLocalSize(isrow[i],&nrow[i]);CHKERRQ(ierr); 700 ierr = ISGetLocalSize(iscol[i],&ncol[i]);CHKERRQ(ierr); 701 } 702 703 /* evaluate communication - mesg to who,length of mesg,and buffer space 704 required. Based on this, buffers are allocated, and data copied into them*/ 705 ierr = PetscMalloc(size*2*sizeof(PetscInt) + size*2*sizeof(PetscMPIInt),&w1);CHKERRQ(ierr); /* mesg size */ 706 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 707 w3 = (PetscInt*)(w2 + size); /* no of IS that needs to be sent to proc i */ 708 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 709 ierr = PetscMemzero(w1,size*sizeof(PetscInt)+2*size*sizeof(PetscMPIInt));CHKERRQ(ierr); /* initialise work vector*/ 710 for (i=0; i<ismax; i++) { 711 ierr = PetscMemzero(w4,size*sizeof(PetscInt));CHKERRQ(ierr); /* initialise work vector*/ 712 jmax = nrow[i]; 713 irow_i = irow[i]; 714 for (j=0; j<jmax; j++) { 715 row = irow_i[j]; 716 #if defined (PETSC_USE_CTABLE) 717 ierr = PetscGetProc(row,size,c->rowners,&proc);CHKERRQ(ierr); 718 #else 719 proc = rtable[row]; 720 #endif 721 w4[proc]++; 722 } 723 for (j=0; j<size; j++) { 724 if (w4[j]) { w1[j] += w4[j]; w3[j]++;} 725 } 726 } 727 728 nrqs = 0; /* no of outgoing messages */ 729 msz = 0; /* total mesg length for all proc */ 730 w1[rank] = 0; /* no mesg sent to intself */ 731 w3[rank] = 0; 732 for (i=0; i<size; i++) { 733 if (w1[i]) { w2[i] = 1; nrqs++;} /* there exists a message to proc i */ 734 } 735 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt),&pa);CHKERRQ(ierr); /*(proc -array)*/ 736 for (i=0,j=0; i<size; i++) { 737 if (w1[i]) { pa[j] = i; j++; } 738 } 739 740 /* Each message would have a header = 1 + 2*(no of IS) + data */ 741 for (i=0; i<nrqs; i++) { 742 j = pa[i]; 743 w1[j] += w2[j] + 2* w3[j]; 744 msz += w1[j]; 745 } 746 747 /* Determine the number of messages to expect, their lengths, from from-ids */ 748 ierr = PetscGatherNumberOfMessages(comm,w2,w1,&nrqr);CHKERRQ(ierr); 749 ierr = PetscGatherMessageLengths(comm,nrqs,nrqr,w1,&onodes1,&olengths1);CHKERRQ(ierr); 750 751 /* Now post the Irecvs corresponding to these messages */ 752 ierr = PetscPostIrecvInt(comm,tag0,nrqr,onodes1,olengths1,&rbuf1,&r_waits1);CHKERRQ(ierr); 753 754 ierr = PetscFree(onodes1);CHKERRQ(ierr); 755 ierr = PetscFree(olengths1);CHKERRQ(ierr); 756 757 /* Allocate Memory for outgoing messages */ 758 len = 2*size*sizeof(PetscInt*) + 2*msz*sizeof(PetscInt) + size*sizeof(PetscInt); 759 ierr = PetscMalloc(len,&sbuf1);CHKERRQ(ierr); 760 ptr = sbuf1 + size; /* Pointers to the data in outgoing buffers */ 761 ierr = PetscMemzero(sbuf1,2*size*sizeof(PetscInt*));CHKERRQ(ierr); 762 /* allocate memory for outgoing data + buf to receive the first reply */ 763 tmp = (PetscInt*)(ptr + size); 764 ctr = tmp + 2*msz; 765 766 { 767 PetscInt *iptr = tmp,ict = 0; 768 for (i=0; i<nrqs; i++) { 769 j = pa[i]; 770 iptr += ict; 771 sbuf1[j] = iptr; 772 ict = w1[j]; 773 } 774 } 775 776 /* Form the outgoing messages */ 777 /* Initialise the header space */ 778 for (i=0; i<nrqs; i++) { 779 j = pa[i]; 780 sbuf1[j][0] = 0; 781 ierr = PetscMemzero(sbuf1[j]+1,2*w3[j]*sizeof(PetscInt));CHKERRQ(ierr); 782 ptr[j] = sbuf1[j] + 2*w3[j] + 1; 783 } 784 785 /* Parse the isrow and copy data into outbuf */ 786 for (i=0; i<ismax; i++) { 787 ierr = PetscMemzero(ctr,size*sizeof(PetscInt));CHKERRQ(ierr); 788 irow_i = irow[i]; 789 jmax = nrow[i]; 790 for (j=0; j<jmax; j++) { /* parse the indices of each IS */ 791 row = irow_i[j]; 792 #if defined (PETSC_USE_CTABLE) 793 ierr = PetscGetProc(row,size,c->rowners,&proc);CHKERRQ(ierr); 794 #else 795 proc = rtable[row]; 796 #endif 797 if (proc != rank) { /* copy to the outgoing buf*/ 798 ctr[proc]++; 799 *ptr[proc] = row; 800 ptr[proc]++; 801 } 802 } 803 /* Update the headers for the current IS */ 804 for (j=0; j<size; j++) { /* Can Optimise this loop too */ 805 if ((ctr_j = ctr[j])) { 806 sbuf1_j = sbuf1[j]; 807 k = ++sbuf1_j[0]; 808 sbuf1_j[2*k] = ctr_j; 809 sbuf1_j[2*k-1] = i; 810 } 811 } 812 } 813 814 /* Now post the sends */ 815 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&s_waits1);CHKERRQ(ierr); 816 for (i=0; i<nrqs; ++i) { 817 j = pa[i]; 818 ierr = MPI_Isend(sbuf1[j],w1[j],MPIU_INT,j,tag0,comm,s_waits1+i);CHKERRQ(ierr); 819 } 820 821 /* Post Recieves to capture the buffer size */ 822 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits2);CHKERRQ(ierr); 823 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt*),&rbuf2);CHKERRQ(ierr); 824 rbuf2[0] = tmp + msz; 825 for (i=1; i<nrqs; ++i) { 826 j = pa[i]; 827 rbuf2[i] = rbuf2[i-1]+w1[pa[i-1]]; 828 } 829 for (i=0; i<nrqs; ++i) { 830 j = pa[i]; 831 ierr = MPI_Irecv(rbuf2[i],w1[j],MPIU_INT,j,tag1,comm,r_waits2+i);CHKERRQ(ierr); 832 } 833 834 /* Send to other procs the buf size they should allocate */ 835 836 /* Receive messages*/ 837 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits2);CHKERRQ(ierr); 838 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&r_status1);CHKERRQ(ierr); 839 len = 2*nrqr*sizeof(PetscInt) + (nrqr+1)*sizeof(PetscInt*); 840 ierr = PetscMalloc(len,&sbuf2);CHKERRQ(ierr); 841 req_size = (PetscInt*)(sbuf2 + nrqr); 842 req_source = req_size + nrqr; 843 844 { 845 Mat_SeqBAIJ *sA = (Mat_SeqBAIJ*)c->A->data,*sB = (Mat_SeqBAIJ*)c->B->data; 846 PetscInt *sAi = sA->i,*sBi = sB->i,id,*sbuf2_i; 847 848 for (i=0; i<nrqr; ++i) { 849 ierr = MPI_Waitany(nrqr,r_waits1,&idex,r_status1+i);CHKERRQ(ierr); 850 req_size[idex] = 0; 851 rbuf1_i = rbuf1[idex]; 852 start = 2*rbuf1_i[0] + 1; 853 ierr = MPI_Get_count(r_status1+i,MPIU_INT,&end);CHKERRQ(ierr); 854 ierr = PetscMalloc(end*sizeof(PetscInt),&sbuf2[idex]);CHKERRQ(ierr); 855 sbuf2_i = sbuf2[idex]; 856 for (j=start; j<end; j++) { 857 id = rbuf1_i[j] - rstart; 858 ncols = sAi[id+1] - sAi[id] + sBi[id+1] - sBi[id]; 859 sbuf2_i[j] = ncols; 860 req_size[idex] += ncols; 861 } 862 req_source[idex] = r_status1[i].MPI_SOURCE; 863 /* form the header */ 864 sbuf2_i[0] = req_size[idex]; 865 for (j=1; j<start; j++) { sbuf2_i[j] = rbuf1_i[j]; } 866 ierr = MPI_Isend(sbuf2_i,end,MPIU_INT,req_source[idex],tag1,comm,s_waits2+i);CHKERRQ(ierr); 867 } 868 } 869 ierr = PetscFree(r_status1);CHKERRQ(ierr); 870 ierr = PetscFree(r_waits1);CHKERRQ(ierr); 871 872 /* recv buffer sizes */ 873 /* Receive messages*/ 874 875 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt*),&rbuf3);CHKERRQ(ierr); 876 ierr = PetscMalloc((nrqs+1)*sizeof(MatScalar*),&rbuf4);CHKERRQ(ierr); 877 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits3);CHKERRQ(ierr); 878 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits4);CHKERRQ(ierr); 879 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status2);CHKERRQ(ierr); 880 881 for (i=0; i<nrqs; ++i) { 882 ierr = MPI_Waitany(nrqs,r_waits2,&idex,r_status2+i);CHKERRQ(ierr); 883 ierr = PetscMalloc(rbuf2[idex][0]*sizeof(PetscInt),&rbuf3[idex]);CHKERRQ(ierr); 884 ierr = PetscMalloc(rbuf2[idex][0]*bs2*sizeof(MatScalar),&rbuf4[idex]);CHKERRQ(ierr); 885 ierr = MPI_Irecv(rbuf3[idex],rbuf2[idex][0],MPIU_INT,r_status2[i].MPI_SOURCE,tag2,comm,r_waits3+idex);CHKERRQ(ierr); 886 ierr = MPI_Irecv(rbuf4[idex],rbuf2[idex][0]*bs2,MPIU_MATSCALAR,r_status2[i].MPI_SOURCE,tag3,comm,r_waits4+idex);CHKERRQ(ierr); 887 } 888 ierr = PetscFree(r_status2);CHKERRQ(ierr); 889 ierr = PetscFree(r_waits2);CHKERRQ(ierr); 890 891 /* Wait on sends1 and sends2 */ 892 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&s_status1);CHKERRQ(ierr); 893 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status2);CHKERRQ(ierr); 894 895 ierr = MPI_Waitall(nrqs,s_waits1,s_status1);CHKERRQ(ierr); 896 ierr = MPI_Waitall(nrqr,s_waits2,s_status2);CHKERRQ(ierr); 897 ierr = PetscFree(s_status1);CHKERRQ(ierr); 898 ierr = PetscFree(s_status2);CHKERRQ(ierr); 899 ierr = PetscFree(s_waits1);CHKERRQ(ierr); 900 ierr = PetscFree(s_waits2);CHKERRQ(ierr); 901 902 /* Now allocate buffers for a->j, and send them off */ 903 ierr = PetscMalloc((nrqr+1)*sizeof(PetscInt*),&sbuf_aj);CHKERRQ(ierr); 904 for (i=0,j=0; i<nrqr; i++) j += req_size[i]; 905 ierr = PetscMalloc((j+1)*sizeof(PetscInt),&sbuf_aj[0]);CHKERRQ(ierr); 906 for (i=1; i<nrqr; i++) sbuf_aj[i] = sbuf_aj[i-1] + req_size[i-1]; 907 908 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits3);CHKERRQ(ierr); 909 { 910 for (i=0; i<nrqr; i++) { 911 rbuf1_i = rbuf1[i]; 912 sbuf_aj_i = sbuf_aj[i]; 913 ct1 = 2*rbuf1_i[0] + 1; 914 ct2 = 0; 915 for (j=1,max1=rbuf1_i[0]; j<=max1; j++) { 916 kmax = rbuf1[i][2*j]; 917 for (k=0; k<kmax; k++,ct1++) { 918 row = rbuf1_i[ct1] - rstart; 919 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 920 ncols = nzA + nzB; 921 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 922 923 /* load the column indices for this row into cols*/ 924 cols = sbuf_aj_i + ct2; 925 for (l=0; l<nzB; l++) { 926 if ((ctmp = bmap[cworkB[l]]) < cstart) cols[l] = ctmp; 927 else break; 928 } 929 imark = l; 930 for (l=0; l<nzA; l++) cols[imark+l] = cstart + cworkA[l]; 931 for (l=imark; l<nzB; l++) cols[nzA+l] = bmap[cworkB[l]]; 932 ct2 += ncols; 933 } 934 } 935 ierr = MPI_Isend(sbuf_aj_i,req_size[i],MPIU_INT,req_source[i],tag2,comm,s_waits3+i);CHKERRQ(ierr); 936 } 937 } 938 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status3);CHKERRQ(ierr); 939 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status3);CHKERRQ(ierr); 940 941 /* Allocate buffers for a->a, and send them off */ 942 ierr = PetscMalloc((nrqr+1)*sizeof(MatScalar *),&sbuf_aa);CHKERRQ(ierr); 943 for (i=0,j=0; i<nrqr; i++) j += req_size[i]; 944 ierr = PetscMalloc((j+1)*bs2*sizeof(MatScalar),&sbuf_aa[0]);CHKERRQ(ierr); 945 for (i=1; i<nrqr; i++) sbuf_aa[i] = sbuf_aa[i-1] + req_size[i-1]*bs2; 946 947 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits4);CHKERRQ(ierr); 948 { 949 for (i=0; i<nrqr; i++) { 950 rbuf1_i = rbuf1[i]; 951 sbuf_aa_i = sbuf_aa[i]; 952 ct1 = 2*rbuf1_i[0]+1; 953 ct2 = 0; 954 for (j=1,max1=rbuf1_i[0]; j<=max1; j++) { 955 kmax = rbuf1_i[2*j]; 956 for (k=0; k<kmax; k++,ct1++) { 957 row = rbuf1_i[ct1] - rstart; 958 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 959 ncols = nzA + nzB; 960 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 961 vworkA = a_a + a_i[row]*bs2; vworkB = b_a + b_i[row]*bs2; 962 963 /* load the column values for this row into vals*/ 964 vals = sbuf_aa_i+ct2*bs2; 965 for (l=0; l<nzB; l++) { 966 if ((bmap[cworkB[l]]) < cstart) { 967 ierr = PetscMemcpy(vals+l*bs2,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 968 } 969 else break; 970 } 971 imark = l; 972 for (l=0; l<nzA; l++) { 973 ierr = PetscMemcpy(vals+(imark+l)*bs2,vworkA+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 974 } 975 for (l=imark; l<nzB; l++) { 976 ierr = PetscMemcpy(vals+(nzA+l)*bs2,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 977 } 978 ct2 += ncols; 979 } 980 } 981 ierr = MPI_Isend(sbuf_aa_i,req_size[i]*bs2,MPIU_MATSCALAR,req_source[i],tag3,comm,s_waits4+i);CHKERRQ(ierr); 982 } 983 } 984 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status4);CHKERRQ(ierr); 985 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status4);CHKERRQ(ierr); 986 ierr = PetscFree(rbuf1);CHKERRQ(ierr); 987 988 /* Form the matrix */ 989 /* create col map */ 990 { 991 PetscInt *icol_i; 992 #if defined (PETSC_USE_CTABLE) 993 /* Create row map*/ 994 ierr = PetscMalloc((1+ismax)*sizeof(PetscTable),&colmaps);CHKERRQ(ierr); 995 for (i=0; i<ismax; i++) { 996 ierr = PetscTableCreate(ncol[i]+1,&colmaps[i]);CHKERRQ(ierr); 997 } 998 #else 999 len = (1+ismax)*sizeof(PetscInt*)+ ismax*c->Nbs*sizeof(PetscInt); 1000 ierr = PetscMalloc(len,&cmap);CHKERRQ(ierr); 1001 cmap[0] = (PetscInt*)(cmap + ismax); 1002 ierr = PetscMemzero(cmap[0],(1+ismax*c->Nbs)*sizeof(PetscInt));CHKERRQ(ierr); 1003 for (i=1; i<ismax; i++) { cmap[i] = cmap[i-1] + c->Nbs; } 1004 #endif 1005 for (i=0; i<ismax; i++) { 1006 jmax = ncol[i]; 1007 icol_i = icol[i]; 1008 #if defined (PETSC_USE_CTABLE) 1009 lcol1_gcol1 = colmaps[i]; 1010 for (j=0; j<jmax; j++) { 1011 ierr = PetscTableAdd(lcol1_gcol1,icol_i[j]+1,j+1);CHKERRQ(ierr); 1012 } 1013 #else 1014 cmap_i = cmap[i]; 1015 for (j=0; j<jmax; j++) { 1016 cmap_i[icol_i[j]] = j+1; 1017 } 1018 #endif 1019 } 1020 } 1021 1022 /* Create lens which is required for MatCreate... */ 1023 for (i=0,j=0; i<ismax; i++) { j += nrow[i]; } 1024 len = (1+ismax)*sizeof(PetscInt*)+ j*sizeof(PetscInt); 1025 ierr = PetscMalloc(len,&lens);CHKERRQ(ierr); 1026 lens[0] = (PetscInt*)(lens + ismax); 1027 ierr = PetscMemzero(lens[0],j*sizeof(PetscInt));CHKERRQ(ierr); 1028 for (i=1; i<ismax; i++) { lens[i] = lens[i-1] + nrow[i-1]; } 1029 1030 /* Update lens from local data */ 1031 for (i=0; i<ismax; i++) { 1032 jmax = nrow[i]; 1033 #if defined (PETSC_USE_CTABLE) 1034 lcol1_gcol1 = colmaps[i]; 1035 #else 1036 cmap_i = cmap[i]; 1037 #endif 1038 irow_i = irow[i]; 1039 lens_i = lens[i]; 1040 for (j=0; j<jmax; j++) { 1041 row = irow_i[j]; 1042 #if defined (PETSC_USE_CTABLE) 1043 ierr = PetscGetProc(row,size,c->rowners,&proc);CHKERRQ(ierr); 1044 #else 1045 proc = rtable[row]; 1046 #endif 1047 if (proc == rank) { 1048 /* Get indices from matA and then from matB */ 1049 row = row - rstart; 1050 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 1051 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 1052 #if defined (PETSC_USE_CTABLE) 1053 for (k=0; k<nzA; k++) { 1054 ierr = PetscTableFind(lcol1_gcol1,cstart+cworkA[k]+1,&tt);CHKERRQ(ierr); 1055 if (tt) { lens_i[j]++; } 1056 } 1057 for (k=0; k<nzB; k++) { 1058 ierr = PetscTableFind(lcol1_gcol1,bmap[cworkB[k]]+1,&tt);CHKERRQ(ierr); 1059 if (tt) { lens_i[j]++; } 1060 } 1061 #else 1062 for (k=0; k<nzA; k++) { 1063 if (cmap_i[cstart + cworkA[k]]) { lens_i[j]++; } 1064 } 1065 for (k=0; k<nzB; k++) { 1066 if (cmap_i[bmap[cworkB[k]]]) { lens_i[j]++; } 1067 } 1068 #endif 1069 } 1070 } 1071 } 1072 #if defined (PETSC_USE_CTABLE) 1073 /* Create row map*/ 1074 ierr = PetscMalloc((1+ismax)*sizeof(PetscTable),&rowmaps);CHKERRQ(ierr); 1075 for (i=0; i<ismax; i++){ 1076 ierr = PetscTableCreate(nrow[i]+1,&rowmaps[i]);CHKERRQ(ierr); 1077 } 1078 #else 1079 /* Create row map*/ 1080 len = (1+ismax)*sizeof(PetscInt*)+ ismax*Mbs*sizeof(PetscInt); 1081 ierr = PetscMalloc(len,&rmap);CHKERRQ(ierr); 1082 rmap[0] = (PetscInt*)(rmap + ismax); 1083 ierr = PetscMemzero(rmap[0],ismax*Mbs*sizeof(PetscInt));CHKERRQ(ierr); 1084 for (i=1; i<ismax; i++) { rmap[i] = rmap[i-1] + Mbs;} 1085 #endif 1086 for (i=0; i<ismax; i++) { 1087 irow_i = irow[i]; 1088 jmax = nrow[i]; 1089 #if defined (PETSC_USE_CTABLE) 1090 lrow1_grow1 = rowmaps[i]; 1091 for (j=0; j<jmax; j++) { 1092 ierr = PetscTableAdd(lrow1_grow1,irow_i[j]+1,j+1);CHKERRQ(ierr); 1093 } 1094 #else 1095 rmap_i = rmap[i]; 1096 for (j=0; j<jmax; j++) { 1097 rmap_i[irow_i[j]] = j; 1098 } 1099 #endif 1100 } 1101 1102 /* Update lens from offproc data */ 1103 { 1104 PetscInt *rbuf2_i,*rbuf3_i,*sbuf1_i; 1105 PetscMPIInt ii; 1106 1107 for (tmp2=0; tmp2<nrqs; tmp2++) { 1108 ierr = MPI_Waitany(nrqs,r_waits3,&ii,r_status3+tmp2);CHKERRQ(ierr); 1109 idex = pa[ii]; 1110 sbuf1_i = sbuf1[idex]; 1111 jmax = sbuf1_i[0]; 1112 ct1 = 2*jmax+1; 1113 ct2 = 0; 1114 rbuf2_i = rbuf2[ii]; 1115 rbuf3_i = rbuf3[ii]; 1116 for (j=1; j<=jmax; j++) { 1117 is_no = sbuf1_i[2*j-1]; 1118 max1 = sbuf1_i[2*j]; 1119 lens_i = lens[is_no]; 1120 #if defined (PETSC_USE_CTABLE) 1121 lcol1_gcol1 = colmaps[is_no]; 1122 lrow1_grow1 = rowmaps[is_no]; 1123 #else 1124 cmap_i = cmap[is_no]; 1125 rmap_i = rmap[is_no]; 1126 #endif 1127 for (k=0; k<max1; k++,ct1++) { 1128 #if defined (PETSC_USE_CTABLE) 1129 ierr = PetscTableFind(lrow1_grow1,sbuf1_i[ct1]+1,&row);CHKERRQ(ierr); 1130 row--; 1131 if (row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1132 #else 1133 row = rmap_i[sbuf1_i[ct1]]; /* the val in the new matrix to be */ 1134 #endif 1135 max2 = rbuf2_i[ct1]; 1136 for (l=0; l<max2; l++,ct2++) { 1137 #if defined (PETSC_USE_CTABLE) 1138 ierr = PetscTableFind(lcol1_gcol1,rbuf3_i[ct2]+1,&tt);CHKERRQ(ierr); 1139 if (tt) { 1140 lens_i[row]++; 1141 } 1142 #else 1143 if (cmap_i[rbuf3_i[ct2]]) { 1144 lens_i[row]++; 1145 } 1146 #endif 1147 } 1148 } 1149 } 1150 } 1151 } 1152 ierr = PetscFree(r_status3);CHKERRQ(ierr); 1153 ierr = PetscFree(r_waits3);CHKERRQ(ierr); 1154 ierr = MPI_Waitall(nrqr,s_waits3,s_status3);CHKERRQ(ierr); 1155 ierr = PetscFree(s_status3);CHKERRQ(ierr); 1156 ierr = PetscFree(s_waits3);CHKERRQ(ierr); 1157 1158 /* Create the submatrices */ 1159 if (scall == MAT_REUSE_MATRIX) { 1160 /* 1161 Assumes new rows are same length as the old rows, hence bug! 1162 */ 1163 for (i=0; i<ismax; i++) { 1164 mat = (Mat_SeqBAIJ *)(submats[i]->data); 1165 if ((mat->mbs != nrow[i]) || (mat->nbs != ncol[i] || mat->bs != bs)) { 1166 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1167 } 1168 ierr = PetscMemcmp(mat->ilen,lens[i],mat->mbs *sizeof(PetscInt),&flag);CHKERRQ(ierr); 1169 if (flag == PETSC_FALSE) { 1170 SETERRQ(PETSC_ERR_ARG_INCOMP,"Cannot reuse matrix. wrong no of nonzeros"); 1171 } 1172 /* Initial matrix as if empty */ 1173 ierr = PetscMemzero(mat->ilen,mat->mbs*sizeof(PetscInt));CHKERRQ(ierr); 1174 submats[i]->factor = C->factor; 1175 } 1176 } else { 1177 for (i=0; i<ismax; i++) { 1178 ierr = MatCreate(PETSC_COMM_SELF,nrow[i]*bs,ncol[i]*bs,nrow[i]*bs,ncol[i]*bs,submats+i);CHKERRQ(ierr); 1179 ierr = MatSetType(submats[i],A->type_name);CHKERRQ(ierr); 1180 ierr = MatSeqBAIJSetPreallocation(submats[i],a->bs,0,lens[i]);CHKERRQ(ierr); 1181 #if !defined(PETSC_USE_64BIT_INT) 1182 ierr = MatSeqSBAIJSetPreallocation(submats[i],a->bs,0,lens[i]);CHKERRQ(ierr); 1183 #endif 1184 } 1185 } 1186 1187 /* Assemble the matrices */ 1188 /* First assemble the local rows */ 1189 { 1190 PetscInt ilen_row,*imat_ilen,*imat_j,*imat_i; 1191 MatScalar *imat_a; 1192 1193 for (i=0; i<ismax; i++) { 1194 mat = (Mat_SeqBAIJ*)submats[i]->data; 1195 imat_ilen = mat->ilen; 1196 imat_j = mat->j; 1197 imat_i = mat->i; 1198 imat_a = mat->a; 1199 1200 #if defined (PETSC_USE_CTABLE) 1201 lcol1_gcol1 = colmaps[i]; 1202 lrow1_grow1 = rowmaps[i]; 1203 #else 1204 cmap_i = cmap[i]; 1205 rmap_i = rmap[i]; 1206 #endif 1207 irow_i = irow[i]; 1208 jmax = nrow[i]; 1209 for (j=0; j<jmax; j++) { 1210 row = irow_i[j]; 1211 #if defined (PETSC_USE_CTABLE) 1212 ierr = PetscGetProc(row,size,c->rowners,&proc);CHKERRQ(ierr); 1213 #else 1214 proc = rtable[row]; 1215 #endif 1216 if (proc == rank) { 1217 row = row - rstart; 1218 nzA = a_i[row+1] - a_i[row]; 1219 nzB = b_i[row+1] - b_i[row]; 1220 cworkA = a_j + a_i[row]; 1221 cworkB = b_j + b_i[row]; 1222 vworkA = a_a + a_i[row]*bs2; 1223 vworkB = b_a + b_i[row]*bs2; 1224 #if defined (PETSC_USE_CTABLE) 1225 ierr = PetscTableFind(lrow1_grow1,row+rstart+1,&row);CHKERRQ(ierr); 1226 row--; 1227 if (row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1228 #else 1229 row = rmap_i[row + rstart]; 1230 #endif 1231 mat_i = imat_i[row]; 1232 mat_a = imat_a + mat_i*bs2; 1233 mat_j = imat_j + mat_i; 1234 ilen_row = imat_ilen[row]; 1235 1236 /* load the column indices for this row into cols*/ 1237 for (l=0; l<nzB; l++) { 1238 if ((ctmp = bmap[cworkB[l]]) < cstart) { 1239 #if defined (PETSC_USE_CTABLE) 1240 ierr = PetscTableFind(lcol1_gcol1,ctmp+1,&tcol);CHKERRQ(ierr); 1241 if (tcol) { 1242 #else 1243 if ((tcol = cmap_i[ctmp])) { 1244 #endif 1245 *mat_j++ = tcol - 1; 1246 ierr = PetscMemcpy(mat_a,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1247 mat_a += bs2; 1248 ilen_row++; 1249 } 1250 } else break; 1251 } 1252 imark = l; 1253 for (l=0; l<nzA; l++) { 1254 #if defined (PETSC_USE_CTABLE) 1255 ierr = PetscTableFind(lcol1_gcol1,cstart+cworkA[l]+1,&tcol);CHKERRQ(ierr); 1256 if (tcol) { 1257 #else 1258 if ((tcol = cmap_i[cstart + cworkA[l]])) { 1259 #endif 1260 *mat_j++ = tcol - 1; 1261 ierr = PetscMemcpy(mat_a,vworkA+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1262 mat_a += bs2; 1263 ilen_row++; 1264 } 1265 } 1266 for (l=imark; l<nzB; l++) { 1267 #if defined (PETSC_USE_CTABLE) 1268 ierr = PetscTableFind(lcol1_gcol1,bmap[cworkB[l]]+1,&tcol);CHKERRQ(ierr); 1269 if (tcol) { 1270 #else 1271 if ((tcol = cmap_i[bmap[cworkB[l]]])) { 1272 #endif 1273 *mat_j++ = tcol - 1; 1274 ierr = PetscMemcpy(mat_a,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1275 mat_a += bs2; 1276 ilen_row++; 1277 } 1278 } 1279 imat_ilen[row] = ilen_row; 1280 } 1281 } 1282 1283 } 1284 } 1285 1286 /* Now assemble the off proc rows*/ 1287 { 1288 PetscInt *sbuf1_i,*rbuf2_i,*rbuf3_i,*imat_ilen,ilen; 1289 PetscInt *imat_j,*imat_i; 1290 MatScalar *imat_a,*rbuf4_i; 1291 PetscMPIInt ii; 1292 1293 for (tmp2=0; tmp2<nrqs; tmp2++) { 1294 ierr = MPI_Waitany(nrqs,r_waits4,&ii,r_status4+tmp2);CHKERRQ(ierr); 1295 idex = pa[ii]; 1296 sbuf1_i = sbuf1[idex]; 1297 jmax = sbuf1_i[0]; 1298 ct1 = 2*jmax + 1; 1299 ct2 = 0; 1300 rbuf2_i = rbuf2[ii]; 1301 rbuf3_i = rbuf3[ii]; 1302 rbuf4_i = rbuf4[ii]; 1303 for (j=1; j<=jmax; j++) { 1304 is_no = sbuf1_i[2*j-1]; 1305 #if defined (PETSC_USE_CTABLE) 1306 lrow1_grow1 = rowmaps[is_no]; 1307 lcol1_gcol1 = colmaps[is_no]; 1308 #else 1309 rmap_i = rmap[is_no]; 1310 cmap_i = cmap[is_no]; 1311 #endif 1312 mat = (Mat_SeqBAIJ*)submats[is_no]->data; 1313 imat_ilen = mat->ilen; 1314 imat_j = mat->j; 1315 imat_i = mat->i; 1316 imat_a = mat->a; 1317 max1 = sbuf1_i[2*j]; 1318 for (k=0; k<max1; k++,ct1++) { 1319 row = sbuf1_i[ct1]; 1320 #if defined (PETSC_USE_CTABLE) 1321 ierr = PetscTableFind(lrow1_grow1,row+1,&row);CHKERRQ(ierr); 1322 row--; 1323 if(row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1324 #else 1325 row = rmap_i[row]; 1326 #endif 1327 ilen = imat_ilen[row]; 1328 mat_i = imat_i[row]; 1329 mat_a = imat_a + mat_i*bs2; 1330 mat_j = imat_j + mat_i; 1331 max2 = rbuf2_i[ct1]; 1332 for (l=0; l<max2; l++,ct2++) { 1333 #if defined (PETSC_USE_CTABLE) 1334 ierr = PetscTableFind(lcol1_gcol1,rbuf3_i[ct2]+1,&tcol);CHKERRQ(ierr); 1335 if (tcol) { 1336 #else 1337 if ((tcol = cmap_i[rbuf3_i[ct2]])) { 1338 #endif 1339 *mat_j++ = tcol - 1; 1340 /* *mat_a++= rbuf4_i[ct2]; */ 1341 ierr = PetscMemcpy(mat_a,rbuf4_i+ct2*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1342 mat_a += bs2; 1343 ilen++; 1344 } 1345 } 1346 imat_ilen[row] = ilen; 1347 } 1348 } 1349 } 1350 } 1351 ierr = PetscFree(r_status4);CHKERRQ(ierr); 1352 ierr = PetscFree(r_waits4);CHKERRQ(ierr); 1353 ierr = MPI_Waitall(nrqr,s_waits4,s_status4);CHKERRQ(ierr); 1354 ierr = PetscFree(s_waits4);CHKERRQ(ierr); 1355 ierr = PetscFree(s_status4);CHKERRQ(ierr); 1356 1357 /* Restore the indices */ 1358 for (i=0; i<ismax; i++) { 1359 ierr = ISRestoreIndices(isrow[i],irow+i);CHKERRQ(ierr); 1360 ierr = ISRestoreIndices(iscol[i],icol+i);CHKERRQ(ierr); 1361 } 1362 1363 /* Destroy allocated memory */ 1364 ierr = PetscFree(irow);CHKERRQ(ierr); 1365 ierr = PetscFree(w1);CHKERRQ(ierr); 1366 ierr = PetscFree(pa);CHKERRQ(ierr); 1367 1368 ierr = PetscFree(sbuf1);CHKERRQ(ierr); 1369 ierr = PetscFree(rbuf2);CHKERRQ(ierr); 1370 for (i=0; i<nrqr; ++i) { 1371 ierr = PetscFree(sbuf2[i]);CHKERRQ(ierr); 1372 } 1373 for (i=0; i<nrqs; ++i) { 1374 ierr = PetscFree(rbuf3[i]);CHKERRQ(ierr); 1375 ierr = PetscFree(rbuf4[i]);CHKERRQ(ierr); 1376 } 1377 1378 ierr = PetscFree(sbuf2);CHKERRQ(ierr); 1379 ierr = PetscFree(rbuf3);CHKERRQ(ierr); 1380 ierr = PetscFree(rbuf4);CHKERRQ(ierr); 1381 ierr = PetscFree(sbuf_aj[0]);CHKERRQ(ierr); 1382 ierr = PetscFree(sbuf_aj);CHKERRQ(ierr); 1383 ierr = PetscFree(sbuf_aa[0]);CHKERRQ(ierr); 1384 ierr = PetscFree(sbuf_aa);CHKERRQ(ierr); 1385 1386 #if defined (PETSC_USE_CTABLE) 1387 for (i=0; i<ismax; i++){ 1388 ierr = PetscTableDelete(rowmaps[i]);CHKERRQ(ierr); 1389 ierr = PetscTableDelete(colmaps[i]);CHKERRQ(ierr); 1390 } 1391 ierr = PetscFree(colmaps);CHKERRQ(ierr); 1392 ierr = PetscFree(rowmaps);CHKERRQ(ierr); 1393 #else 1394 ierr = PetscFree(rmap);CHKERRQ(ierr); 1395 ierr = PetscFree(cmap);CHKERRQ(ierr); 1396 #endif 1397 ierr = PetscFree(lens);CHKERRQ(ierr); 1398 1399 for (i=0; i<ismax; i++) { 1400 ierr = MatAssemblyBegin(submats[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1401 ierr = MatAssemblyEnd(submats[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1402 } 1403 1404 PetscFunctionReturn(0); 1405 } 1406 1407