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