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>
6: #include <petscbt.h>
8: static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Local(Mat, PetscInt, PetscBT *, PetscInt *, PetscInt **);
9: static PetscErrorCode MatIncreaseOverlap_MPIBAIJ_Receive(Mat, PetscInt, PetscInt **, PetscInt **, PetscInt *);
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;
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: }
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: -----------
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: PetscByte *t_p;
71: PetscFunctionBegin;
72: PetscCall(PetscObjectGetComm((PetscObject)C, &comm));
73: size = c->size;
74: rank = c->rank;
75: Mbs = c->Mbs;
77: PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag1));
78: PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag2));
80: PetscCall(PetscMalloc2(imax, (PetscInt ***)&idx, imax, &n));
82: for (PetscInt i = 0; i < imax; i++) {
83: PetscCall(ISGetIndices(is[i], &idx[i]));
84: PetscCall(ISGetLocalSize(is[i], &n[i]));
85: }
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: }
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: }
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: }
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));
138: /* Now post the Irecvs corresponding to these messages */
139: PetscCall(PetscPostIrecvInt(comm, tag1, nrqr, onodes1, olengths1, &rbuf, &r_waits1));
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: }
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: }
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));
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: }
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;
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: }
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: }
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));
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: }
227: /* Do Local work*/
228: PetscCall(MatIncreaseOverlap_MPIBAIJ_Local(C, imax, table, isz, data));
230: /* Receive messages*/
231: PetscCallMPI(MPI_Waitall(nrqr, r_waits1, MPI_STATUSES_IGNORE));
232: PetscCallMPI(MPI_Waitall(nrqs, s_waits1, MPI_STATUSES_IGNORE));
234: /* Phase 1 sends are complete - deallocate buffers */
235: PetscCall(PetscFree4(outdat, ptr, tmp, ctr));
236: PetscCall(PetscFree4(w1, w2, w3, w4));
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: }
246: /* Send the data back*/
247: /* Do a global reduction to know the buffer space req for incoming messages*/
248: {
249: PetscMPIInt *rw1;
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: }
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));
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: }
271: PetscCall(PetscFree(onodes1));
272: PetscCall(PetscFree(olengths1));
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;
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: }
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: }
307: PetscCall(PetscFree(iscomms));
308: PetscCall(PetscFree(onodes2));
309: PetscCall(PetscFree(olengths2));
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: }
329: /*
330: MatIncreaseOverlap_MPIBAIJ_Local - Called by MatincreaseOverlap, to do
331: the work on the local processor.
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.
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;
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;
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: }
385: /*
386: MatIncreaseOverlap_MPIBAIJ_Receive - Process the received messages,
387: and return the output
389: Input:
390: C - the matrix
391: nrqr - no of messages being processed.
392: rbuf - an array of pointers to the received requests
394: Output:
395: xdata - array of messages to be sent back
396: isz1 - size of each message
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.
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;
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;
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: }
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));
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: }
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, 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];
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));
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;
563: if (scall == MAT_INITIAL_MATRIX) {
564: nstages = ismax / nmax + ((ismax % nmax) ? 1 : 0); /* local nstages */
566: /* Make sure every processor loops through the nstages */
567: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &nstages, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)C)));
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: }
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;
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: }
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: }
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: }
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 PetscDefined(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: PetscInt nzA, nzB, *a_i = a->i, *b_i = b->i, *a_j = a->j, *b_j = b->j, ctmp, imark, *cworkA, *cworkB;
644: PetscScalar *vworkA = NULL, *vworkB = NULL, *a_a = a->a, *b_a = b->a;
645: PetscInt cstart = c->cstartbs, *bmap = c->garray;
646: PetscBool *allrows, *allcolumns;
648: PetscFunctionBegin;
649: PetscCall(PetscObjectGetComm((PetscObject)C, &comm));
650: size = c->size;
651: rank = c->rank;
653: PetscCall(PetscMalloc5(ismax, &row2proc, ismax, &cmap, ismax, &rmap, ismax + 1, &allcolumns, ismax, &allrows));
654: PetscCall(PetscMalloc5(ismax, (PetscInt ***)&irow, ismax, (PetscInt ***)&icol, ismax, &nrow, ismax, &ncol, ismax, &issorted));
656: for (PetscInt i = 0; i < ismax; i++) {
657: PetscCall(ISSorted(iscol[i], &issorted[i]));
658: if (!issorted[i]) iscsorted = issorted[i]; /* columns are not sorted! */
659: PetscCall(ISSorted(isrow[i], &issorted[i]));
661: /* Check for special case: allcolumns */
662: PetscCall(ISIdentity(iscol[i], &colflag));
663: PetscCall(ISGetLocalSize(iscol[i], &ncol[i]));
665: if (colflag && ncol[i] == c->Nbs) {
666: allcolumns[i] = PETSC_TRUE;
667: icol[i] = NULL;
668: } else {
669: allcolumns[i] = PETSC_FALSE;
670: PetscCall(ISGetIndices(iscol[i], &icol[i]));
671: }
673: /* Check for special case: allrows */
674: PetscCall(ISIdentity(isrow[i], &colflag));
675: PetscCall(ISGetLocalSize(isrow[i], &nrow[i]));
676: if (colflag && nrow[i] == c->Mbs) {
677: allrows[i] = PETSC_TRUE;
678: irow[i] = NULL;
679: } else {
680: allrows[i] = PETSC_FALSE;
681: PetscCall(ISGetIndices(isrow[i], &irow[i]));
682: }
683: }
685: if (scall == MAT_REUSE_MATRIX) {
686: /* Assumes new rows are same length as the old rows */
687: for (PetscInt i = 0; i < ismax; i++) {
688: subc = (Mat_SeqBAIJ *)submats[i]->data;
690: PetscCheck(subc->mbs == nrow[i] && subc->nbs == ncol[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong size");
692: /* Initial matrix as if empty */
693: PetscCall(PetscArrayzero(subc->ilen, subc->mbs));
695: /* Initial matrix as if empty */
696: submats[i]->factortype = C->factortype;
698: smat_i = subc->submatis1;
700: nrqs = smat_i->nrqs;
701: nrqr = smat_i->nrqr;
702: rbuf1 = smat_i->rbuf1;
703: rbuf2 = smat_i->rbuf2;
704: rbuf3 = smat_i->rbuf3;
705: req_source2 = smat_i->req_source2;
707: sbuf1 = smat_i->sbuf1;
708: sbuf2 = smat_i->sbuf2;
709: ptr = smat_i->ptr;
710: tmp = smat_i->tmp;
711: ctr = smat_i->ctr;
713: pa = smat_i->pa;
714: req_size = smat_i->req_size;
715: req_source1 = smat_i->req_source1;
717: allcolumns[i] = smat_i->allcolumns;
718: allrows[i] = smat_i->allrows;
719: row2proc[i] = smat_i->row2proc;
720: rmap[i] = smat_i->rmap;
721: cmap[i] = smat_i->cmap;
722: }
724: if (!ismax) { /* Get dummy submatrices and retrieve struct submatis1 */
725: PetscCheck(submats[0], PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "submats are null, cannot reuse");
726: smat_i = (Mat_SubSppt *)submats[0]->data;
728: nrqs = smat_i->nrqs;
729: nrqr = smat_i->nrqr;
730: rbuf1 = smat_i->rbuf1;
731: rbuf2 = smat_i->rbuf2;
732: rbuf3 = smat_i->rbuf3;
733: req_source2 = smat_i->req_source2;
735: sbuf1 = smat_i->sbuf1;
736: sbuf2 = smat_i->sbuf2;
737: ptr = smat_i->ptr;
738: tmp = smat_i->tmp;
739: ctr = smat_i->ctr;
741: pa = smat_i->pa;
742: req_size = smat_i->req_size;
743: req_source1 = smat_i->req_source1;
745: allcolumns[0] = PETSC_FALSE;
746: }
747: } else { /* scall == MAT_INITIAL_MATRIX */
748: /* Get some new tags to keep the communication clean */
749: PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag2));
750: PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag3));
752: /* evaluate communication - mesg to who, length of mesg, and buffer space
753: required. Based on this, buffers are allocated, and data copied into them*/
754: PetscCall(PetscCalloc4(size, &w1, size, &w2, size, &w3, size, &w4)); /* mesg size, initialize work vectors */
756: for (PetscInt i = 0; i < ismax; i++) {
757: jmax = nrow[i];
758: irow_i = irow[i];
760: PetscCall(PetscMalloc1(jmax, &row2proc_i));
761: row2proc[i] = row2proc_i;
763: if (issorted[i]) proc = 0;
764: for (PetscInt j = 0; j < jmax; j++) {
765: if (!issorted[i]) proc = 0;
766: if (allrows[i]) row = j;
767: else row = irow_i[j];
769: while (row >= c->rangebs[proc + 1]) proc++;
770: w4[proc]++;
771: row2proc_i[j] = proc; /* map row index to proc */
772: }
773: for (PetscMPIInt j = 0; j < size; j++) {
774: if (w4[j]) {
775: w1[j] += w4[j];
776: w3[j]++;
777: w4[j] = 0;
778: }
779: }
780: }
782: nrqs = 0; /* no of outgoing messages */
783: msz = 0; /* total mesg length (for all procs) */
784: w1[rank] = 0; /* no mesg sent to self */
785: w3[rank] = 0;
786: for (PetscMPIInt i = 0; i < size; i++) {
787: if (w1[i]) {
788: w2[i] = 1;
789: nrqs++;
790: } /* there exists a message to proc i */
791: }
792: PetscCall(PetscMalloc1(nrqs, &pa)); /*(proc -array)*/
793: for (PetscMPIInt i = 0, j = 0; i < size; i++) {
794: if (w1[i]) pa[j++] = i;
795: }
797: /* Each message would have a header = 1 + 2*(no of IS) + data */
798: for (PetscMPIInt i = 0; i < nrqs; i++) {
799: PetscMPIInt j = pa[i];
800: w1[j] += w2[j] + 2 * w3[j];
801: msz += w1[j];
802: }
803: PetscCall(PetscInfo(0, "Number of outgoing messages %d Total message length %" PetscInt_FMT "\n", nrqs, msz));
805: /* Determine the number of messages to expect, their lengths, from from-ids */
806: PetscCall(PetscGatherNumberOfMessages(comm, w2, w1, &nrqr));
807: PetscCall(PetscGatherMessageLengths(comm, nrqs, nrqr, w1, &onodes1, &olengths1));
809: /* Now post the Irecvs corresponding to these messages */
810: PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag0));
811: PetscCall(PetscPostIrecvInt(comm, tag0, nrqr, onodes1, olengths1, &rbuf1, &r_waits1));
813: /* Allocate Memory for outgoing messages */
814: PetscCall(PetscMalloc4(size, &sbuf1, size, &ptr, 2 * msz, &tmp, size, &ctr));
815: PetscCall(PetscArrayzero(sbuf1, size));
816: PetscCall(PetscArrayzero(ptr, size));
818: {
819: PetscInt *iptr = tmp;
820: PetscMPIInt k = 0;
821: for (PetscMPIInt i = 0; i < nrqs; i++) {
822: PetscMPIInt j = pa[i];
823: iptr += k;
824: sbuf1[j] = iptr;
825: k = w1[j];
826: }
827: }
829: /* Form the outgoing messages. Initialize the header space */
830: for (PetscMPIInt i = 0; i < nrqs; i++) {
831: PetscMPIInt j = pa[i];
833: sbuf1[j][0] = 0;
834: PetscCall(PetscArrayzero(sbuf1[j] + 1, 2 * w3[j]));
835: ptr[j] = sbuf1[j] + 2 * w3[j] + 1;
836: }
838: /* Parse the isrow and copy data into outbuf */
839: for (PetscInt i = 0; i < ismax; i++) {
840: row2proc_i = row2proc[i];
841: PetscCall(PetscArrayzero(ctr, size));
842: irow_i = irow[i];
843: jmax = nrow[i];
844: for (PetscInt j = 0; j < jmax; j++) { /* parse the indices of each IS */
845: proc = row2proc_i[j];
846: if (allrows[i]) row = j;
847: else row = irow_i[j];
849: if (proc != rank) { /* copy to the outgoing buf*/
850: ctr[proc]++;
851: *ptr[proc] = row;
852: ptr[proc]++;
853: }
854: }
855: /* Update the headers for the current IS */
856: for (PetscMPIInt j = 0; j < size; j++) { /* Can Optimise this loop too */
857: if ((ctr_j = ctr[j])) {
858: PetscInt k;
860: sbuf1_j = sbuf1[j];
861: k = ++sbuf1_j[0];
862: sbuf1_j[2 * k] = ctr_j;
863: sbuf1_j[2 * k - 1] = i;
864: }
865: }
866: }
868: /* Now post the sends */
869: PetscCall(PetscMalloc1(nrqs, &s_waits1));
870: for (PetscMPIInt i = 0; i < nrqs; ++i) {
871: PetscMPIInt j = pa[i];
872: PetscCallMPI(MPIU_Isend(sbuf1[j], w1[j], MPIU_INT, j, tag0, comm, s_waits1 + i));
873: }
875: /* Post Receives to capture the buffer size */
876: PetscCall(PetscMalloc1(nrqs, &r_waits2));
877: PetscCall(PetscMalloc3(nrqs, &req_source2, nrqs, &rbuf2, nrqs, &rbuf3));
878: if (nrqs) rbuf2[0] = tmp + msz;
879: for (PetscMPIInt i = 1; i < nrqs; ++i) rbuf2[i] = rbuf2[i - 1] + w1[pa[i - 1]];
880: for (PetscMPIInt i = 0; i < nrqs; ++i) {
881: PetscMPIInt j = pa[i];
882: PetscCallMPI(MPIU_Irecv(rbuf2[i], w1[j], MPIU_INT, j, tag2, comm, r_waits2 + i));
883: }
885: /* Send to other procs the buf size they should allocate */
886: /* Receive messages*/
887: PetscCall(PetscMalloc1(nrqr, &s_waits2));
888: PetscCall(PetscMalloc3(nrqr, &sbuf2, nrqr, &req_size, nrqr, &req_source1));
890: PetscCallMPI(MPI_Waitall(nrqr, r_waits1, MPI_STATUSES_IGNORE));
891: for (PetscMPIInt i = 0; i < nrqr; ++i) {
892: req_size[i] = 0;
893: rbuf1_i = rbuf1[i];
894: start = 2 * rbuf1_i[0] + 1;
895: end = olengths1[i];
896: PetscCall(PetscMalloc1(end, &sbuf2[i]));
897: sbuf2_i = sbuf2[i];
898: for (PetscInt j = start; j < end; j++) {
899: row = rbuf1_i[j] - rstart;
900: ncols = a_i[row + 1] - a_i[row] + b_i[row + 1] - b_i[row];
901: sbuf2_i[j] = ncols;
902: req_size[i] += ncols;
903: }
904: req_source1[i] = onodes1[i];
905: /* form the header */
906: sbuf2_i[0] = req_size[i];
907: for (PetscInt j = 1; j < start; j++) sbuf2_i[j] = rbuf1_i[j];
909: PetscCallMPI(MPIU_Isend(sbuf2_i, end, MPIU_INT, req_source1[i], tag2, comm, s_waits2 + i));
910: }
911: PetscCall(PetscFree(onodes1));
912: PetscCall(PetscFree(olengths1));
914: PetscCall(PetscFree(r_waits1));
915: PetscCall(PetscFree4(w1, w2, w3, w4));
917: /* Receive messages*/
918: PetscCall(PetscMalloc1(nrqs, &r_waits3));
920: PetscCallMPI(MPI_Waitall(nrqs, r_waits2, MPI_STATUSES_IGNORE));
921: for (PetscMPIInt i = 0; i < nrqs; ++i) {
922: PetscCall(PetscMalloc1(rbuf2[i][0], &rbuf3[i]));
923: req_source2[i] = pa[i];
924: PetscCallMPI(MPIU_Irecv(rbuf3[i], rbuf2[i][0], MPIU_INT, req_source2[i], tag3, comm, r_waits3 + i));
925: }
926: PetscCall(PetscFree(r_waits2));
928: /* Wait on sends1 and sends2 */
929: PetscCallMPI(MPI_Waitall(nrqs, s_waits1, MPI_STATUSES_IGNORE));
930: PetscCallMPI(MPI_Waitall(nrqr, s_waits2, MPI_STATUSES_IGNORE));
931: PetscCall(PetscFree(s_waits1));
932: PetscCall(PetscFree(s_waits2));
934: /* Now allocate sending buffers for a->j, and send them off */
935: PetscCall(PetscMalloc1(nrqr, &sbuf_aj));
936: jcnt = 0;
937: for (PetscMPIInt i = 0; i < nrqr; i++) jcnt += req_size[i];
938: if (nrqr) PetscCall(PetscMalloc1(jcnt, &sbuf_aj[0]));
939: for (PetscMPIInt i = 1; i < nrqr; i++) sbuf_aj[i] = sbuf_aj[i - 1] + req_size[i - 1];
941: PetscCall(PetscMalloc1(nrqr, &s_waits3));
942: {
943: for (PetscMPIInt i = 0; i < nrqr; i++) {
944: rbuf1_i = rbuf1[i];
945: sbuf_aj_i = sbuf_aj[i];
946: ct1 = 2 * rbuf1_i[0] + 1;
947: ct2 = 0;
948: for (PetscInt j = 1, max1 = rbuf1_i[0]; j <= max1; j++) {
949: kmax = rbuf1[i][2 * j];
950: for (PetscInt k = 0; k < kmax; k++, ct1++) {
951: PetscInt l;
952: row = rbuf1_i[ct1] - rstart;
953: nzA = a_i[row + 1] - a_i[row];
954: nzB = b_i[row + 1] - b_i[row];
955: ncols = nzA + nzB;
956: cworkA = PetscSafePointerPlusOffset(a_j, a_i[row]);
957: cworkB = PetscSafePointerPlusOffset(b_j, b_i[row]);
959: /* load the column indices for this row into cols */
960: cols = sbuf_aj_i + ct2;
961: for (l = 0; l < nzB; l++) {
962: if ((ctmp = bmap[cworkB[l]]) < cstart) cols[l] = ctmp;
963: else break;
964: }
965: imark = l;
966: for (l = 0; l < nzA; l++) cols[imark + l] = cstart + cworkA[l];
967: for (l = imark; l < nzB; l++) cols[nzA + l] = bmap[cworkB[l]];
968: ct2 += ncols;
969: }
970: }
971: PetscCallMPI(MPIU_Isend(sbuf_aj_i, req_size[i], MPIU_INT, req_source1[i], tag3, comm, s_waits3 + i));
972: }
973: }
975: /* create col map: global col of C -> local col of submatrices */
976: #if PetscDefined(USE_CTABLE)
977: for (PetscInt i = 0; i < ismax; i++) {
978: if (!allcolumns[i]) {
979: PetscCall(PetscHMapICreateWithSize(ncol[i], cmap + i));
981: jmax = ncol[i];
982: icol_i = icol[i];
983: cmap_i = cmap[i];
984: for (PetscInt j = 0; j < jmax; j++) PetscCall(PetscHMapISet(cmap[i], icol_i[j] + 1, j + 1));
985: } else cmap[i] = NULL;
986: }
987: #else
988: for (PetscInt i = 0; i < ismax; i++) {
989: if (!allcolumns[i]) {
990: PetscCall(PetscCalloc1(c->Nbs, &cmap[i]));
991: jmax = ncol[i];
992: icol_i = icol[i];
993: cmap_i = cmap[i];
994: for (PetscInt j = 0; j < jmax; j++) cmap_i[icol_i[j]] = j + 1;
995: } else cmap[i] = NULL;
996: }
997: #endif
999: /* Create lens which is required for MatCreate... */
1000: jcnt = 0;
1001: for (PetscInt i = 0; i < ismax; i++) jcnt += nrow[i];
1002: PetscCall(PetscMalloc1(ismax, &lens));
1003: if (ismax) PetscCall(PetscCalloc1(jcnt, &lens[0]));
1004: for (PetscInt i = 1; i < ismax; i++) lens[i] = PetscSafePointerPlusOffset(lens[i - 1], nrow[i - 1]);
1006: /* Update lens from local data */
1007: for (PetscInt i = 0; i < ismax; i++) {
1008: row2proc_i = row2proc[i];
1009: jmax = nrow[i];
1010: if (!allcolumns[i]) cmap_i = cmap[i];
1011: irow_i = irow[i];
1012: lens_i = lens[i];
1013: for (PetscInt j = 0; j < jmax; j++) {
1014: if (allrows[i]) row = j;
1015: else row = irow_i[j]; /* global blocked row of C */
1017: proc = row2proc_i[j];
1018: if (proc == rank) {
1019: /* Get indices from matA and then from matB */
1020: #if PetscDefined(USE_CTABLE)
1021: PetscInt tt;
1022: #endif
1023: row = row - rstart;
1024: nzA = a_i[row + 1] - a_i[row];
1025: nzB = b_i[row + 1] - b_i[row];
1026: cworkA = a_j + a_i[row];
1027: cworkB = PetscSafePointerPlusOffset(b_j, b_i[row]);
1029: if (!allcolumns[i]) {
1030: #if PetscDefined(USE_CTABLE)
1031: for (PetscInt k = 0; k < nzA; k++) {
1032: PetscCall(PetscHMapIGetWithDefault(cmap_i, cstart + cworkA[k] + 1, 0, &tt));
1033: if (tt) lens_i[j]++;
1034: }
1035: for (PetscInt k = 0; k < nzB; k++) {
1036: PetscCall(PetscHMapIGetWithDefault(cmap_i, bmap[cworkB[k]] + 1, 0, &tt));
1037: if (tt) lens_i[j]++;
1038: }
1040: #else
1041: for (PetscInt k = 0; k < nzA; k++) {
1042: if (cmap_i[cstart + cworkA[k]]) lens_i[j]++;
1043: }
1044: for (PetscInt k = 0; k < nzB; k++) {
1045: if (cmap_i[bmap[cworkB[k]]]) lens_i[j]++;
1046: }
1047: #endif
1048: } else { /* allcolumns */
1049: lens_i[j] = nzA + nzB;
1050: }
1051: }
1052: }
1053: }
1055: /* Create row map: global row of C -> local row of submatrices */
1056: for (PetscInt i = 0; i < ismax; i++) {
1057: if (!allrows[i]) {
1058: #if PetscDefined(USE_CTABLE)
1059: PetscCall(PetscHMapICreateWithSize(nrow[i], rmap + i));
1060: irow_i = irow[i];
1061: jmax = nrow[i];
1062: for (PetscInt j = 0; j < jmax; j++) {
1063: if (allrows[i]) {
1064: PetscCall(PetscHMapISet(rmap[i], j + 1, j + 1));
1065: } else {
1066: PetscCall(PetscHMapISet(rmap[i], irow_i[j] + 1, j + 1));
1067: }
1068: }
1069: #else
1070: PetscCall(PetscCalloc1(c->Mbs, &rmap[i]));
1071: rmap_i = rmap[i];
1072: irow_i = irow[i];
1073: jmax = nrow[i];
1074: for (PetscInt j = 0; j < jmax; j++) {
1075: if (allrows[i]) rmap_i[j] = j;
1076: else rmap_i[irow_i[j]] = j;
1077: }
1078: #endif
1079: } else rmap[i] = NULL;
1080: }
1082: /* Update lens from offproc data */
1083: {
1084: PetscInt *rbuf2_i, *rbuf3_i, *sbuf1_i;
1086: PetscCallMPI(MPI_Waitall(nrqs, r_waits3, MPI_STATUSES_IGNORE));
1087: for (tmp2 = 0; tmp2 < nrqs; tmp2++) {
1088: sbuf1_i = sbuf1[pa[tmp2]];
1089: jmax = sbuf1_i[0];
1090: ct1 = 2 * jmax + 1;
1091: ct2 = 0;
1092: rbuf2_i = rbuf2[tmp2];
1093: rbuf3_i = rbuf3[tmp2];
1094: for (PetscInt j = 1; j <= jmax; j++) {
1095: is_no = sbuf1_i[2 * j - 1];
1096: max1 = sbuf1_i[2 * j];
1097: lens_i = lens[is_no];
1098: if (!allcolumns[is_no]) cmap_i = cmap[is_no];
1099: rmap_i = rmap[is_no];
1100: for (PetscInt k = 0; k < max1; k++, ct1++) {
1101: if (allrows[is_no]) {
1102: row = sbuf1_i[ct1];
1103: } else {
1104: #if PetscDefined(USE_CTABLE)
1105: PetscCall(PetscHMapIGetWithDefault(rmap_i, sbuf1_i[ct1] + 1, 0, &row));
1106: row--;
1107: PetscCheck(row >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "row not found in table");
1108: #else
1109: row = rmap_i[sbuf1_i[ct1]]; /* the val in the new matrix to be */
1110: #endif
1111: }
1112: max2 = rbuf2_i[ct1];
1113: for (PetscInt l = 0; l < max2; l++, ct2++) {
1114: if (!allcolumns[is_no]) {
1115: #if PetscDefined(USE_CTABLE)
1116: PetscCall(PetscHMapIGetWithDefault(cmap_i, rbuf3_i[ct2] + 1, 0, &tcol));
1117: #else
1118: tcol = cmap_i[rbuf3_i[ct2]];
1119: #endif
1120: if (tcol) lens_i[row]++;
1121: } else { /* allcolumns */
1122: lens_i[row]++; /* lens_i[row] += max2 ? */
1123: }
1124: }
1125: }
1126: }
1127: }
1128: }
1129: PetscCall(PetscFree(r_waits3));
1130: PetscCallMPI(MPI_Waitall(nrqr, s_waits3, MPI_STATUSES_IGNORE));
1131: PetscCall(PetscFree(s_waits3));
1133: /* Create the submatrices */
1134: for (PetscInt i = 0; i < ismax; i++) {
1135: PetscCall(MatCreate(PETSC_COMM_SELF, submats + i));
1136: PetscCall(MatSetSizes(submats[i], nrow[i] * bs, ncol[i] * bs, PETSC_DETERMINE, PETSC_DETERMINE));
1138: PetscCall(MatSetType(submats[i], sym ? ((PetscObject)A)->type_name : MATSEQBAIJ));
1139: PetscCall(MatSeqBAIJSetPreallocation(submats[i], bs, 0, lens[i]));
1140: PetscCall(MatSeqSBAIJSetPreallocation(submats[i], bs, 0, lens[i])); /* this subroutine is used by SBAIJ routines */
1142: /* create struct Mat_SubSppt and attached it to submat */
1143: PetscCall(PetscNew(&smat_i));
1144: subc = (Mat_SeqBAIJ *)submats[i]->data;
1145: subc->submatis1 = smat_i;
1147: smat_i->destroy = submats[i]->ops->destroy;
1148: submats[i]->ops->destroy = MatDestroySubMatrix_SeqBAIJ;
1149: submats[i]->factortype = C->factortype;
1151: smat_i->id = i;
1152: smat_i->nrqs = nrqs;
1153: smat_i->nrqr = nrqr;
1154: smat_i->rbuf1 = rbuf1;
1155: smat_i->rbuf2 = rbuf2;
1156: smat_i->rbuf3 = rbuf3;
1157: smat_i->sbuf2 = sbuf2;
1158: smat_i->req_source2 = req_source2;
1160: smat_i->sbuf1 = sbuf1;
1161: smat_i->ptr = ptr;
1162: smat_i->tmp = tmp;
1163: smat_i->ctr = ctr;
1165: smat_i->pa = pa;
1166: smat_i->req_size = req_size;
1167: smat_i->req_source1 = req_source1;
1169: smat_i->allcolumns = allcolumns[i];
1170: smat_i->allrows = allrows[i];
1171: smat_i->singleis = PETSC_FALSE;
1172: smat_i->row2proc = row2proc[i];
1173: smat_i->rmap = rmap[i];
1174: smat_i->cmap = cmap[i];
1175: }
1177: if (!ismax) { /* Create dummy submats[0] for reuse struct subc */
1178: PetscCall(MatCreate(PETSC_COMM_SELF, &submats[0]));
1179: PetscCall(MatSetSizes(submats[0], 0, 0, PETSC_DETERMINE, PETSC_DETERMINE));
1180: PetscCall(MatSetType(submats[0], MATDUMMY));
1182: /* create struct Mat_SubSppt and attached it to submat */
1183: PetscCall(PetscNew(&smat_i));
1184: submats[0]->data = (void *)smat_i;
1186: smat_i->destroy = submats[0]->ops->destroy;
1187: submats[0]->ops->destroy = MatDestroySubMatrix_Dummy;
1188: submats[0]->factortype = C->factortype;
1190: smat_i->id = 0;
1191: smat_i->nrqs = nrqs;
1192: smat_i->nrqr = nrqr;
1193: smat_i->rbuf1 = rbuf1;
1194: smat_i->rbuf2 = rbuf2;
1195: smat_i->rbuf3 = rbuf3;
1196: smat_i->sbuf2 = sbuf2;
1197: smat_i->req_source2 = req_source2;
1199: smat_i->sbuf1 = sbuf1;
1200: smat_i->ptr = ptr;
1201: smat_i->tmp = tmp;
1202: smat_i->ctr = ctr;
1204: smat_i->pa = pa;
1205: smat_i->req_size = req_size;
1206: smat_i->req_source1 = req_source1;
1208: smat_i->allcolumns = PETSC_FALSE;
1209: smat_i->singleis = PETSC_FALSE;
1210: smat_i->row2proc = NULL;
1211: smat_i->rmap = NULL;
1212: smat_i->cmap = NULL;
1213: }
1215: if (ismax) PetscCall(PetscFree(lens[0]));
1216: PetscCall(PetscFree(lens));
1217: if (sbuf_aj) {
1218: PetscCall(PetscFree(sbuf_aj[0]));
1219: PetscCall(PetscFree(sbuf_aj));
1220: }
1222: } /* endof scall == MAT_INITIAL_MATRIX */
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++) {
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 = 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));
1270: ct2 += ncols;
1271: }
1272: }
1273: PetscCallMPI(MPIU_Isend(sbuf_aa_i, req_size[i] * bs2, MPIU_SCALAR, req_source1[i], tag4, comm, s_waits4 + i));
1274: }
1276: /* Assemble the matrices */
1277: /* First assemble the local rows */
1278: for (PetscInt i = 0; i < ismax; i++) {
1279: row2proc_i = row2proc[i];
1280: subc = (Mat_SeqBAIJ *)submats[i]->data;
1281: imat_ilen = subc->ilen;
1282: imat_j = subc->j;
1283: imat_i = subc->i;
1284: imat_a = subc->a;
1286: if (!allcolumns[i]) cmap_i = cmap[i];
1287: rmap_i = rmap[i];
1288: irow_i = irow[i];
1289: jmax = nrow[i];
1290: for (PetscInt j = 0; j < jmax; j++) {
1291: if (allrows[i]) row = j;
1292: else row = irow_i[j];
1293: proc = row2proc_i[j];
1295: if (proc == rank) {
1296: row = row - rstart;
1297: nzA = a_i[row + 1] - a_i[row];
1298: nzB = b_i[row + 1] - b_i[row];
1299: cworkA = a_j + a_i[row];
1300: cworkB = PetscSafePointerPlusOffset(b_j, b_i[row]);
1301: vworkA = a_a + a_i[row] * bs2;
1302: vworkB = PetscSafePointerPlusOffset(b_a, b_i[row] * bs2);
1304: if (allrows[i]) {
1305: row = row + rstart;
1306: } else {
1307: #if PetscDefined(USE_CTABLE)
1308: PetscCall(PetscHMapIGetWithDefault(rmap_i, row + rstart + 1, 0, &row));
1309: row--;
1311: PetscCheck(row >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "row not found in table");
1312: #else
1313: row = rmap_i[row + rstart];
1314: #endif
1315: }
1316: mat_i = imat_i[row];
1317: mat_a = PetscSafePointerPlusOffset(imat_a, mat_i * bs2);
1318: mat_j = PetscSafePointerPlusOffset(imat_j, mat_i);
1319: ilen = imat_ilen[row];
1321: /* load the column indices for this row into cols*/
1322: if (!allcolumns[i]) {
1323: PetscInt l;
1325: for (l = 0; l < nzB; l++) {
1326: if ((ctmp = bmap[cworkB[l]]) < cstart) {
1327: #if PetscDefined(USE_CTABLE)
1328: PetscCall(PetscHMapIGetWithDefault(cmap_i, ctmp + 1, 0, &tcol));
1329: if (tcol) {
1330: #else
1331: if ((tcol = cmap_i[ctmp])) {
1332: #endif
1333: *mat_j++ = tcol - 1;
1334: PetscCall(PetscArraycpy(mat_a, vworkB + l * bs2, bs2));
1335: mat_a += bs2;
1336: ilen++;
1337: }
1338: } else break;
1339: }
1340: imark = l;
1341: for (PetscInt l = 0; l < nzA; l++) {
1342: #if PetscDefined(USE_CTABLE)
1343: PetscCall(PetscHMapIGetWithDefault(cmap_i, cstart + cworkA[l] + 1, 0, &tcol));
1344: if (tcol) {
1345: #else
1346: if ((tcol = cmap_i[cstart + cworkA[l]])) {
1347: #endif
1348: *mat_j++ = tcol - 1;
1349: PetscCall(PetscArraycpy(mat_a, vworkA + l * bs2, bs2));
1350: mat_a += bs2;
1351: ilen++;
1352: }
1353: }
1354: for (l = imark; l < nzB; l++) {
1355: #if PetscDefined(USE_CTABLE)
1356: PetscCall(PetscHMapIGetWithDefault(cmap_i, bmap[cworkB[l]] + 1, 0, &tcol));
1357: if (tcol) {
1358: #else
1359: if ((tcol = cmap_i[bmap[cworkB[l]]])) {
1360: #endif
1361: *mat_j++ = tcol - 1;
1362: PetscCall(PetscArraycpy(mat_a, vworkB + l * bs2, bs2));
1363: mat_a += bs2;
1364: ilen++;
1365: }
1366: }
1367: } else { /* allcolumns */
1368: PetscInt l;
1369: for (l = 0; l < nzB; l++) {
1370: if ((ctmp = bmap[cworkB[l]]) < cstart) {
1371: *mat_j++ = ctmp;
1372: PetscCall(PetscArraycpy(mat_a, vworkB + l * bs2, bs2));
1373: mat_a += bs2;
1374: ilen++;
1375: } else break;
1376: }
1377: imark = l;
1378: for (l = 0; l < nzA; l++) {
1379: *mat_j++ = cstart + cworkA[l];
1380: PetscCall(PetscArraycpy(mat_a, vworkA + l * bs2, bs2));
1381: mat_a += bs2;
1382: ilen++;
1383: }
1384: for (l = imark; l < nzB; l++) {
1385: *mat_j++ = bmap[cworkB[l]];
1386: PetscCall(PetscArraycpy(mat_a, vworkB + l * bs2, bs2));
1387: mat_a += bs2;
1388: ilen++;
1389: }
1390: }
1391: imat_ilen[row] = ilen;
1392: }
1393: }
1394: }
1396: /* Now assemble the off proc rows */
1397: PetscCallMPI(MPI_Waitall(nrqs, r_waits4, MPI_STATUSES_IGNORE));
1398: for (tmp2 = 0; tmp2 < nrqs; tmp2++) {
1399: sbuf1_i = sbuf1[pa[tmp2]];
1400: jmax = sbuf1_i[0];
1401: ct1 = 2 * jmax + 1;
1402: ct2 = 0;
1403: rbuf2_i = rbuf2[tmp2];
1404: rbuf3_i = rbuf3[tmp2];
1405: rbuf4_i = rbuf4[tmp2];
1406: for (PetscInt j = 1; j <= jmax; j++) {
1407: is_no = sbuf1_i[2 * j - 1];
1408: rmap_i = rmap[is_no];
1409: if (!allcolumns[is_no]) cmap_i = cmap[is_no];
1410: subc = (Mat_SeqBAIJ *)submats[is_no]->data;
1411: imat_ilen = subc->ilen;
1412: imat_j = subc->j;
1413: imat_i = subc->i;
1414: imat_a = subc->a;
1415: max1 = sbuf1_i[2 * j];
1416: for (PetscInt k = 0; k < max1; k++, ct1++) { /* for each recved block row */
1417: row = sbuf1_i[ct1];
1419: if (allrows[is_no]) {
1420: row = sbuf1_i[ct1];
1421: } else {
1422: #if PetscDefined(USE_CTABLE)
1423: PetscCall(PetscHMapIGetWithDefault(rmap_i, row + 1, 0, &row));
1424: row--;
1425: PetscCheck(row >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "row not found in table");
1426: #else
1427: row = rmap_i[row];
1428: #endif
1429: }
1430: ilen = imat_ilen[row];
1431: mat_i = imat_i[row];
1432: mat_a = imat_a + mat_i * bs2;
1433: mat_j = imat_j + mat_i;
1434: max2 = rbuf2_i[ct1];
1435: if (!allcolumns[is_no]) {
1436: for (PetscInt l = 0; l < max2; l++, ct2++) {
1437: #if PetscDefined(USE_CTABLE)
1438: PetscCall(PetscHMapIGetWithDefault(cmap_i, rbuf3_i[ct2] + 1, 0, &tcol));
1439: #else
1440: tcol = cmap_i[rbuf3_i[ct2]];
1441: #endif
1442: if (tcol) {
1443: *mat_j++ = tcol - 1;
1444: PetscCall(PetscArraycpy(mat_a, rbuf4_i + ct2 * bs2, bs2));
1445: mat_a += bs2;
1446: ilen++;
1447: }
1448: }
1449: } else { /* allcolumns */
1450: for (PetscInt l = 0; l < max2; l++, ct2++, ilen++) {
1451: *mat_j++ = rbuf3_i[ct2]; /* same global column index of C */
1452: PetscCall(PetscArraycpy(mat_a, rbuf4_i + ct2 * bs2, bs2));
1453: mat_a += bs2;
1454: }
1455: }
1456: imat_ilen[row] = ilen;
1457: }
1458: }
1459: }
1461: if (!iscsorted) { /* sort column indices of the rows */
1462: MatScalar *work;
1464: PetscCall(PetscMalloc1(bs2, &work));
1465: for (PetscInt i = 0; i < ismax; i++) {
1466: subc = (Mat_SeqBAIJ *)submats[i]->data;
1467: imat_ilen = subc->ilen;
1468: imat_j = subc->j;
1469: imat_i = subc->i;
1470: imat_a = subc->a;
1471: if (allcolumns[i]) continue;
1473: jmax = nrow[i];
1474: for (PetscInt j = 0; j < jmax; j++) {
1475: mat_i = imat_i[j];
1476: mat_j = imat_j + mat_i;
1477: ilen = imat_ilen[j];
1478: mat_a = imat_a + mat_i * bs2;
1479: PetscCall(PetscSortIntWithDataArray(ilen, mat_j, mat_a, bs2 * sizeof(MatScalar), work));
1480: }
1481: }
1482: PetscCall(PetscFree(work));
1483: }
1485: PetscCall(PetscFree(r_waits4));
1486: PetscCallMPI(MPI_Waitall(nrqr, s_waits4, MPI_STATUSES_IGNORE));
1487: PetscCall(PetscFree(s_waits4));
1489: /* Restore the indices */
1490: for (PetscInt i = 0; i < ismax; i++) {
1491: if (!allrows[i]) PetscCall(ISRestoreIndices(isrow[i], irow + i));
1492: if (!allcolumns[i]) PetscCall(ISRestoreIndices(iscol[i], icol + i));
1493: }
1495: for (PetscInt i = 0; i < ismax; i++) {
1496: PetscCall(MatAssemblyBegin(submats[i], MAT_FINAL_ASSEMBLY));
1497: PetscCall(MatAssemblyEnd(submats[i], MAT_FINAL_ASSEMBLY));
1498: }
1500: PetscCall(PetscFree5(*(PetscInt ***)&irow, *(PetscInt ***)&icol, nrow, ncol, issorted));
1501: PetscCall(PetscFree5(row2proc, cmap, rmap, allcolumns, allrows));
1503: if (sbuf_aa) {
1504: PetscCall(PetscFree(sbuf_aa[0]));
1505: PetscCall(PetscFree(sbuf_aa));
1506: }
1508: for (PetscMPIInt i = 0; i < nrqs; ++i) PetscCall(PetscFree(rbuf4[i]));
1509: PetscCall(PetscFree(rbuf4));
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }