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