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->cmap->N, bs=C->rmap->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 const PetscInt **idx,*idx_i; 68 PetscInt *n,*w3,*w4,*rtable,**data,len; 69 PetscErrorCode ierr; 70 PetscMPIInt size,rank,tag1,tag2,*w2,*w1,nrqr; 71 PetscInt Mbs,i,j,k,**rbuf,row,proc,nrqs,msz,**outdat,**ptr; 72 PetscInt *ctr,*pa,*tmp,*isz,*isz1,**xdata,**rbuf2; 73 PetscMPIInt *onodes1,*olengths1,*onodes2,*olengths2; 74 PetscBT *table; 75 MPI_Comm comm; 76 MPI_Request *s_waits1,*r_waits1,*s_waits2,*r_waits2; 77 MPI_Status *s_status,*recv_status; 78 79 PetscFunctionBegin; 80 comm = ((PetscObject)C)->comm; 81 size = c->size; 82 rank = c->rank; 83 Mbs = c->Mbs; 84 85 ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr); 86 ierr = PetscObjectGetNewTag((PetscObject)C,&tag2);CHKERRQ(ierr); 87 88 ierr = PetscMalloc((imax+1)*sizeof(PetscInt*)+ (imax + Mbs)*sizeof(PetscInt),&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->rangebs[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 ierr = PetscMalloc(2*size*sizeof(PetscInt*) + (size+msz)*sizeof(PetscInt),&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 ierr = PetscMalloc((imax)*(sizeof(PetscBT) + sizeof(PetscInt*)+ sizeof(PetscInt)) + 188 (Mbs)*imax*sizeof(PetscInt) + (Mbs/PETSC_BITS_PER_BYTE+1)*imax*sizeof(char),&table);CHKERRQ(ierr); 189 ierr = PetscMemzero(table,(imax)*(sizeof(PetscBT) + sizeof(PetscInt*)+ sizeof(PetscInt)) + 190 (Mbs)*imax*sizeof(PetscInt) + (Mbs/PETSC_BITS_PER_BYTE+1)*imax*sizeof(char));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 if (nrqr) {ierr = MPI_Waitall(nrqr,r_waits1,recv_status);CHKERRQ(ierr);} 261 262 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&s_status);CHKERRQ(ierr); 263 if (nrqs) {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 ierr = PetscFree(rw1);CHKERRQ(ierr); 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 if (nrqr) {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->rstartbs; 387 cstart = c->cstartbs; 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 separately 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->rstartbs; 455 cstart = c->cstartbs; 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 ierr = PetscInfo3(C,"Allocated %D bytes, required %D, no of mallocs = %D\n",mem_estimate,ct3,no_malloc);CHKERRQ(ierr); 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->cmap->N,bs=C->rmap->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 = ISCompressIndicesGeneral(N,bs,ismax,isrow,isrow_new);CHKERRQ(ierr); 575 ierr = ISCompressIndicesGeneral(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,((PetscObject)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 = PetscMPIIntCast( (PetscInt)(((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 const PetscInt **irow,**icol,*irow_i; 633 PetscInt *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; 640 PetscInt ctr_j,*sbuf1_j,*sbuf_aj_i,*rbuf1_i,kmax,*lens_i; 641 PetscInt bs=C->rmap->bs,bs2=c->bs2,*a_j=a->j,*b_j=b->j,*cworkA,*cworkB; 642 PetscInt cstart = c->cstartbs,nzA,nzB,*a_i=a->i,*b_i=b->i,imark; 643 PetscInt *bmap = c->garray,ctmp,rstart=c->rstartbs; 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 = ((PetscObject)C)->comm; 663 tag0 = ((PetscObject)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 #if defined(PETSC_USE_CTABLE) 673 ierr = PetscMalloc((2*ismax+1)*(sizeof(PetscInt*)+ sizeof(PetscInt)),&irow);CHKERRQ(ierr); 674 #else 675 ierr = PetscMalloc((2*ismax+1)*(sizeof(PetscInt*)+ sizeof(PetscInt)) + (Mbs+1)*sizeof(PetscInt),&irow);CHKERRQ(ierr); 676 #endif 677 678 icol = irow + ismax; 679 nrow = (PetscInt*)(icol + ismax); 680 ncol = nrow + ismax; 681 #if !defined (PETSC_USE_CTABLE) 682 rtable = ncol + ismax; 683 /* Create hash table for the mapping :row -> proc*/ 684 for (i=0,j=0; i<size; i++) { 685 jmax = c->rowners[i+1]; 686 for (; j<jmax; j++) { 687 rtable[j] = i; 688 } 689 } 690 #endif 691 692 for (i=0; i<ismax; i++) { 693 ierr = ISGetIndices(isrow[i],&irow[i]);CHKERRQ(ierr); 694 ierr = ISGetIndices(iscol[i],&icol[i]);CHKERRQ(ierr); 695 ierr = ISGetLocalSize(isrow[i],&nrow[i]);CHKERRQ(ierr); 696 ierr = ISGetLocalSize(iscol[i],&ncol[i]);CHKERRQ(ierr); 697 } 698 699 /* evaluate communication - mesg to who,length of mesg,and buffer space 700 required. Based on this, buffers are allocated, and data copied into them*/ 701 ierr = PetscMalloc(size*2*sizeof(PetscInt) + size*2*sizeof(PetscMPIInt),&w1);CHKERRQ(ierr); /* mesg size */ 702 w2 = w1 + size; /* if w2[i] marked, then a message to proc i*/ 703 w3 = (PetscInt*)(w2 + size); /* no of IS that needs to be sent to proc i */ 704 w4 = w3 + size; /* temp work space used in determining w1, w2, w3 */ 705 ierr = PetscMemzero(w1,size*sizeof(PetscInt)+2*size*sizeof(PetscMPIInt));CHKERRQ(ierr); /* initialise work vector*/ 706 for (i=0; i<ismax; i++) { 707 ierr = PetscMemzero(w4,size*sizeof(PetscInt));CHKERRQ(ierr); /* initialise work vector*/ 708 jmax = nrow[i]; 709 irow_i = irow[i]; 710 for (j=0; j<jmax; j++) { 711 row = irow_i[j]; 712 #if defined (PETSC_USE_CTABLE) 713 ierr = PetscGetProc(row,size,c->rangebs,&proc);CHKERRQ(ierr); 714 #else 715 proc = rtable[row]; 716 #endif 717 w4[proc]++; 718 } 719 for (j=0; j<size; j++) { 720 if (w4[j]) { w1[j] += w4[j]; w3[j]++;} 721 } 722 } 723 724 nrqs = 0; /* no of outgoing messages */ 725 msz = 0; /* total mesg length for all proc */ 726 w1[rank] = 0; /* no mesg sent to intself */ 727 w3[rank] = 0; 728 for (i=0; i<size; i++) { 729 if (w1[i]) { w2[i] = 1; nrqs++;} /* there exists a message to proc i */ 730 } 731 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt),&pa);CHKERRQ(ierr); /*(proc -array)*/ 732 for (i=0,j=0; i<size; i++) { 733 if (w1[i]) { pa[j] = i; j++; } 734 } 735 736 /* Each message would have a header = 1 + 2*(no of IS) + data */ 737 for (i=0; i<nrqs; i++) { 738 j = pa[i]; 739 w1[j] += w2[j] + 2* w3[j]; 740 msz += w1[j]; 741 } 742 743 /* Determine the number of messages to expect, their lengths, from from-ids */ 744 ierr = PetscGatherNumberOfMessages(comm,w2,w1,&nrqr);CHKERRQ(ierr); 745 ierr = PetscGatherMessageLengths(comm,nrqs,nrqr,w1,&onodes1,&olengths1);CHKERRQ(ierr); 746 747 /* Now post the Irecvs corresponding to these messages */ 748 ierr = PetscPostIrecvInt(comm,tag0,nrqr,onodes1,olengths1,&rbuf1,&r_waits1);CHKERRQ(ierr); 749 750 ierr = PetscFree(onodes1);CHKERRQ(ierr); 751 ierr = PetscFree(olengths1);CHKERRQ(ierr); 752 753 /* Allocate Memory for outgoing messages */ 754 ierr = PetscMalloc(2*size*sizeof(PetscInt*) + 2*msz*sizeof(PetscInt) + size*sizeof(PetscInt),&sbuf1);CHKERRQ(ierr); 755 ptr = sbuf1 + size; /* Pointers to the data in outgoing buffers */ 756 ierr = PetscMemzero(sbuf1,2*size*sizeof(PetscInt*));CHKERRQ(ierr); 757 /* allocate memory for outgoing data + buf to receive the first reply */ 758 tmp = (PetscInt*)(ptr + size); 759 ctr = tmp + 2*msz; 760 761 { 762 PetscInt *iptr = tmp,ict = 0; 763 for (i=0; i<nrqs; i++) { 764 j = pa[i]; 765 iptr += ict; 766 sbuf1[j] = iptr; 767 ict = w1[j]; 768 } 769 } 770 771 /* Form the outgoing messages */ 772 /* Initialise the header space */ 773 for (i=0; i<nrqs; i++) { 774 j = pa[i]; 775 sbuf1[j][0] = 0; 776 ierr = PetscMemzero(sbuf1[j]+1,2*w3[j]*sizeof(PetscInt));CHKERRQ(ierr); 777 ptr[j] = sbuf1[j] + 2*w3[j] + 1; 778 } 779 780 /* Parse the isrow and copy data into outbuf */ 781 for (i=0; i<ismax; i++) { 782 ierr = PetscMemzero(ctr,size*sizeof(PetscInt));CHKERRQ(ierr); 783 irow_i = irow[i]; 784 jmax = nrow[i]; 785 for (j=0; j<jmax; j++) { /* parse the indices of each IS */ 786 row = irow_i[j]; 787 #if defined (PETSC_USE_CTABLE) 788 ierr = PetscGetProc(row,size,c->rangebs,&proc);CHKERRQ(ierr); 789 #else 790 proc = rtable[row]; 791 #endif 792 if (proc != rank) { /* copy to the outgoing buf*/ 793 ctr[proc]++; 794 *ptr[proc] = row; 795 ptr[proc]++; 796 } 797 } 798 /* Update the headers for the current IS */ 799 for (j=0; j<size; j++) { /* Can Optimise this loop too */ 800 if ((ctr_j = ctr[j])) { 801 sbuf1_j = sbuf1[j]; 802 k = ++sbuf1_j[0]; 803 sbuf1_j[2*k] = ctr_j; 804 sbuf1_j[2*k-1] = i; 805 } 806 } 807 } 808 809 /* Now post the sends */ 810 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&s_waits1);CHKERRQ(ierr); 811 for (i=0; i<nrqs; ++i) { 812 j = pa[i]; 813 ierr = MPI_Isend(sbuf1[j],w1[j],MPIU_INT,j,tag0,comm,s_waits1+i);CHKERRQ(ierr); 814 } 815 816 /* Post Recieves to capture the buffer size */ 817 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits2);CHKERRQ(ierr); 818 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt*),&rbuf2);CHKERRQ(ierr); 819 rbuf2[0] = tmp + msz; 820 for (i=1; i<nrqs; ++i) { 821 j = pa[i]; 822 rbuf2[i] = rbuf2[i-1]+w1[pa[i-1]]; 823 } 824 for (i=0; i<nrqs; ++i) { 825 j = pa[i]; 826 ierr = MPI_Irecv(rbuf2[i],w1[j],MPIU_INT,j,tag1,comm,r_waits2+i);CHKERRQ(ierr); 827 } 828 829 /* Send to other procs the buf size they should allocate */ 830 831 /* Receive messages*/ 832 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits2);CHKERRQ(ierr); 833 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&r_status1);CHKERRQ(ierr); 834 ierr = PetscMalloc(2*nrqr*sizeof(PetscInt) + (nrqr+1)*sizeof(PetscInt*),&sbuf2);CHKERRQ(ierr); 835 req_size = (PetscInt*)(sbuf2 + nrqr); 836 req_source = req_size + nrqr; 837 838 { 839 Mat_SeqBAIJ *sA = (Mat_SeqBAIJ*)c->A->data,*sB = (Mat_SeqBAIJ*)c->B->data; 840 PetscInt *sAi = sA->i,*sBi = sB->i,id,*sbuf2_i; 841 842 for (i=0; i<nrqr; ++i) { 843 ierr = MPI_Waitany(nrqr,r_waits1,&idex,r_status1+i);CHKERRQ(ierr); 844 req_size[idex] = 0; 845 rbuf1_i = rbuf1[idex]; 846 start = 2*rbuf1_i[0] + 1; 847 ierr = MPI_Get_count(r_status1+i,MPIU_INT,&end);CHKERRQ(ierr); 848 ierr = PetscMalloc(end*sizeof(PetscInt),&sbuf2[idex]);CHKERRQ(ierr); 849 sbuf2_i = sbuf2[idex]; 850 for (j=start; j<end; j++) { 851 id = rbuf1_i[j] - rstart; 852 ncols = sAi[id+1] - sAi[id] + sBi[id+1] - sBi[id]; 853 sbuf2_i[j] = ncols; 854 req_size[idex] += ncols; 855 } 856 req_source[idex] = r_status1[i].MPI_SOURCE; 857 /* form the header */ 858 sbuf2_i[0] = req_size[idex]; 859 for (j=1; j<start; j++) { sbuf2_i[j] = rbuf1_i[j]; } 860 ierr = MPI_Isend(sbuf2_i,end,MPIU_INT,req_source[idex],tag1,comm,s_waits2+i);CHKERRQ(ierr); 861 } 862 } 863 ierr = PetscFree(r_status1);CHKERRQ(ierr); 864 ierr = PetscFree(r_waits1);CHKERRQ(ierr); 865 866 /* recv buffer sizes */ 867 /* Receive messages*/ 868 869 ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt*),&rbuf3);CHKERRQ(ierr); 870 ierr = PetscMalloc((nrqs+1)*sizeof(MatScalar*),&rbuf4);CHKERRQ(ierr); 871 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits3);CHKERRQ(ierr); 872 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Request),&r_waits4);CHKERRQ(ierr); 873 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status2);CHKERRQ(ierr); 874 875 for (i=0; i<nrqs; ++i) { 876 ierr = MPI_Waitany(nrqs,r_waits2,&idex,r_status2+i);CHKERRQ(ierr); 877 ierr = PetscMalloc(rbuf2[idex][0]*sizeof(PetscInt),&rbuf3[idex]);CHKERRQ(ierr); 878 ierr = PetscMalloc(rbuf2[idex][0]*bs2*sizeof(MatScalar),&rbuf4[idex]);CHKERRQ(ierr); 879 ierr = MPI_Irecv(rbuf3[idex],rbuf2[idex][0],MPIU_INT,r_status2[i].MPI_SOURCE,tag2,comm,r_waits3+idex);CHKERRQ(ierr); 880 ierr = MPI_Irecv(rbuf4[idex],rbuf2[idex][0]*bs2,MPIU_MATSCALAR,r_status2[i].MPI_SOURCE,tag3,comm,r_waits4+idex);CHKERRQ(ierr); 881 } 882 ierr = PetscFree(r_status2);CHKERRQ(ierr); 883 ierr = PetscFree(r_waits2);CHKERRQ(ierr); 884 885 /* Wait on sends1 and sends2 */ 886 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&s_status1);CHKERRQ(ierr); 887 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status2);CHKERRQ(ierr); 888 889 if (nrqs) {ierr = MPI_Waitall(nrqs,s_waits1,s_status1);CHKERRQ(ierr);} 890 if (nrqr) {ierr = MPI_Waitall(nrqr,s_waits2,s_status2);CHKERRQ(ierr);} 891 ierr = PetscFree(s_status1);CHKERRQ(ierr); 892 ierr = PetscFree(s_status2);CHKERRQ(ierr); 893 ierr = PetscFree(s_waits1);CHKERRQ(ierr); 894 ierr = PetscFree(s_waits2);CHKERRQ(ierr); 895 896 /* Now allocate buffers for a->j, and send them off */ 897 ierr = PetscMalloc((nrqr+1)*sizeof(PetscInt*),&sbuf_aj);CHKERRQ(ierr); 898 for (i=0,j=0; i<nrqr; i++) j += req_size[i]; 899 ierr = PetscMalloc((j+1)*sizeof(PetscInt),&sbuf_aj[0]);CHKERRQ(ierr); 900 for (i=1; i<nrqr; i++) sbuf_aj[i] = sbuf_aj[i-1] + req_size[i-1]; 901 902 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits3);CHKERRQ(ierr); 903 { 904 for (i=0; i<nrqr; i++) { 905 rbuf1_i = rbuf1[i]; 906 sbuf_aj_i = sbuf_aj[i]; 907 ct1 = 2*rbuf1_i[0] + 1; 908 ct2 = 0; 909 for (j=1,max1=rbuf1_i[0]; j<=max1; j++) { 910 kmax = rbuf1[i][2*j]; 911 for (k=0; k<kmax; k++,ct1++) { 912 row = rbuf1_i[ct1] - rstart; 913 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 914 ncols = nzA + nzB; 915 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 916 917 /* load the column indices for this row into cols*/ 918 cols = sbuf_aj_i + ct2; 919 for (l=0; l<nzB; l++) { 920 if ((ctmp = bmap[cworkB[l]]) < cstart) cols[l] = ctmp; 921 else break; 922 } 923 imark = l; 924 for (l=0; l<nzA; l++) cols[imark+l] = cstart + cworkA[l]; 925 for (l=imark; l<nzB; l++) cols[nzA+l] = bmap[cworkB[l]]; 926 ct2 += ncols; 927 } 928 } 929 ierr = MPI_Isend(sbuf_aj_i,req_size[i],MPIU_INT,req_source[i],tag2,comm,s_waits3+i);CHKERRQ(ierr); 930 } 931 } 932 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status3);CHKERRQ(ierr); 933 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status3);CHKERRQ(ierr); 934 935 /* Allocate buffers for a->a, and send them off */ 936 ierr = PetscMalloc((nrqr+1)*sizeof(MatScalar *),&sbuf_aa);CHKERRQ(ierr); 937 for (i=0,j=0; i<nrqr; i++) j += req_size[i]; 938 ierr = PetscMalloc((j+1)*bs2*sizeof(MatScalar),&sbuf_aa[0]);CHKERRQ(ierr); 939 for (i=1; i<nrqr; i++) sbuf_aa[i] = sbuf_aa[i-1] + req_size[i-1]*bs2; 940 941 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&s_waits4);CHKERRQ(ierr); 942 { 943 for (i=0; i<nrqr; i++) { 944 rbuf1_i = rbuf1[i]; 945 sbuf_aa_i = sbuf_aa[i]; 946 ct1 = 2*rbuf1_i[0]+1; 947 ct2 = 0; 948 for (j=1,max1=rbuf1_i[0]; j<=max1; j++) { 949 kmax = rbuf1_i[2*j]; 950 for (k=0; k<kmax; k++,ct1++) { 951 row = rbuf1_i[ct1] - rstart; 952 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 953 ncols = nzA + nzB; 954 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 955 vworkA = a_a + a_i[row]*bs2; vworkB = b_a + b_i[row]*bs2; 956 957 /* load the column values for this row into vals*/ 958 vals = sbuf_aa_i+ct2*bs2; 959 for (l=0; l<nzB; l++) { 960 if ((bmap[cworkB[l]]) < cstart) { 961 ierr = PetscMemcpy(vals+l*bs2,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 962 } 963 else break; 964 } 965 imark = l; 966 for (l=0; l<nzA; l++) { 967 ierr = PetscMemcpy(vals+(imark+l)*bs2,vworkA+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 968 } 969 for (l=imark; l<nzB; l++) { 970 ierr = PetscMemcpy(vals+(nzA+l)*bs2,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 971 } 972 ct2 += ncols; 973 } 974 } 975 ierr = MPI_Isend(sbuf_aa_i,req_size[i]*bs2,MPIU_MATSCALAR,req_source[i],tag3,comm,s_waits4+i);CHKERRQ(ierr); 976 } 977 } 978 ierr = PetscMalloc((nrqs+1)*sizeof(MPI_Status),&r_status4);CHKERRQ(ierr); 979 ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Status),&s_status4);CHKERRQ(ierr); 980 ierr = PetscFree(rbuf1);CHKERRQ(ierr); 981 982 /* Form the matrix */ 983 /* create col map */ 984 { 985 const PetscInt *icol_i; 986 #if defined (PETSC_USE_CTABLE) 987 /* Create row map*/ 988 ierr = PetscMalloc((1+ismax)*sizeof(PetscTable),&colmaps);CHKERRQ(ierr); 989 for (i=0; i<ismax; i++) { 990 ierr = PetscTableCreate(ncol[i]+1,&colmaps[i]);CHKERRQ(ierr); 991 } 992 #else 993 ierr = PetscMalloc((1+ismax)*sizeof(PetscInt*)+ ismax*c->Nbs*sizeof(PetscInt),&cmap);CHKERRQ(ierr); 994 cmap[0] = (PetscInt*)(cmap + ismax); 995 ierr = PetscMemzero(cmap[0],(1+ismax*c->Nbs)*sizeof(PetscInt));CHKERRQ(ierr); 996 for (i=1; i<ismax; i++) { cmap[i] = cmap[i-1] + c->Nbs; } 997 #endif 998 for (i=0; i<ismax; i++) { 999 jmax = ncol[i]; 1000 icol_i = icol[i]; 1001 #if defined (PETSC_USE_CTABLE) 1002 lcol1_gcol1 = colmaps[i]; 1003 for (j=0; j<jmax; j++) { 1004 ierr = PetscTableAdd(lcol1_gcol1,icol_i[j]+1,j+1);CHKERRQ(ierr); 1005 } 1006 #else 1007 cmap_i = cmap[i]; 1008 for (j=0; j<jmax; j++) { 1009 cmap_i[icol_i[j]] = j+1; 1010 } 1011 #endif 1012 } 1013 } 1014 1015 /* Create lens which is required for MatCreate... */ 1016 for (i=0,j=0; i<ismax; i++) { j += nrow[i]; } 1017 ierr = PetscMalloc((1+ismax)*sizeof(PetscInt*)+ j*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1018 lens[0] = (PetscInt*)(lens + ismax); 1019 ierr = PetscMemzero(lens[0],j*sizeof(PetscInt));CHKERRQ(ierr); 1020 for (i=1; i<ismax; i++) { lens[i] = lens[i-1] + nrow[i-1]; } 1021 1022 /* Update lens from local data */ 1023 for (i=0; i<ismax; i++) { 1024 jmax = nrow[i]; 1025 #if defined (PETSC_USE_CTABLE) 1026 lcol1_gcol1 = colmaps[i]; 1027 #else 1028 cmap_i = cmap[i]; 1029 #endif 1030 irow_i = irow[i]; 1031 lens_i = lens[i]; 1032 for (j=0; j<jmax; j++) { 1033 row = irow_i[j]; 1034 #if defined (PETSC_USE_CTABLE) 1035 ierr = PetscGetProc(row,size,c->rangebs,&proc);CHKERRQ(ierr); 1036 #else 1037 proc = rtable[row]; 1038 #endif 1039 if (proc == rank) { 1040 /* Get indices from matA and then from matB */ 1041 row = row - rstart; 1042 nzA = a_i[row+1] - a_i[row]; nzB = b_i[row+1] - b_i[row]; 1043 cworkA = a_j + a_i[row]; cworkB = b_j + b_i[row]; 1044 #if defined (PETSC_USE_CTABLE) 1045 for (k=0; k<nzA; k++) { 1046 ierr = PetscTableFind(lcol1_gcol1,cstart+cworkA[k]+1,&tt);CHKERRQ(ierr); 1047 if (tt) { lens_i[j]++; } 1048 } 1049 for (k=0; k<nzB; k++) { 1050 ierr = PetscTableFind(lcol1_gcol1,bmap[cworkB[k]]+1,&tt);CHKERRQ(ierr); 1051 if (tt) { lens_i[j]++; } 1052 } 1053 #else 1054 for (k=0; k<nzA; k++) { 1055 if (cmap_i[cstart + cworkA[k]]) { lens_i[j]++; } 1056 } 1057 for (k=0; k<nzB; k++) { 1058 if (cmap_i[bmap[cworkB[k]]]) { lens_i[j]++; } 1059 } 1060 #endif 1061 } 1062 } 1063 } 1064 #if defined (PETSC_USE_CTABLE) 1065 /* Create row map*/ 1066 ierr = PetscMalloc((1+ismax)*sizeof(PetscTable),&rowmaps);CHKERRQ(ierr); 1067 for (i=0; i<ismax; i++){ 1068 ierr = PetscTableCreate(nrow[i]+1,&rowmaps[i]);CHKERRQ(ierr); 1069 } 1070 #else 1071 /* Create row map*/ 1072 ierr = PetscMalloc((1+ismax)*sizeof(PetscInt*)+ ismax*Mbs*sizeof(PetscInt),&rmap);CHKERRQ(ierr); 1073 rmap[0] = (PetscInt*)(rmap + ismax); 1074 ierr = PetscMemzero(rmap[0],ismax*Mbs*sizeof(PetscInt));CHKERRQ(ierr); 1075 for (i=1; i<ismax; i++) { rmap[i] = rmap[i-1] + Mbs;} 1076 #endif 1077 for (i=0; i<ismax; i++) { 1078 irow_i = irow[i]; 1079 jmax = nrow[i]; 1080 #if defined (PETSC_USE_CTABLE) 1081 lrow1_grow1 = rowmaps[i]; 1082 for (j=0; j<jmax; j++) { 1083 ierr = PetscTableAdd(lrow1_grow1,irow_i[j]+1,j+1);CHKERRQ(ierr); 1084 } 1085 #else 1086 rmap_i = rmap[i]; 1087 for (j=0; j<jmax; j++) { 1088 rmap_i[irow_i[j]] = j; 1089 } 1090 #endif 1091 } 1092 1093 /* Update lens from offproc data */ 1094 { 1095 PetscInt *rbuf2_i,*rbuf3_i,*sbuf1_i; 1096 PetscMPIInt ii; 1097 1098 for (tmp2=0; tmp2<nrqs; tmp2++) { 1099 ierr = MPI_Waitany(nrqs,r_waits3,&ii,r_status3+tmp2);CHKERRQ(ierr); 1100 idex = pa[ii]; 1101 sbuf1_i = sbuf1[idex]; 1102 jmax = sbuf1_i[0]; 1103 ct1 = 2*jmax+1; 1104 ct2 = 0; 1105 rbuf2_i = rbuf2[ii]; 1106 rbuf3_i = rbuf3[ii]; 1107 for (j=1; j<=jmax; j++) { 1108 is_no = sbuf1_i[2*j-1]; 1109 max1 = sbuf1_i[2*j]; 1110 lens_i = lens[is_no]; 1111 #if defined (PETSC_USE_CTABLE) 1112 lcol1_gcol1 = colmaps[is_no]; 1113 lrow1_grow1 = rowmaps[is_no]; 1114 #else 1115 cmap_i = cmap[is_no]; 1116 rmap_i = rmap[is_no]; 1117 #endif 1118 for (k=0; k<max1; k++,ct1++) { 1119 #if defined (PETSC_USE_CTABLE) 1120 ierr = PetscTableFind(lrow1_grow1,sbuf1_i[ct1]+1,&row);CHKERRQ(ierr); 1121 row--; 1122 if (row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1123 #else 1124 row = rmap_i[sbuf1_i[ct1]]; /* the val in the new matrix to be */ 1125 #endif 1126 max2 = rbuf2_i[ct1]; 1127 for (l=0; l<max2; l++,ct2++) { 1128 #if defined (PETSC_USE_CTABLE) 1129 ierr = PetscTableFind(lcol1_gcol1,rbuf3_i[ct2]+1,&tt);CHKERRQ(ierr); 1130 if (tt) { 1131 lens_i[row]++; 1132 } 1133 #else 1134 if (cmap_i[rbuf3_i[ct2]]) { 1135 lens_i[row]++; 1136 } 1137 #endif 1138 } 1139 } 1140 } 1141 } 1142 } 1143 ierr = PetscFree(r_status3);CHKERRQ(ierr); 1144 ierr = PetscFree(r_waits3);CHKERRQ(ierr); 1145 if (nrqr) {ierr = MPI_Waitall(nrqr,s_waits3,s_status3);CHKERRQ(ierr);} 1146 ierr = PetscFree(s_status3);CHKERRQ(ierr); 1147 ierr = PetscFree(s_waits3);CHKERRQ(ierr); 1148 1149 /* Create the submatrices */ 1150 if (scall == MAT_REUSE_MATRIX) { 1151 /* 1152 Assumes new rows are same length as the old rows, hence bug! 1153 */ 1154 for (i=0; i<ismax; i++) { 1155 mat = (Mat_SeqBAIJ *)(submats[i]->data); 1156 if ((mat->mbs != nrow[i]) || (mat->nbs != ncol[i] || C->rmap->bs != bs)) { 1157 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1158 } 1159 ierr = PetscMemcmp(mat->ilen,lens[i],mat->mbs *sizeof(PetscInt),&flag);CHKERRQ(ierr); 1160 if (!flag) { 1161 SETERRQ(PETSC_ERR_ARG_INCOMP,"Cannot reuse matrix. wrong no of nonzeros"); 1162 } 1163 /* Initial matrix as if empty */ 1164 ierr = PetscMemzero(mat->ilen,mat->mbs*sizeof(PetscInt));CHKERRQ(ierr); 1165 submats[i]->factor = C->factor; 1166 } 1167 } else { 1168 for (i=0; i<ismax; i++) { 1169 ierr = MatCreate(PETSC_COMM_SELF,submats+i);CHKERRQ(ierr); 1170 ierr = MatSetSizes(submats[i],nrow[i]*bs,ncol[i]*bs,nrow[i]*bs,ncol[i]*bs);CHKERRQ(ierr); 1171 ierr = MatSetType(submats[i],((PetscObject)A)->type_name);CHKERRQ(ierr); 1172 ierr = MatSeqBAIJSetPreallocation(submats[i],C->rmap->bs,0,lens[i]);CHKERRQ(ierr); 1173 ierr = MatSeqSBAIJSetPreallocation(submats[i],C->rmap->bs,0,lens[i]);CHKERRQ(ierr); 1174 } 1175 } 1176 1177 /* Assemble the matrices */ 1178 /* First assemble the local rows */ 1179 { 1180 PetscInt ilen_row,*imat_ilen,*imat_j,*imat_i; 1181 MatScalar *imat_a; 1182 1183 for (i=0; i<ismax; i++) { 1184 mat = (Mat_SeqBAIJ*)submats[i]->data; 1185 imat_ilen = mat->ilen; 1186 imat_j = mat->j; 1187 imat_i = mat->i; 1188 imat_a = mat->a; 1189 1190 #if defined (PETSC_USE_CTABLE) 1191 lcol1_gcol1 = colmaps[i]; 1192 lrow1_grow1 = rowmaps[i]; 1193 #else 1194 cmap_i = cmap[i]; 1195 rmap_i = rmap[i]; 1196 #endif 1197 irow_i = irow[i]; 1198 jmax = nrow[i]; 1199 for (j=0; j<jmax; j++) { 1200 row = irow_i[j]; 1201 #if defined (PETSC_USE_CTABLE) 1202 ierr = PetscGetProc(row,size,c->rangebs,&proc);CHKERRQ(ierr); 1203 #else 1204 proc = rtable[row]; 1205 #endif 1206 if (proc == rank) { 1207 row = row - rstart; 1208 nzA = a_i[row+1] - a_i[row]; 1209 nzB = b_i[row+1] - b_i[row]; 1210 cworkA = a_j + a_i[row]; 1211 cworkB = b_j + b_i[row]; 1212 vworkA = a_a + a_i[row]*bs2; 1213 vworkB = b_a + b_i[row]*bs2; 1214 #if defined (PETSC_USE_CTABLE) 1215 ierr = PetscTableFind(lrow1_grow1,row+rstart+1,&row);CHKERRQ(ierr); 1216 row--; 1217 if (row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1218 #else 1219 row = rmap_i[row + rstart]; 1220 #endif 1221 mat_i = imat_i[row]; 1222 mat_a = imat_a + mat_i*bs2; 1223 mat_j = imat_j + mat_i; 1224 ilen_row = imat_ilen[row]; 1225 1226 /* load the column indices for this row into cols*/ 1227 for (l=0; l<nzB; l++) { 1228 if ((ctmp = bmap[cworkB[l]]) < cstart) { 1229 #if defined (PETSC_USE_CTABLE) 1230 ierr = PetscTableFind(lcol1_gcol1,ctmp+1,&tcol);CHKERRQ(ierr); 1231 if (tcol) { 1232 #else 1233 if ((tcol = cmap_i[ctmp])) { 1234 #endif 1235 *mat_j++ = tcol - 1; 1236 ierr = PetscMemcpy(mat_a,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1237 mat_a += bs2; 1238 ilen_row++; 1239 } 1240 } else break; 1241 } 1242 imark = l; 1243 for (l=0; l<nzA; l++) { 1244 #if defined (PETSC_USE_CTABLE) 1245 ierr = PetscTableFind(lcol1_gcol1,cstart+cworkA[l]+1,&tcol);CHKERRQ(ierr); 1246 if (tcol) { 1247 #else 1248 if ((tcol = cmap_i[cstart + cworkA[l]])) { 1249 #endif 1250 *mat_j++ = tcol - 1; 1251 ierr = PetscMemcpy(mat_a,vworkA+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1252 mat_a += bs2; 1253 ilen_row++; 1254 } 1255 } 1256 for (l=imark; l<nzB; l++) { 1257 #if defined (PETSC_USE_CTABLE) 1258 ierr = PetscTableFind(lcol1_gcol1,bmap[cworkB[l]]+1,&tcol);CHKERRQ(ierr); 1259 if (tcol) { 1260 #else 1261 if ((tcol = cmap_i[bmap[cworkB[l]]])) { 1262 #endif 1263 *mat_j++ = tcol - 1; 1264 ierr = PetscMemcpy(mat_a,vworkB+l*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1265 mat_a += bs2; 1266 ilen_row++; 1267 } 1268 } 1269 imat_ilen[row] = ilen_row; 1270 } 1271 } 1272 1273 } 1274 } 1275 1276 /* Now assemble the off proc rows*/ 1277 { 1278 PetscInt *sbuf1_i,*rbuf2_i,*rbuf3_i,*imat_ilen,ilen; 1279 PetscInt *imat_j,*imat_i; 1280 MatScalar *imat_a,*rbuf4_i; 1281 PetscMPIInt ii; 1282 1283 for (tmp2=0; tmp2<nrqs; tmp2++) { 1284 ierr = MPI_Waitany(nrqs,r_waits4,&ii,r_status4+tmp2);CHKERRQ(ierr); 1285 idex = pa[ii]; 1286 sbuf1_i = sbuf1[idex]; 1287 jmax = sbuf1_i[0]; 1288 ct1 = 2*jmax + 1; 1289 ct2 = 0; 1290 rbuf2_i = rbuf2[ii]; 1291 rbuf3_i = rbuf3[ii]; 1292 rbuf4_i = rbuf4[ii]; 1293 for (j=1; j<=jmax; j++) { 1294 is_no = sbuf1_i[2*j-1]; 1295 #if defined (PETSC_USE_CTABLE) 1296 lrow1_grow1 = rowmaps[is_no]; 1297 lcol1_gcol1 = colmaps[is_no]; 1298 #else 1299 rmap_i = rmap[is_no]; 1300 cmap_i = cmap[is_no]; 1301 #endif 1302 mat = (Mat_SeqBAIJ*)submats[is_no]->data; 1303 imat_ilen = mat->ilen; 1304 imat_j = mat->j; 1305 imat_i = mat->i; 1306 imat_a = mat->a; 1307 max1 = sbuf1_i[2*j]; 1308 for (k=0; k<max1; k++,ct1++) { 1309 row = sbuf1_i[ct1]; 1310 #if defined (PETSC_USE_CTABLE) 1311 ierr = PetscTableFind(lrow1_grow1,row+1,&row);CHKERRQ(ierr); 1312 row--; 1313 if(row < 0) { SETERRQ(PETSC_ERR_PLIB,"row not found in table"); } 1314 #else 1315 row = rmap_i[row]; 1316 #endif 1317 ilen = imat_ilen[row]; 1318 mat_i = imat_i[row]; 1319 mat_a = imat_a + mat_i*bs2; 1320 mat_j = imat_j + mat_i; 1321 max2 = rbuf2_i[ct1]; 1322 for (l=0; l<max2; l++,ct2++) { 1323 #if defined (PETSC_USE_CTABLE) 1324 ierr = PetscTableFind(lcol1_gcol1,rbuf3_i[ct2]+1,&tcol);CHKERRQ(ierr); 1325 if (tcol) { 1326 #else 1327 if ((tcol = cmap_i[rbuf3_i[ct2]])) { 1328 #endif 1329 *mat_j++ = tcol - 1; 1330 /* *mat_a++= rbuf4_i[ct2]; */ 1331 ierr = PetscMemcpy(mat_a,rbuf4_i+ct2*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1332 mat_a += bs2; 1333 ilen++; 1334 } 1335 } 1336 imat_ilen[row] = ilen; 1337 } 1338 } 1339 } 1340 } 1341 ierr = PetscFree(r_status4);CHKERRQ(ierr); 1342 ierr = PetscFree(r_waits4);CHKERRQ(ierr); 1343 if (nrqr) {ierr = MPI_Waitall(nrqr,s_waits4,s_status4);CHKERRQ(ierr);} 1344 ierr = PetscFree(s_waits4);CHKERRQ(ierr); 1345 ierr = PetscFree(s_status4);CHKERRQ(ierr); 1346 1347 /* Restore the indices */ 1348 for (i=0; i<ismax; i++) { 1349 ierr = ISRestoreIndices(isrow[i],irow+i);CHKERRQ(ierr); 1350 ierr = ISRestoreIndices(iscol[i],icol+i);CHKERRQ(ierr); 1351 } 1352 1353 /* Destroy allocated memory */ 1354 ierr = PetscFree(irow);CHKERRQ(ierr); 1355 ierr = PetscFree(w1);CHKERRQ(ierr); 1356 ierr = PetscFree(pa);CHKERRQ(ierr); 1357 1358 ierr = PetscFree(sbuf1);CHKERRQ(ierr); 1359 ierr = PetscFree(rbuf2);CHKERRQ(ierr); 1360 for (i=0; i<nrqr; ++i) { 1361 ierr = PetscFree(sbuf2[i]);CHKERRQ(ierr); 1362 } 1363 for (i=0; i<nrqs; ++i) { 1364 ierr = PetscFree(rbuf3[i]);CHKERRQ(ierr); 1365 ierr = PetscFree(rbuf4[i]);CHKERRQ(ierr); 1366 } 1367 1368 ierr = PetscFree(sbuf2);CHKERRQ(ierr); 1369 ierr = PetscFree(rbuf3);CHKERRQ(ierr); 1370 ierr = PetscFree(rbuf4);CHKERRQ(ierr); 1371 ierr = PetscFree(sbuf_aj[0]);CHKERRQ(ierr); 1372 ierr = PetscFree(sbuf_aj);CHKERRQ(ierr); 1373 ierr = PetscFree(sbuf_aa[0]);CHKERRQ(ierr); 1374 ierr = PetscFree(sbuf_aa);CHKERRQ(ierr); 1375 1376 #if defined (PETSC_USE_CTABLE) 1377 for (i=0; i<ismax; i++){ 1378 ierr = PetscTableDestroy(rowmaps[i]);CHKERRQ(ierr); 1379 ierr = PetscTableDestroy(colmaps[i]);CHKERRQ(ierr); 1380 } 1381 ierr = PetscFree(colmaps);CHKERRQ(ierr); 1382 ierr = PetscFree(rowmaps);CHKERRQ(ierr); 1383 #else 1384 ierr = PetscFree(rmap);CHKERRQ(ierr); 1385 ierr = PetscFree(cmap);CHKERRQ(ierr); 1386 #endif 1387 ierr = PetscFree(lens);CHKERRQ(ierr); 1388 1389 for (i=0; i<ismax; i++) { 1390 ierr = MatAssemblyBegin(submats[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1391 ierr = MatAssemblyEnd(submats[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1392 } 1393 PetscFunctionReturn(0); 1394 } 1395 1396