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