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