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