Actual source code: matij.c
petsc-3.3-p7 2013-05-11
1: #define PETSCMAT_DLL
2: #include <petsc-private/matimpl.h>
3: #include <../src/mat/impls/ij/stashij.h>
4: #include <../src/sys/utils/hash.h>
6: /*MC
7: MATIJ - MATIJ = "ij".
8: A matrix class encoding a PseudoGraph -- a directed graph that admits multiple edges between its vertices.
9: The underlying pseudograph, and therefore the matrix, can be interpreted as a multiset-valued or array-valued
10: map from vertices to vertices: each vertex v is mapped to the multiset or array of the vertices that terminate
11: the edges originating at v.
13: Vertices, edges, and local sizes:
14: Pseudograph vertices are identified with indices -- nonnegative integers of type PetscInt:
15: - domain indices, from which the edges emanate
16: - range or codomain indices, at which the edges terminate
17: Each processor owns the domain indices falling within the local ownership range (see MatGetOwnershipRange()).
19: Edges emanating from a local domain index correspond to the matrix entries in the corresponding local row.
20: Indices terminating the local edges can have any value in [0,N) (where N is Mat's global column size).
21: Since any global index can be the target of any local edge, or even of multiple local edges with the same
22: source index, the matrix column size does not reflect row sizes. In particular, the number of edges with the
23: same local source can be greater than N (where n is the global column size). As with MatMPIADJ, there is no
24: particular distinction attached to the local column size n.
27: Map, support, image(s):
28: The interpretation as an array-valued map allows MATIJ to define its action on indices or indexed arrays.
29: An array of indices with entries within the local ownership range can be mapped to the index array obtained by
30: a concatenation of the images of all of the input indices. Likewise, an indexed array of weights -- scalars,
31: integers or integer-scalar pairs -- can be mapped to a similar indexed array with the indices replaced by
32: their images, and the weights duplicated, if necessary.
34: Using the above map interpretation of MATIJ, the indices within the local ownership range and nonempty
35: images constitute the local support of the Mat -- an array of size m0 <= m. The indices that belong to any of
36: the images of the locally-supported indices constitute the local image of size n0 <= N.
38: Level: advanced
39: M*/
41: typedef struct {
42: PetscBool multivalued; /* Whether the underlying pseudograph is not a graph. */
43: /* The following data structures are using for stashing. */
44: MatStashMPIIJ stash;
45: /* The following data structures are used for mapping. */
46: PetscHashI hsupp; /* local support in a hash table */
47: PetscInt m,n;
48: PetscInt *ij; /* concatenated local images of ALL local input elements (i.e., all indices from the local ownership range), sorted within each image */
49: PetscInt *ijlen; /* image segment boundaries for each local input index */
50: PetscInt *image; /* local image */
51: PetscInt minijlen, maxijlen; /* min and max image segment size */
52: /* The following data structures are used for binning. */
53: PetscInt *binoffsets, *binsizes;
54: } Mat_IJ;
57: #define MatIJCheckAssembled(mat,needassembled,arg) \
58: do { \
59: if(!((mat)->assembled) && (needassembled)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, \
60: "MatIJ not assembled"); \
61: if(((mat)->assembled) && !(needassembled)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, \
62: "Mat already assembled"); \
63: } while(0)
67: #define MatIJGetSuppIndex_Private(A,mode,i,ii) \
68: do \
69: { \
70: (ii) = -1; \
71: if((mode) == MATIJ_LOCAL) { \
72: (ii) = i; \
73: } \
74: else { \
75: Mat_IJ *_13_ij = (Mat_IJ*)((A)->data); \
76: if(!_13_ij->hsupp) { \
77: if((ii) < (A)->rmap->rend && (ii) >= (i) - A->rmap->rstart) \
78: (ii) = (i) - (A)->rmap->rstart; \
79: else \
80: (ii) = -1; \
81: } \
82: else { \
83: PetscHashIMap(_13_ij->hsupp,(i),(ii)); \
84: } \
85: } \
86: } \
87: while(0)
89: #define MatIJGetIndexImage_Private(A,mode,i,ii) \
90: { \
91: if((mode) == MATIJ_LOCAL) { \
92: /* Assume image has been "localized". */ \
93: ii = i; \
94: } \
95: else { \
96: ii = !((Mat_IJ*)((A)->data))->image?i:(((Mat_IJ*)((A)->data))->image)[i]; \
97: } \
98: } \
103: static PetscErrorCode MatIJLocalizeImage_Private(Mat);
105: /*@C
106: MatIJMap - map an array of global indices (inidxi) with index (inidxj) and scalar (inval) weights, by pushing
107: the indices along the edges of the underlying pseudograph (see MATIJ).
108: Each locally-owned global index i from inidxi is replaced by the array of global indices terminating
109: the Mat's pseudograph edges that emanate from i, in the order the edges were provided to
110: MatIJSetEdges() or MatIJSetEdgesIS(); the individual image arrays are concatenated. inidxi ndices
111: outside the local ownership range or the local support are silently ignored -- replaced by
112: empty arrays. The weights from the domain indices are attached to the corresponding images, with
113: duplication, if necessary.
115: Not collective.
117: Input Parameters:
118: + A - pseudograph
119: . intype - (MATIJ_LOCAL | MATIJ_GLOBAL) meaning of inidxi: local support numbers or global indices
120: . insize - size of the input index and weight arrays; PETSC_NULL indicates _all_ support indices
121: . inidxi - array (of size insize) of global indices
122: . inidxj - array (of size insize) of index weights
123: . inval - array (of size insize) of scalar weights
124: - outtype - (MATIJ_LOCAL | MATIJ_GLOBAL) desired meaning of outdxi: local support numbers or global indices
126: Output Parameters:
127: + outsize - size of the output index and weight arrays
128: . outidxi - array (of size outsize) of the global indices adjacent to the indices in inidxi
129: . outidxj - array (of size outsize) of the index weights inherited by outidxi from inidxi
130: . outval - array (of size outsize) of the scalar weights inherited by outidxi from inidxi
131: - imgsizes - array (of size insize) of the sizes of image segments within outidxi for each i from inidxi
133: Level: advanced
134: .seealso: MatIJBin(), MatIJBinMap(), MatIJGetSupport()
135: @*/
138: PetscErrorCode MatIJMap(Mat A, MatIJIndexType intype, PetscInt insize, const PetscInt *inidxi, const PetscInt *inidxj, const PetscScalar *inval, MatIJIndexType outtype, PetscInt *outsize, PetscInt **outidxi, PetscInt **outidxj, PetscScalar **outval, PetscInt **outsizes)
139: {
141: Mat_IJ *pg = (Mat_IJ*)A->data;
142: PetscInt i,j,k,indi=0,indj,outsize_ = -1,*outidxi_ = PETSC_NULL, *outidxj_ = PETSC_NULL, *outsizes_ = PETSC_NULL;
143: PetscScalar *outval_ = PETSC_NULL;
146: if((outidxi && !*outidxi) || (outidxj && !*outidxj) || (outval && !*outval)) {
147: MatIJMap(A,intype,insize,inidxi,inidxj,inval,outtype,&outsize_,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
148: }
149: if(insize == PETSC_DETERMINE)
150: inidxi = PETSC_NULL;
151: else if(insize < 0)
152: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid input array size: %D", insize);
153: if(outidxi) {
154: if(!*outidxi) {
155: PetscMalloc(sizeof(PetscInt)*outsize_, outidxi);
156: }
157: outidxi_ = *outidxi;
158: }
159: if(outidxj) {
160: if(!*outidxj) {
161: PetscMalloc(sizeof(PetscInt)*outsize_, outidxj);
162: }
163: outidxj_= *outidxj;
164: }
165: if(outval) {
166: if(!*outval) {
167: PetscMalloc(sizeof(PetscScalar)*outsize_, outval);
168: }
169: outval_ = *outval;
170: }
171: if(outsizes_) {
172: if(!*outsizes) {
173: PetscMalloc(sizeof(PetscInt)*insize, outsizes);
174: }
175: outsizes_ = *outsizes;
176: }
178: if(intype == MATIJ_LOCAL && !pg->image) {
179: MatIJLocalizeImage_Private(A);
180: }
181: j = 0;
182: for(i = 0; i < insize; ++i) {
183: if(!inidxi) {
184: indi = i;
185: }
186: else {
187: /* Convert to local. */
188: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
189: if((indi < 0 || indi >= pg->m)){
190: /* drop */
191: if(outsizes_) outsizes_[i] = 0;
192: continue;
193: }
194: }
195: if(outidxi_ || (inval && outval_) || (inidxj && outidxj_) ) {
196: for(k = pg->ijlen[indi]; k < pg->ijlen[indi+1]; ++k) {
197: MatIJGetIndexImage_Private(A,outtype,pg->ij[k],indj);
198: if(outidxi_) outidxi_[j] = indj;
199: if(inidxj&&outidxj_) outidxj_[j] = inidxj[i];
200: if(inval&&outval_) outval_[j] = inval[i];
201: ++j;
202: }
203: }
204: else {
205: j += pg->ijlen[indi+1]-pg->ijlen[indi];
206: }
207: if(outsizes_) outsizes_[i] = (pg->ijlen[indi+1]-pg->ijlen[indi]);
208: }/* for(i = 0; i < len; ++i) */
209: if(outsize) *outsize = j;
210: return(0);
211: }
215: /*@C
216: MatIJBin - bin an array of global indices (inidxi) along with index (inidxj) and scalar (inval) weights by pushing the indices
217: along the edges of the underlying pseudograph (see MATIJ).
218: Each locally-owned global index i from inidxi is put in the arrays corresponding to the global indices
219: terminating the Mat's pseudograph edges that emanate from i. The bin arrays are ordered by the terminating
220: index. inidxi ndices outside the local ownership range or the local support are silently ignored --
221: contribute to no bins. The index weights in inidxj and inval are arranged into bins of their own, exactly mirroring
222: the binning of inidxi.
223:
225: Not collective.
227: Input Parameters:
228: + A - pseudograph
229: . intype - (MATIJ_LOCAL | MATIJ_GLOBAL) meaning of inidxi: local support numbers or global indices
230: . insize - size of the input index and weight arrays; PETSC_NULL indicates _all_ support indices
231: . inidxi - array (of size insize) of global indices
232: . inidxj - array (of size insize) of index weights
233: - inval - array (of size insize) of scalar weights
236: Output Parameters:
237: + outsize - size of the array of concatenated bins
238: . outidxi - array (of size outsize) containing the binned indices from inidxi
239: . outidxj - array (of size outsize) containing the binned index weights from inidxj
240: . outval - array (of size outsize) containing the binned scalar weights from inval
241: - binsizes - array (of size n) of bin sizes
243: Note: n0 is the local image size -- the number of indices terminating the locally-supported indices
244: (see MATIJ) -- and can be obtained with MatIJGetImageSize().
246: Level: advanced
247: .seealso: MatIJMap(), MatIJBinMap(), MatIJGetSupport(), MatIJGetImageSize()
248: @*/
251: PetscErrorCode MatIJBin(Mat A, MatIJIndexType intype, PetscInt insize, const PetscInt *inidxi, const PetscInt *inidxj, const PetscScalar *inval, PetscInt *outsize, PetscInt **outidxi, PetscInt **outidxj, PetscScalar **outval, PetscInt **binsizes)
252: {
254: Mat_IJ *pg = (Mat_IJ*)A->data;
255: PetscInt i,j,k,indi=0,outsize_ = -1,*outidxi_ = PETSC_NULL, *outidxj_ = PETSC_NULL, *binsizes_ = PETSC_NULL;
256: PetscScalar *outval_ = PETSC_NULL;
259: /* Binning requires a localized image. */
260: MatIJLocalizeImage_Private(A);
261: if((outidxi && !*outidxi) || (outidxj && !*outidxj) || (outval && !*outval)) {
262: MatIJBin(A,intype,insize,inidxi,PETSC_NULL,PETSC_NULL,&outsize_,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
263: }
264: if(insize == PETSC_DETERMINE){
265: insize = pg->m;
266: inidxi = PETSC_NULL;
267: }
268: else if(insize < 0)
269: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid input array size: %D", insize);
270: if(outidxi) {
271: if(!*outidxi) {
272: PetscMalloc(sizeof(PetscInt)*outsize_,outidxi);
273: }
274: outidxi_ = *outidxi;
275: }
276: if(outidxj) {
277: if(!*outidxj) {
278: PetscMalloc(sizeof(PetscInt)*outsize_,outidxj);
279: }
280: outidxj_ = *outidxj;
281: }
282: if(outval) {
283: if(!*outval) {
284: PetscMalloc(sizeof(PetscScalar)*outsize_,outval);
285: }
286: outval_ = *outval;
287: }
288: if(binsizes) {
289: if(!*binsizes) {
290: PetscMalloc(sizeof(PetscInt)*(pg->n), binsizes);
291: }
292: binsizes_ = *binsizes;
293: }
295: /* We'll need to count the contributions to each "bin" and the offset of each bin in outidx. */
296: /* Allocate the bin offset array, if necessary. */
297: if(!pg->binoffsets) {
298: PetscMalloc((pg->n+1)*sizeof(PetscInt), &(pg->binoffsets));
299: }
300: /* Initialize bin offsets */
301: for(j = 0; j <= pg->n; ++j) {
302: pg->binoffsets[j] = 0;
303: }
304: for(i = 0; i < insize; ++i) {
305: if(!inidxi) {
306: indi = i;
307: }
308: else {
309: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
310: if((indi < 0 || indi >=pg->m)){
311: /* drop */
312: continue;
313: }
314: }
315: for(k = pg->ijlen[indi]; k < pg->ijlen[indi+1]; ++k) {
316: ++(pg->binoffsets[pg->ij[k]+1]);
317: }
318: }/* for(i = 0; i < insize; ++i) */
319: /* Convert bin sizes into bin offsets. */
320: for(j = 0; j < pg->n; ++j) {
321: pg->binoffsets[j+1] += pg->binoffsets[j];
322: }
323: /* Now bin the input indices and values. */
324: if(outidxi_ || (inval && outval) || (inidxj && outidxj_) ) {
325: /* Allocate the bin size array, if necessary. */
326: if(!binsizes_) {
327: if(!pg->binsizes) {
328: PetscMalloc((pg->n)*sizeof(PetscInt), &(pg->binsizes));
329: }
330: binsizes_ = pg->binsizes;
331: }
332: /* Initialize bin sizes to zero. */
333: for(j = 0; j < pg->n; ++j) {
334: binsizes_[j] = 0;
335: }
336: for(i = 0; i < insize; ++i) {
337: if(!inidxi) {
338: indi = i;
339: }
340: else {
341: /* Convert to local. */
342: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
343: }
344: if((indi < 0 || indi >= pg->m)){
345: /* drop */
346: continue;
347: }
348: for(k = pg->ijlen[indi]; k < pg->ijlen[indi+1]; ++k) {
349: j = pg->ij[k];
350: if(outidxi_) outidxi_[pg->binoffsets[j]+binsizes_[j]] = inidxi[i];
351: if(outval_ && inval) outval_ [pg->binoffsets[j]+binsizes_[j]] = inval[i];
352: if(outidxj_ && inidxj) outidxj_[pg->binoffsets[j]+binsizes_[j]] = inidxj[i];
353: ++binsizes[j];
354: }
355: }/* for(i = 0; i < insize; ++i) */
356: }/* if(outidxi_ || (inval && outval_) || (inidxj && outidxj_) ) */
357: if(outsize) *outsize = pg->binoffsets[pg->n];
358: return(0);
359: }
362: /*@C
363: MatIJBinMap - simultaneously bin and map an array of indices (inidxi) along with index (inidxj) and scalar (inval) weights
364: by pushing the indices along the edges of two pseudographs (see MATIJ, MatIJMap(), MatIJBin()).
365: Each locally-supported index i from inidxi is assigned to the arrays (bins) corresponding to the global
366: indices terminating the A's pseudograph edges that emanate from i. i's location in each bin is occupied
367: by the index terminating the corresponding pseudograph edge that emanate from i in B. Thus, A and B must be
368: compatible in the following sense: they must have the same local suppors and local image sizes.
369: inidxi indices outside the local support are silently ignored -- contribute to no bins. The index (inidxj)
370: and scalar (inval) weights are arranged in bins of their own, exactly mirroring the binning of inidxi.
372: Not collective.
374: Input Parameters:
375: + A - binning pseudograph
376: . B - mapping pseudograph
377: . intype - (MATIJ_LOCAL | MATIJ_GLOBAL) meaning of inidxi: local support numbers or global indices
378: . insize - size of the input index and weight arrays
379: . inidxi - array (of size insize) of global indices
380: . inidxj - array (of size insize) of index weights
381: . inval - array (of size insize) of scalar weights
382: - outtype - (MATIJ_LOCAL | MATIJ_GLOBAL) desired meaning of inidxi: local image numbers or global indices
385: Output Parameters:
386: + outsize - size of the array of concatenated bins
387: . outidxi - array (of size outsize) containing the binned images of the indices from inidxi
388: . outidxj - array (of size outsize) containing the binned index weights from inidxj
389: . outval - array (of size outsize) containing the binned scalar weights from inval
390: - binsizes - array (of size n0) of bin sizes
392: Note:
393: + The idea behind MatIJBinMap is that the binning is done by A, while what is actually binned comes from B.
394: Pseudographs A and B are structurally isomorphic. Moreover, they only differ in the terminating indices:
395: edges i-e->jA and i-eB->jB are in a one-to-one correspondence, if their positions in the ordering of i's
396: images are the same. Each source index i is assigned to every bin jA labeled by each of the indices
397: attached to i in A by some eA, but the binned value is jB -- the index attached to i in B by eB, which
398: corresponds to eA.
399: - Another way of viewing the pseudograph pair A and B is as a single pseudograph (B), whose edges are
400: colored (by A's terminating indices), and ordered on the color within each originating index's image.
403: Level: advanced
404: .seealso: MATIJ, MatIJBin(), MatIJMap(), MatIJGetSupport(), MatIJGetImage(), MatIJGetRowSizes()
405: @*/
408: PetscErrorCode MatIJBinMap(Mat A, Mat B, MatIJIndexType intype, PetscInt insize, const PetscInt *inidxi, const PetscInt *inidxj, const PetscScalar *inval, MatIJIndexType outtype, PetscInt *outsize, PetscInt **outidxi, PetscInt **outidxj, PetscScalar **outval, PetscInt **binsizes)
409: {
411: Mat_IJ *pga = (Mat_IJ*)A->data;
412: Mat_IJ *pgb = (Mat_IJ*)B->data;
413: PetscBool isij;
414: PetscInt indi = -1, indj, i,j,k,outsize_ = -1,*outidxi_ = PETSC_NULL, *outidxj_ = PETSC_NULL, *binsizes_ = PETSC_NULL;
416: PetscScalar *outval_ = PETSC_NULL;
420: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
421: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix 1 not of type MATIJ: %s", ((PetscObject)A)->type);
423: PetscObjectTypeCompare((PetscObject)B,MATIJ,&isij);
424: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix 2 not of type MATIJ: %s", ((PetscObject)B)->type);
427: if(A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible local row sizes: %D and %D", A->rmap->n, B->rmap->n);
428: if(A->rmap->N != B->rmap->N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible global row sizes: %D and %D", A->rmap->N, B->rmap->N);
429: if(A->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible local column sizes: %D and %D", A->cmap->n, B->cmap->n);
430: if(A->cmap->N != B->cmap->N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible global column sizes: %D and %D", A->cmap->N, B->cmap->N);
432: if(pga->m != pgb->m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible local support sizes: %D and %D", pga->m, pgb->m);
433: if(pga->n != pgb->n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible local image sizes: %D and %D", pga->n, pgb->n);
435: /* Binning requires a localized image. */
436: MatIJLocalizeImage_Private(A);
437: if((outidxi && !*outidxi) || (outidxj && !*outidxj) || (outval && !*outval)) {
438: MatIJBinMap(A,B,intype,insize,inidxi,inidxj,inval,outtype,&outsize_,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
439: }
440: if(insize == PETSC_DETERMINE){
441: insize = pga->m;
442: inidxi = PETSC_NULL;
443: }
444: else if(insize < 0)
445: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid input array size: %D", insize);
446: if(outidxi) {
447: if(!*outidxi) {
448: PetscMalloc(sizeof(PetscInt)*outsize_, outidxi);
449: }
450: outidxi_ = *outidxi;
451: }
452: if(outidxj) {
453: if(!*outidxj) {
454: PetscMalloc(sizeof(PetscInt)*outsize_, outidxj);
455: }
456: outidxj_ = *outidxj;
457: }
458: if(outval) {
459: if(!*outval) {
460: PetscMalloc(sizeof(PetscScalar)*outsize_, outval);
461: }
462: outval_ = *outval;
463: }
464: if(binsizes) {
465: if(!*binsizes) {
466: PetscMalloc(sizeof(PetscInt)*(pga->n), binsizes);
467: }
468: binsizes_ = *binsizes;
469: }
471: /* We'll need to count the contributions to each "bin" and the offset of each bin in outidx_. */
472: /* Allocate the bin offset array, if necessary. */
473: if(!pga->binoffsets) {
474: PetscMalloc((pga->n+1)*sizeof(PetscInt), &(pga->binoffsets));
475: }
476: /* Allocate the bin size array, if necessary. */
477: if(!binsizes_) {
478: if(!pga->binsizes) {
479: PetscMalloc((pga->n)*sizeof(PetscInt), &(pga->binsizes));
480: }
481: binsizes_ = pga->binsizes;
482: }
483: /* Initialize bin offsets */
484: for(j = 0; j <= pga->n; ++j) {
485: pga->binoffsets[j] = 0;
486: }
487: for(i = 0; i < insize; ++i) {
488: if(!inidxi) {
489: indi = i;
490: }
491: else {
492: /* Convert to local. */
493: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
494: if((indi < 0 || indi >= pga->m)){
495: /* drop */
496: continue;
497: }
498: }
499: if(pga->ijlen[indi] != pgb->ijlen[indi])
500: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Image sizes different for local index %D = indi: %D and %D", indi, pga->ijlen[indi], pgb->ijlen[indi]);
501: for(k = pga->ijlen[indi]; k < pga->ijlen[indi+1]; ++k) {
502: ++(pga->binoffsets[pga->ij[k]+1]);
503: }
504: }/* for(i = 0; i < insize; ++i) */
505: /* Convert bin sizes into bin offsets. */
506: for(j = 0; j < pga->n; ++j) {
507: pga->binoffsets[j+1] += pga->binoffsets[j];
508: }
509: /* Now bin the input indices and values. */
510: if(outidxi_ || (inval && outval_) || (inidxj && outidxj_) ) {
511: /* Initialize bin sizes to zero. */
512: for(j = 0; j < pga->n; ++j) {
513: binsizes_[j] = 0;
514: }
515: for(i = 0; i < insize; ++i) {
516: if(!inidxi) {
517: indi = inidxi[i];
518: }
519: else {
520: /* Convert to local. */
521: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
522: if((indi < 0 || indi >= pga->m)){
523: /* drop */
524: continue;
525: }
526: }
527: for(k = pga->ijlen[indi]; k < pga->ijlen[indi+1]; ++k) {
528: j = pga->ij[k];
529: MatIJGetIndexImage_Private(A,outtype,pgb->ij[k],indj);
530: if(outidxi_) outidxi_[pga->binoffsets[j]+binsizes_[j]] = indj;
531: if(outval_ && inval) outval_[pga->binoffsets[j] +binsizes_[j]] = inval[i];
532: if(outidxj_ && inidxj) outidxj_[pga->binoffsets[j]+binsizes_[j]] = inidxj[i];
533: ++binsizes_[j];
534: }
535: }/* for(i = 0; i < insize; ++i) */
536: }/* if(outidxi_ || (inval && outval_) || (inidxj && outidxj_) ) */
537: return(0);
538: }
542: PetscErrorCode MatGetRow_IJ(Mat A, PetscInt row, PetscInt *rowsize, PetscInt *cols[], PetscScalar *vals[]) {
543: PetscInt off,len,i,r;
545: Mat_IJ *pg = (Mat_IJ*)A->data;
548: /* It is easy to implement this, but will only be done, if there is demand. */
549: if(rowsize) *rowsize = 0;
550: if(cols) *cols = PETSC_NULL;
551: if(vals) *vals = PETSC_NULL;
553: /* Convert to local. */
554: MatIJGetSuppIndex_Private(A,MATIJ_GLOBAL,row,r);
555: if((r >= 0 && r < pg->m)){
556: off = pg->ijlen[r];
557: len = pg->ijlen[r+1]-pg->ijlen[r];
558: if(cols) *cols = pg->ij+off;
559: if(rowsize) *rowsize = len;
560: if(vals) {
561: PetscMalloc(sizeof(PetscScalar)*len, vals);
562: for(i = 0; i < len; ++i) {
563: (*vals)[i] = (PetscScalar)1.0;
564: }
565: }
566: }
567: return(0);
568: }
572: PetscErrorCode MatRestoreRow_IJ(Mat A, PetscInt row, PetscInt *rowsize, PetscInt *cols[], PetscScalar *vals[]) {
574: PetscInt r;
576: Mat_IJ *pg = (Mat_IJ*)A->data;
580: /* Convert to local. */
581: MatIJGetSuppIndex_Private(A,MATIJ_GLOBAL,row,r);
582: if((r >= 0 && r < pg->m)){
583: if(vals) {
584: PetscFree(*vals);
585: }
586: }
587: return(0);
588: }
591: /*@C
592: MatIJSetMultivalued - indicates whether the underlying pseudograph is a multivalued or not (a graph).
593:
594: Not collective.
596: Input arguments:
597: + A - pseudograph
598: - multivalued - whether the matrix encodes a multivalued (pseudo)graph.
600: Level: advanced
602: .seealso: MatIJGetMultivalued(), MatIJSetEdges(), MatIJGetEdges(), MatIJGetSupport(), MatIJGetImage()
603: @*/
604: #undef __FUNCT__
606: PetscErrorCode MatIJSetMultivalued(Mat A, PetscBool multivalued)
607: {
609: Mat_IJ *pg = (Mat_IJ *)(A->data);
610: PetscBool isij;
613: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
614: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
615: MatIJCheckAssembled(A,PETSC_FALSE,1);
616: MatStashMPIIJSetMultivalued_Private(pg->stash,multivalued);
617: pg->multivalued = multivalued;
618: return(0);
619: }
621: /*@C
622: MatIJGetMultivalued - return a flag indicating whether the underlying pseudograph is a multivalued or not (a graph).
623:
624: Not collective.
626: Input arguments:
627: . A - pseudograph
629: Output arguments:
630: . multivalued - whether the matrix encodes a multivalued (pseudo)graph.
632: Level: advanced
634: .seealso: MatIJSetMultivalued(), MatIJSetEdges(), MatIJGetEdges(), MatIJGetSupport(), MatIJGetImage()
635: @*/
636: #undef __FUNCT__
638: PetscErrorCode MatIJGetMultivalued(Mat A, PetscBool *multivalued)
639: {
641: Mat_IJ *pg = (Mat_IJ *)(A->data);
642: PetscBool isij;
645: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
646: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
648: *multivalued = pg->multivalued;
649: return(0);
650: }
653: /* Clears everything, but the stash and multivalued. */
656: static PetscErrorCode MatIJClear_Private(Mat mat)
657: {
658: Mat_IJ *pg = (Mat_IJ*)(mat->data);
661: if(pg->hsupp) {
662: PetscHashIDestroy((pg->hsupp));
663: }
664: if(pg->image) {
665: PetscFree(pg->image);
666: }
667: if(pg->ij) {
668: PetscFree(pg->ij);
669: }
670: if(pg->ijlen) {
671: PetscFree(pg->ijlen);
672: }
673: if(pg->binoffsets) {
674: PetscFree(pg->binoffsets);
675: }
676: pg->m = pg->minijlen = pg->maxijlen = 0;
677: pg->n = PETSC_DETERMINE;
678: return(0);
679: }
681: /*@C
682: MatIJSetEdgesIS - sets the edges in the pseudograph matrix.
683: The edges are specified as two index sets of vertices of equal length:
684: outgoing and incoming vertices (ix -> iy).
685:
686: Not collective
688: Input parameters:
689: + A - pseudograph
690: . ix - list of outgoing vertices
691: - iy - list of incoming vertices
693: Note:
694: + This will cause the matrix to be be put in an unassembled state.
695: . Edges are assembled during MatAssembly -- moved to the processor owning the outgoing vertex.
696: - Communicators of the IS objects must match that of MatIJ.
698: Level: intermediate
700: .seealso: MATIJ, MatIJSetEdges(), MatIJGetEdgesIS(), MatIJGetSupportIS(), MatIJGetImageIS()
701: @*/
702: #undef __FUNCT__
704: PetscErrorCode MatIJSetEdgesIS(Mat A, IS ix, IS iy)
705: {
706: IS iix, iiy;
707: PetscInt nix, niy;
708: const PetscInt *ixidx, *iyidx;
710: PetscBool isij;
714: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
715: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
723: if(!ix) {
724: ISCreateStride(((PetscObject)A)->comm, A->rmap->n, A->rmap->rstart, 1, &(iix));
725: nix = A->rmap->n;
726: }
727: else
728: iix = ix;
729: if(!iy) {
730: ISCreateStride(((PetscObject)A)->comm, A->cmap->n, A->cmap->rstart, 1, &(iiy));
731: niy = A->cmap->n;
732: }
733: else
734: iiy = iy;
735: ISGetLocalSize(iix,&nix);
736: ISGetLocalSize(iiy,&niy);
737: if(nix != niy) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %D and %D", nix, niy);
739: ISGetIndices(iix, &ixidx);
740: ISGetIndices(iiy, &iyidx);
741: MatIJSetEdges(A,nix,ixidx, iyidx);
742: ISRestoreIndices(iix, &ixidx);
743: ISRestoreIndices(iiy, &iyidx);
744: if(!ix) {
745: ISDestroy(&iix);
746: }
747: if(!iy) {
748: ISDestroy(&iiy);
749: }
750: return(0);
751: }
753: /*@C
754: MatIJSetEdges - sets the edges in the pseudograph matrix.
755: The edges are specified as two integer arrays of vertices of equal length:
756: outgoing and incoming vertices (ix -> iy).
757:
758: Not collective
760: Input parameters:
761: + A - pseudograph
762: . len - length of vertex arrays
763: . ixidx - list of outgoing vertices
764: - iyidx - list of incoming vertices
766: Note:
767: + This will cause the matrix to be be put in an unassembled state.
768: - Edges are assembled during MatAssembly -- moved to the processor owning the outgoing vertex.
770: Level: intermediate
772: .seealso: MATIJ, MatIJSetEdgesIS(), MatIJGetEdges(), MatIJGetSupport(), MatIJGetImage()
773: @*/
774: #undef __FUNCT__
776: PetscErrorCode MatIJSetEdges(Mat A, PetscInt len, const PetscInt *ixidx, const PetscInt *iyidx)
777: {
778: Mat_IJ *pg = (Mat_IJ*)(A->data);
779: PetscInt *iixidx = PETSC_NULL, *iiyidx = PETSC_NULL, k;
781: PetscBool isij;
785: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
786: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
788: if(len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Negative edge array length: %D", len);
789:
790: if(!ixidx){
791: if(len != A->rmap->n)
792: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The length of an empty source array %D must equal the local row size %D", len, A->rmap->n);
793: PetscMalloc(len*sizeof(PetscInt), &iixidx);
794: for(k = 0; k < len; ++k) {
795: iixidx[k] = A->rmap->rstart + k;
796: }
797: ixidx = iixidx;
798: }
799: if(!iyidx) {
800: if(len != A->cmap->n)
801: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The length of an empty target array %D must equal the local column size %D", len, A->cmap->n);
802: for(k = 0; k < len; ++k) {
803: iiyidx[k] = A->cmap->rstart + k;
804: }
805: iyidx = iiyidx;
806: }
807: MatStashMPIIJExtend_Private(pg->stash, len, ixidx, iyidx);
808: if(!iixidx) {
809: PetscFree(iixidx);
810: }
811: if(!iiyidx) {
812: PetscFree(iiyidx);
813: }
814: A->was_assembled = A->assembled;
815: A->assembled = PETSC_FALSE;
816: return(0);
817: }
819: #undef __FUNCT__
821: static PetscErrorCode MatIJGetAssembledEdges_Private(Mat A, PetscInt *len, PetscInt **ixidx, PetscInt **iyidx)
822: {
824: Mat_IJ *pg = (Mat_IJ *)(A->data);
825: PetscInt len_, *ixidx_ = PETSC_NULL,*iyidx_ = PETSC_NULL, k,i, ii;
826: PetscHashIIter hi;
827: PetscBool isij;
831: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
832: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
834: if(!len && !ixidx && !iyidx) return(0);
836: len_ = pg->ijlen[pg->m];
837: if(len) *len = len_;
838: if(!ixidx && !iyidx) return(0);
839: if(ixidx) {
840: if(!*ixidx) {
841: PetscMalloc(sizeof(PetscInt)*len_, ixidx);
842: }
843: ixidx_ = *ixidx;
844: }
845: if(iyidx) {
846: if(!*iyidx) {
847: PetscMalloc(sizeof(PetscInt)*len_, iyidx);
848: }
849: iyidx_ = *iyidx;
850: }
851: if(pg->hsupp) {
852: PetscHashIIterBegin(pg->hsupp,hi);
853: while(!PetscHashIIterAtEnd(pg->hsupp,hi)){
854: PetscHashIIterGetKeyVal(pg->hsupp,hi,ii,i);
855: for(k = pg->ijlen[i]; k < pg->ijlen[i+1]; ++k) {
856: if(ixidx_) ixidx_[k] = ii;
857: if(iyidx_) iyidx_[k] = pg->image[pg->ij[k]];
858: }
859: PetscHashIIterNext(pg->hsupp,hi);
860: }
861: }
862: else {
863: for(i = 0; i < pg->m; ++i) {
864: for(k = pg->ijlen[i]; k < pg->ijlen[i+1]; ++k) {
865: if(ixidx_) ixidx_[k] = i + A->rmap->rstart;
866: if(iyidx_) iyidx_[k] = pg->image[pg->ij[k]];
867: }
868: }
869: }
870: return(0);
871: }
874: /*@C
875: MatIJGetEdges - retrieves the edges from a pseudograph matrix.
876: The edges are specified as two integer arrays of vertices of equal length:
877: outgoing and incoming vertices (ix -> iy).
878:
879: Not collective
881: Input parameters:
882: . A - pseudograph
884: Output parameters:
885: + len - length of vertex arrays
886: . ixidx - list of outgoing vertices
887: - iyidx - list of incoming vertices
890: Notes:
891: + Both assembled and unassembled edges are returned.
892: - For an assembled matrix the retrieved outgoing vertices are guaranteed to be locally-owned.
894: Level: advanced
896: .seealso: MATIJ, MatIJSetEdges(), MatIJGetEdgesIS(), MatIJGetSupport(), MatIJGetImage()
897: @*/
898: #undef __FUNCT__
900: PetscErrorCode MatIJGetEdges(Mat A, PetscInt *len, PetscInt **ixidx, PetscInt **iyidx)
901: {
903: Mat_IJ *pg = (Mat_IJ *)(A->data);
904: PetscInt len_, lenI, lenII;
905: PetscInt *ixidxI = PETSC_NULL, *iyidxI = PETSC_NULL, *ixidxII = PETSC_NULL, *iyidxII = PETSC_NULL;
906: PetscBool isij;
910: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
911: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
913: if(!len && !ixidx && !iyidx) return(0);
915: MatIJGetAssembledEdges_Private(A, &lenI, PETSC_NULL, PETSC_NULL);
916: MatStashMPIIJGetIndicesMerged_Private(pg->stash, &lenII, PETSC_NULL, PETSC_NULL);
917:
918: len_ = lenI + lenII;
919: if(len) *len = len_;
921: if(!len_ || (!ixidx && !iyidx)) return(0);
923: if(!ixidx || !iyidx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Vertex array pointers must be null or non-null together");
925: if(!lenI) {
926: /* Only stash indices need to be returned. */
927: MatStashMPIIJGetIndicesMerged_Private(pg->stash, len, ixidx, iyidx);
928: }
929: else if(!lenII) {
930: /* Only assembled edges must be returned. */
931: MatIJGetAssembledEdges_Private(A, len, ixidx, iyidx);
932: }
933: else {
934: /* Retrieve the two sets of indices. */
935: MatIJGetAssembledEdges_Private(A, &lenI, &ixidxI, &iyidxI);
936: MatStashMPIIJGetIndicesMerged_Private(pg->stash, &lenII, &ixidxII, &iyidxII);
937: /* Merge. */
938: PetscMergeIntArrayPair(lenI,ixidxI,iyidxI,lenII,ixidxII,iyidxII,len,ixidx,iyidx);
939: /* Clean up. */
940: PetscFree(ixidxI);
941: PetscFree(iyidxI);
942: PetscFree(ixidxII);
943: PetscFree(iyidxII);
944: }
945: return(0);
946: }
948: /*@C
949: MatIJGetEdgesIS - retrieves the edges in the boolean matrix graph.
950: The edges are specified as two index sets of vertices of equal length:
951: outgoing and incoming vertices (ix -> iy).
952:
953: Not collective
955: Input parameters:
956: . A - pseudograph
958: Output parameters:
959: + ix - IS of outgoing vertices
960: - iy - IS of incoming vertices
962: Note:
963: + Both assembled and unassembled edges are returned.
964: . For an assembled matrix the retrieved outgoing vertices are guaranteed to be locally-owned.
965: - ix and iy will have the same communicator as MatIJ and will have the same length.
968: Level: intermediate
970: .seealso: MATIJ, MatIJSetEdgesIS(), MatIJGetEdges(), MatIJGetSupportIS(), MatIJGetImageIS()
971: @*/
972: #undef __FUNCT__
974: PetscErrorCode MatIJGetEdgesIS(Mat A, IS *ix, IS *iy)
975: {
977: PetscInt len, *ixidx = PETSC_NULL, *iyidx = PETSC_NULL, **_ixidx = PETSC_NULL, **_iyidx = PETSC_NULL;
978: PetscBool isij;
982: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
983: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
984: if(ix){
985: _ixidx = &ixidx;
986: }
987: if(iy) {
988: _iyidx = &iyidx;
989: }
990: MatIJGetEdges(A, &len, _ixidx, _iyidx);
991: if(ix) {
992: ISCreateGeneral(A->rmap->comm, len, ixidx, PETSC_OWN_POINTER, ix);
993: }
994: if(iy) {
995: ISCreateGeneral(A->cmap->comm, len, iyidx, PETSC_OWN_POINTER, iy);
996: }
997: return(0);
998: }
1000: /*
1001: Sort iy and store the unique in a (temporary) hash table to determine the local image.
1002: Endow the global image indices with a local number, then replace global indices in ij
1003: with the local numbers. Store the global image in an array: each local number will
1004: naturally serve as the index into the array for the corresponding global index.
1005: */
1008: static PetscErrorCode MatIJLocalizeImage_Private(Mat A)
1009: {
1011: Mat_IJ *pg = (Mat_IJ *) A->data;
1012: PetscInt i,j,k,n,totalnij,*image;
1013: PetscHashI himage;
1015: if(pg->image) return(0);
1017: /* (A) Construct image -- the set of unique ij indices. */
1018: /* (B) Endow the image with a local numbering. */
1019: /* (C) Then convert ij to this local numbering. */
1021: /* (A) */
1022: /* Insert all of the ij into himage. */
1023: PetscHashICreate(himage);
1024: for(i = 0; i < pg->m; ++i) {
1025: for(k = pg->ijlen[i]; pg->ijlen[i+1]; ++i) {
1026: PetscHashIAdd(himage,pg->ij[k],0);
1027: }
1028: }
1029: /* (B) */
1030: /* Endow the image with a local numbering: retrieve and sort its elements. */
1031: PetscHashISize(himage,n);
1032: PetscMalloc(n*sizeof(PetscInt), &image);
1033: PetscHashIGetKeys(himage,n,image);
1034: PetscSortInt(n,image);
1035: /* (C) */
1036: /*
1037: Convert ij to local numbering: insert image elements into an emptied and resized himage, mapping them to their local numbers.
1038: Then remap all of ij using himage.
1039: */
1040: PetscHashIClear(himage);
1041: PetscHashIResize(himage,n);
1042: for(j = 0; j < n; ++j) {
1043: PetscHashIAdd(himage,image[j],j);
1044: }
1045: totalnij = 0;
1046: PetscHashIMapArray(himage,pg->ijlen[pg->m],pg->ij,totalnij,pg->ij);
1047: if(totalnij!=pg->ijlen[pg->m]) {
1048: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of image indices before %D and after %D localization do not match", pg->ijlen[pg->m],totalnij);
1049: }
1050: /* Store the newly computed image array. */
1051: pg->image = image;
1052: /* Clean up. */
1053: PetscHashIDestroy(himage);
1054: return(0);
1055: }
1057: /*
1058: Indices are assumed sorted on ix.
1059: If !multivalued, remove iy duplicates from each ix's image segment.
1060: Record the number of images in ijlen.
1061: Store unique ix in a hash table along with their local numbers.
1062: Sort iy and store the unique in a (temporary) hash table to determine the local image.
1063: Note: this routine takes ownership of ("steals the reference to") ixidx and iyidx.
1064: */
1067: static PetscErrorCode MatIJSetEdgesLocal_Private(Mat A, const PetscInt len, PetscInt *ixidx, PetscInt *iyidx)
1068: {
1070: Mat_IJ *pg = (Mat_IJ *) A->data;
1071: PetscInt start, end, totalnij;
1072: PetscInt i,j,minnij, maxnij, nij,m,*ij, *ijlen, *supp;
1073: PetscHashI hsupp = PETSC_NULL;
1076: if(!len) {
1077: return(0);
1078: }
1080: m = 0;
1081: totalnij = 0;
1082: maxnij = 0;
1083: minnij = len;
1084: start = 0;
1085: while (start < len) {
1086: end = start+1;
1087: while (end < len && ixidx[end] == ixidx[start]) ++end;
1088: ++(m); /* count all of ixidx[start:end-1] as a single occurence of an idx index */
1089: nij = 1; /* At least one element in the image. */
1090: if (end - 1 > start) { /* found 2 or more of ixidx[start] in a row */
1091: /* sort the relevant portion of iy */
1092: PetscSortInt(end-start,iyidx+start);
1093: if(pg->multivalued) {
1094: nij = end-start;
1095: }
1096: else {
1097: /* count unique elements in iyidx[start,end-1] */
1098: for(j=start+1; j < end; ++j){
1099: if(iyidx[j] > iyidx[j-1]) ++nij;
1100: }
1101: }
1102: }
1103: totalnij += nij;
1104: minnij = PetscMin(minnij, nij);
1105: maxnij = PetscMax(maxnij, nij);
1106: start = end;
1107: }
1108: /*
1109: Now we know the size of the support -- m, and the total size of concatenated image segments -- totalnij.
1110: Allocate an array for recording the support indices -- supp.
1111: Allocate an array for recording the images of each support index -- ij.
1112: Allocate an array for counting the number of images for each support index -- ijlen.
1113: */
1114: PetscMalloc(sizeof(PetscInt)*(m+1), &ijlen);
1115: PetscMalloc(sizeof(PetscInt)*(totalnij),&ij);
1116: PetscMalloc(sizeof(PetscInt)*m, &supp);
1118: /*
1119: We now record in supp only the unique ixidx indices, and in ij the iyidx indices in each of the image segments.
1120: */
1121: if(m < A->rmap->n)
1122: PetscHashICreate(hsupp);
1123: i = 0;
1124: j = 0;
1125: start = 0;
1126: ijlen[0] = 0;
1127: while (start < len) {
1128: end = start+1;
1129: while (end < len && ixidx[end] == ixidx[start]) ++end;
1130: if(hsupp) {
1131: PetscHashIAdd(hsupp,ixidx[start],i);
1132: }
1133: ++i;
1134: /* the relevant portion of iy is already sorted. */
1135: ij[j++] = iyidx[start++];
1136: while(start < end) {
1137: if(pg->multivalued || iyidx[start] > iyidx[start-1])
1138: ij[j++] = iyidx[start];
1139: ++start;
1140: }
1141: ijlen[i] = j;
1142: }
1144: PetscFree(ixidx);
1145: PetscFree(iyidx);
1146: /* Record the changes. */
1147: pg->hsupp = hsupp;
1148: pg->m = m;
1149: pg->ij = ij;
1150: pg->ijlen = ijlen;
1151: pg->minijlen = minnij;
1152: pg->maxijlen = maxnij;
1154: return(0);
1155: }
1159: PetscErrorCode MatAssemblyBegin_IJ(Mat A, MatAssemblyType type)
1160: {
1161: Mat_IJ *ij = (Mat_IJ*)(A->data);
1162: PetscInt len, *ixidx, *iyidx;
1166: MatStashMPIIJAssemble_Private(ij->stash);
1167: if(type == MAT_FINAL_ASSEMBLY) {
1168: MatIJGetEdges(A, &len, &ixidx, &iyidx);
1169: MatStashMPIIJClear_Private(ij->stash);
1171: MatStashMPIIJSetPreallocation_Private(ij->stash, 0,0);
1172: MatIJClear_Private(A);
1173: MatIJSetEdgesLocal_Private(A, len, ixidx, iyidx);
1174: }
1176: return(0);
1177: }
1182: PetscErrorCode MatAssemblyEnd_IJ(Mat A, MatAssemblyType type)
1183: {
1185: /* Currently a noop */
1186: return(0);
1187: }
1190: /*@C
1191: MatIJGetSupport - retrieves the global indices of the graph's vertices of nonzero outdegree
1192: (i.e., the global indices of this processor's nonzero rows).
1193: If the graph is regarded as a multivalued map on integers, this is
1194: the support of the map (i.e., the set of indices with nonempty images).
1195:
1196:
1197: Not collective
1199: Input parameters:
1200: . A - pseudograph
1202: Output parameters:
1203: + len - the length of the support array
1204: - supp - the support array
1206: Note:
1207: + This operation fails for a nonassembled matrix.
1208: . In general, the returned indices are unsorted; use PetscSortInt, if necessary.
1209: - The caller is responsible for freeing the support array.
1212: Level: intermediate
1214: .seealso: MATIJ, MatIJGetSupportIS(), MatIJGetImage()
1215: @*/
1216: #undef __FUNCT__
1218: PetscErrorCode MatIJGetSupport(Mat A, PetscInt *len, PetscInt **supp)
1219: {
1221: Mat_IJ *pg = (Mat_IJ *)(A->data);
1222: PetscBool isij;
1225: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1226: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1227: MatIJCheckAssembled(A,PETSC_TRUE,1);
1228: if(!len && !supp) return(0);
1229: if(len) *len = pg->m;
1230: if(!supp) return(0);
1231: if(!*supp) {
1232: PetscMalloc(sizeof(PetscInt)*pg->m, supp);
1233: }
1234: if(!pg->hsupp) {
1235: PetscInt i;
1236: for(i = 0; i < pg->m; ++i) {
1237: (*supp)[i] = i + A->rmap->rstart;
1238: }
1239: }
1240: *len = 0;
1241: PetscHashIGetKeys(pg->hsupp, *len, *supp);
1242: return(0);
1243: }
1245: /*@C
1246: MatIJGetSupportIS - retrieves the global indices of the graph's vertices of nonzero outdegree
1247: (i.e., the global indices of this processor's nonzero rows).
1248: If the graph is regarded as a multivalued map on integers, this is
1249: the support of the map (i.e., the set of indices with nonempty images).
1250:
1251: Not collective
1253: Input parameters:
1254: . A - pseudograph
1256: Output parameters:
1257: . supp - the support IS
1259: Note:
1260: + This operation fails for a nonassembled matrix.
1261: - The caller is responsible for destroying the support IS.
1264: Level: intermediate
1266: .seealso: MATIJ, MatIJGetSupport(), MatIJGetSupportISSupported(), MatIJGetImage()
1267: @*/
1268: #undef __FUNCT__
1270: PetscErrorCode MatIJGetSupportIS(Mat A, IS *supp)
1271: {
1273: Mat_IJ *pg = (Mat_IJ *)(A->data);
1274: PetscBool isij;
1275: PetscInt ilen, *isupp;
1278: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1279: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1280: MatIJCheckAssembled(A,PETSC_TRUE,1);
1281: if(!supp) return(0);
1282: if(pg->hsupp) {
1283: MatIJGetSupport(A, &ilen, &isupp);
1284: ISCreateGeneral(A->rmap->comm, ilen, isupp, PETSC_OWN_POINTER, supp);
1285: }
1286: else {
1287: ISCreateStride(A->rmap->comm, A->rmap->n, A->rmap->rstart,1,supp);
1288: }
1289: return(0);
1290: }
1295: /*@C
1296: MatIJGetImage - retrieves the global indices of the graph's vertices of nonzero indegree
1297: on this processor (i.e., the global indices of this processor's nonzero columns).
1298: If the graph is regarded as a multivalued map on integers, this is
1299: the image of the map (the union of the images of this processor's support indices).
1300:
1301: Not collective
1303: Input parameters:
1304: . A - pseudograph
1306: Output parameters:
1307: + len - the length of the image array
1308: - image - the image array
1310: Note:
1311: + This operation fails for a nonassembled matrix.
1312: - The caller is responsible for freeing the image array.
1315: Level: intermediate
1317: .seealso: MATIJ, MatIJGetImageIS(), MatIJGetSupport(), MatIJGetMaxRowSize()
1318: @*/
1319: #undef __FUNCT__
1321: PetscErrorCode MatIJGetImage(Mat A, PetscInt *len, PetscInt **image)
1322: {
1324: Mat_IJ *pg = (Mat_IJ *)(A->data);
1325: PetscBool isij;
1328: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1329: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1330: MatIJCheckAssembled(A,PETSC_TRUE,1);
1331: MatIJLocalizeImage_Private(A);
1332: if(len) *len = pg->n;
1333: if(image) {
1334: PetscMalloc(sizeof(PetscInt)*pg->n, image);
1335: PetscMemcpy(*image, pg->image, sizeof(PetscInt)*pg->n);
1336: }
1337: return(0);
1338: }
1342: /*@C
1343: MatIJGetMaxRowSize - returns the largest number of nonzero columns in all of this processor's rows.
1344: If MatIJ (equivalently, the underlying graph) is regarded as a multivalued
1345: mapping on integers, then the result is the size of the largest set among
1346: the images of this processor's indices.
1347: Not collective.
1349: Input parameters:
1350: . A - pseudograph
1351:
1352: Output parameters:
1353: . maxsize - the size of the largest image set
1355: Level: advanced
1357: Notes:
1358: + This routine is useful for preallocating arrays to hold the images of the local indices:
1359: if an array of the largest image size is allocated, it can be used for repeatedly computing
1360: the images of the local indices.
1361: - This routine will fail for an unassembled matrix.
1363: .seealso: MATIJ, MatIJGetImage(), MatIJGetSupport(), MatIJMapI(), MatIJMapIJ(), MatIJMapIW(), MatIJMapIJW()
1364: @*/
1367: PetscErrorCode MatIJGetMaxRowSize(Mat A, PetscInt *maxsize)
1368: {
1370: Mat_IJ *pg = (Mat_IJ *)(A->data);
1371: PetscBool isij;
1374: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1375: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1376: MatIJCheckAssembled(A,PETSC_TRUE,1);
1377: *maxsize = pg->maxijlen;
1378: return(0);
1379: }
1381: /*@C
1382: MatIJGetMinRowSize - returns the largest number of nonzero columns in all of this processor's rows nonzero rows,
1383: or zero if no local nonzero rows exist.
1384: If MatIJ (equivalently, the underlying graph) is regarded as a multivalued
1385: mapping on integers, then the result is the size of the smallest nonempty set among
1386: the images of this processor's indices. If all images are empty, the result is zero.
1387: Not collective.
1389: Input parameters:
1390: . A - pseudograph
1391:
1392: Output parameters:
1393: . minsize - the size of the smallest nonempty image set
1395: Level: advanced
1397: .seealso: MatIJGetMinRowSize()
1398: @*/
1401: PetscErrorCode MatIJGetMinRowSize(Mat A, PetscInt *minsize)
1402: {
1404: Mat_IJ *pg = (Mat_IJ *)(A->data);
1405: PetscBool isij;
1408: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1409: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1410: MatIJCheckAssembled(A,PETSC_TRUE,1);
1411: *minsize = pg->minijlen;
1412: return(0);
1413: }
1418: #undef __FUNCT__
1420: PetscErrorCode MatDuplicate_IJ(Mat A, MatDuplicateOption op, Mat *B)
1421: {
1423: Mat_IJ* aij = (Mat_IJ*)(A->data), *bij;
1425: MatIJCheckAssembled(A,PETSC_TRUE,1);
1426: MatCreate(((PetscObject)A)->comm, B);
1427: MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
1428: MatSetType(*B, MATIJ);
1429: bij = (Mat_IJ*)((*B)->data);
1430: if(aij->hsupp) PetscHashIDuplicate(aij->hsupp, bij->hsupp);
1431: bij->m = aij->m;
1432: PetscMemcpy(bij->ijlen, aij->ijlen, sizeof(PetscInt)*(bij->m+1));
1433: PetscMemcpy(bij->ij, aij->ij, sizeof(PetscInt)*(bij->ijlen[bij->m]));
1434: bij->n = aij->n;
1435: if(aij->image) {
1436: PetscMemcpy(bij->image, aij->image, sizeof(PetscInt)*bij->n);
1437: }
1438: bij->maxijlen = aij->maxijlen;
1439: bij->minijlen = aij->minijlen;
1440: return(0);
1441: }
1443: /*@C
1444: MatIJGetImageIS - retrieves the global indices of the graph's vertices of nonzero indegree
1445: on this processor (i.e., the global indices of this processor's nonzero columns).
1446: If the graph is regarded as a multivalued map on integers, this is
1447: the image of the map (the union of the images of this processor's support indices).
1448:
1449: Not collective
1451: Input parameters:
1452: . A - pseudograph
1454: Output parameters:
1455: . image - image IS
1457: Note:
1458: + This operation fails for a nonassembled matrix.
1459: - The caller is responsible for freeing the image IS.
1462: Level: intermediate
1464: .seealso: MatIJGetImage(), MatIJGetSupportIS(), MatIJGetMaxRowSize()
1465: @*/
1466: #undef __FUNCT__
1468: PetscErrorCode MatIJGetImageIS(Mat A, IS *image)
1469: {
1471: Mat_IJ *pg = (Mat_IJ *)(A->data);
1472: PetscBool isij;
1475: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1476: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1477: MatIJCheckAssembled(A,PETSC_TRUE,1);
1478: if(image) {
1479: ISCreateGeneral(A->cmap->comm, pg->n, pg->image, PETSC_COPY_VALUES, image);
1480: }
1481: return(0);
1482: }
1485: /*@C
1486: MatIJGetRowSizes - retrieves the numbers of edges emanating from the each of the supplied global indices,
1487: provided they fall into the local ownership range. Other indices result in an error.
1488:
1489: Not collective
1491: Input parameters:
1492: + A - pseudograph
1493: . intype - (MATIJ_LOCAL | MATIJ_GLOBAL) meaning of inidxi: local support numbers or global indices
1494: . len - the length of the index list, or PETSC_DETERMINE to indicate all of the locally supported indices
1495: - inidxi - array (of length len) of global indices whose image sizes are sought; ignored if len == PETSC_DETERMINE,
1496: which is equivalent to using all of the supported indices.
1498: Output parameters:
1499: . sizes - array (of length len) of image sizes of the global indices in inidxi
1501: Note:
1502: + This operation fails for a nonassembled matrix.
1503: . If len is PETSC_DEFAULT, inidxi must be PETSC_DEFAULT, and vice versa.
1504: - The caller is responsible for freeing sizes.
1507: Level: intermediate
1509: .seealso: MatIJ, MatIJGetRowSizesSupported(), MatIJGetImage(), MatIJGetImageIS(), MatIJGetMaxRowSize()
1510: @*/
1511: #undef __FUNCT__
1513: PetscErrorCode MatIJGetRowSizes(Mat A, MatIJIndexType intype, PetscInt len, const PetscInt *inidxi, PetscInt **sizes)
1514: {
1516: PetscBool isij;
1519: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1520: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1521: MatIJCheckAssembled(A,PETSC_TRUE,1);
1522: MatIJMap(A,intype,len,inidxi,PETSC_NULL,PETSC_NULL,MATIJ_GLOBAL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,sizes);
1523: return(0);
1524: }
1526: /*@C
1527: MatIJGetSupportSize - retrieves the total numbers of nonzero outdegree indices in the local ownership range
1528: (the number of nonzero local rows).
1529:
1530: Not collective.
1532: Input parameters:
1533: . A - pseudograph
1535: Output parameters:
1536: . size - local support size
1538: Note:
1539: . This operation fails for a nonassembled matrix.
1542: Level: intermediate
1544: .seealso: MATIJ, MatIJGetRowSizes(), MatIJGetImage(), MatIJGetImageIS(), MatIJGetMinRowSize(), MatIJGetMaxRowSize(), MatIJGetSupportSize()
1545: @*/
1546: #undef __FUNCT__
1548: PetscErrorCode MatIJGetSupportSize(Mat A, PetscInt *size)
1549: {
1551: PetscBool isij;
1552: Mat_IJ *pg = (Mat_IJ*)A->data;
1555: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1556: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1558: MatIJCheckAssembled(A,PETSC_TRUE,1);
1559: *size = pg->m;
1560: return(0);
1561: }
1564: /*@C
1565: MatIJGetImageSize - retrieves the total numbers of target indices adjacent to the source indices in the local ownership
1566: range (the number of nonzero local columns).
1567:
1568: Not collective.
1570: Input parameters:
1571: . A - pseudograph
1573: Output parameters:
1574: . size - local image size
1576: Note:
1577: . This operation fails for a nonassembled matrix.
1580: Level: intermediate
1582: .seealso: MATIJ, MatIJGetSupportSize(), MatIJGetSupport(), MatIJGetSupportIS()
1583: @*/
1584: #undef __FUNCT__
1586: PetscErrorCode MatIJGetImageSize(Mat A, PetscInt *size)
1587: {
1589: PetscBool isij;
1590: Mat_IJ *pg = (Mat_IJ*)A->data;
1593: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
1594: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
1596: MatIJCheckAssembled(A,PETSC_TRUE,1);
1597: MatIJLocalizeImage_Private(A);
1598: *size = pg->n;
1599: return(0);
1600: }
1602: #undef __FUNCT__
1604: PetscErrorCode MatIJBinRenumberLocal_Private(Mat A, MatIJIndexType intype, PetscInt insize, const PetscInt *inidxi, PetscInt *_outsize, PetscInt **_outidxi, PetscInt **_binsizes)
1605: {
1607: PetscInt indi = -1, i,j,k,outsize = -1, *outidxi = PETSC_NULL, *binsizes = PETSC_NULL;
1608: Mat_IJ *pg = (Mat_IJ*)A->data;
1612: /* We bin all of the input, and in the process assign bin-wise local numbers to them. */
1613: /* We'll need to count the contributions to each "bin" and the offset of each bin in outidx. */
1614: /* Binning requires a localized image. */
1615: MatIJLocalizeImage_Private(A);
1616: if((_outidxi && !*_outidxi)) {
1617: MatIJBinRenumberLocal_Private(A,intype,insize,inidxi,&outsize,PETSC_NULL,PETSC_NULL);
1618: }
1619: if(insize == PETSC_DETERMINE){
1620: insize = pg->m;
1621: inidxi = PETSC_NULL;
1622: }
1623: else if(insize < 0)
1624: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid input array size: %D", insize);
1625: if(_outidxi) {
1626: if(!*_outidxi) {
1627: PetscMalloc(sizeof(PetscInt)*outsize, _outidxi);
1628: }
1629: outidxi = *_outidxi;
1630: }
1631: if(_binsizes) {
1632: if(!*_binsizes) {
1633: PetscMalloc(sizeof(PetscInt)*(pg->n), _binsizes);
1634: }
1635: binsizes = *_binsizes;
1636: }
1637: /* Allocate the bin offset array, if necessary. */
1638: if(!pg->binoffsets) {
1639: PetscMalloc((pg->n+1)*sizeof(PetscInt), &(pg->binoffsets));
1640: }
1641: /* Initialize bin offsets */
1642: for(j = 0; j <= pg->n; ++j) {
1643: pg->binoffsets[j] = 0;
1644: }
1646: for(i = 0; i < insize; ++i) {
1647: if(!inidxi) {
1648: indi = i;
1649: }
1650: else {
1651: /* Convert to local. */
1652: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
1653: }
1654: if((indi < 0 || indi > pg->m)){
1655: /* drop */
1656: continue;
1657: }
1658: for(k = pg->ijlen[indi]; k < pg->ijlen[indi+1]; ++k) {
1659: ++(pg->binoffsets[pg->ij[k]+1]);
1660: }
1661: }/* for(i = 0; i < insize; ++i) */
1662: /* Convert bin sizes into bin offsets. */
1663: for(j = 0; j < pg->n; ++j) {
1664: pg->binoffsets[j+1] += pg->binoffsets[j];
1665: }
1666: /* Now bin the input indices and values. */
1667: if(outidxi) {
1668: /* Allocate the bin size array, if necessary. */
1669: if(!binsizes) {
1670: if(!pg->binsizes) {
1671: PetscMalloc((pg->n)*sizeof(PetscInt), &(pg->binsizes));
1672: }
1673: binsizes = pg->binsizes;
1674: }
1675: /* Initialize bin sizes to zero. */
1676: for(j = 0; j < pg->n; ++j) {
1677: binsizes[j] = 0;
1678: }
1679: for(i = 0; i < insize; ++i) {
1680: if(!inidxi) {
1681: indi = i;
1682: }
1683: else {
1684: /* Convert to local. */
1685: MatIJGetSuppIndex_Private(A,intype,inidxi[i],indi);
1686: if((indi < 0 || indi > pg->m)){
1687: /* drop */
1688: continue;
1689: }
1690: }
1691: for(k = pg->ijlen[indi]; k < pg->ijlen[indi+1]; ++k) {
1692: j = pg->ij[k];
1693: outidxi[pg->binoffsets[j]+binsizes[j]] = binsizes[j];
1694: ++binsizes[j];
1695: }
1696: }/* for(i = 0; i < insize; ++i) */
1697: }/* if(outidxi) */
1698: if(_outsize) *_outsize = pg->binoffsets[pg->n];
1699: return(0);
1700: }
1703: /*@C
1704: MatIJBinRenumber - map the support indices to their global numbers within their image bins.
1705: If the image indices are interpreted as colors labeling subdomains, then
1706: each global subdomain is given a new contiguous zero-based numbering
1707: uniquely defined by the following: the new vertex numbers increase with the
1708: owning processor's rank, and within each rank they are arranged according
1709: to their order in the local portion of the bin.
1712: Collective on A.
1714: Input arguments:
1715: . A - pseudograph
1716:
1717: Output arguments:
1718: . B - renumbering pseudograph
1720: Level: advanced
1722: Notes: observe that each local support index might be mapped to the same global index multiple times,
1723: if it happens to have the same number within different bins. In order to decide which color
1724: each of the new numbers refers to, it is useful to use the result B in conjunction with
1725: the original pseudograph A as the second and first argument to MatIJBinMap(), respectively.
1726: By construction, B is compatible to A in the sense of MatIJBinMap().
1728: .keywords: pseudograph, coloring, binning, numbering
1729: .seealso: MatIJBinMap()
1730: @*/
1731: #undef __FUNCT__
1733: PetscErrorCode MatIJBinRenumber(Mat A, Mat *B)
1734: {
1736: Mat_IJ *pg = (Mat_IJ*)A->data;
1737: PetscInt iysize = -1, len, *ixidx, *iyidx = PETSC_NULL, *bsizes = PETSC_NULL, N = 0, g,b,blow,bhigh,j;
1738: PetscMPIInt rank, size, NN, bsize,goffset;
1741: MatIJLocalizeImage_Private(A);
1742: MatIJBinRenumberLocal_Private(A, MATIJ_LOCAL, pg->m, PETSC_NULL, &iysize, &iyidx, &bsizes);
1743: MPI_Comm_size(((PetscObject)A)->comm, &size);
1744: MPI_Comm_rank(((PetscObject)A)->comm, &rank);
1745: /*
1746: Since the new mapping is to the global bin numberings, we need to adjust the numberings further,
1747: by determining the local offsets for each bin: the sizes of each bin on the preceeding
1748: processors within the communicator. This is accomplished by a series of scan ops, one per bin.
1749: */
1750: /*
1751: In addition, we need to know the largest global bin size N, which can be determined by looking at
1752: bin offset of the process with the largest rank in the comm, adding the local bin size to it,
1753: taking the max across the bins and then broadcasting it to the other processes in the comm.
1754: */
1755: /* Loop over the global bins g; count the locally-supported bins b. */
1756: b = 0;
1757: blow = bhigh = 0;
1758: for(g = 0, b = 0; g < A->cmap->N; ++g) {
1759: if(pg->image[b] == g) {
1760: bsize = PetscMPIIntCast(bsizes[b]);
1761: }
1762: else {
1763: bsize = 0;
1764: }
1765: MPI_Scan(&bsize,&goffset,1,MPI_INT, MPI_SUM,((PetscObject)A)->comm);
1766: if(pg->image[b] == g) {
1767: blow = bhigh;
1768: bhigh = blow + bsizes[b];
1769: /* Shift the indices by the goffset. */
1770: for(j = blow; j < bhigh; ++j) iyidx[j] += goffset;
1771: /* Compute the maximum partial global bin size, up to and including this proc's portion. */
1772: NN = PetscMPIIntCast(PetscMax(NN,goffset + bsizes[b]));
1773: ++b;
1774: }
1775: else {
1776: /* Compute the maximum partial global bin size, up to and including this proc's portion. */
1777: NN = PetscMPIIntCast(PetscMax(NN,goffset));
1778: }
1779: }/* Loop over the global bins. */
1780: if(b != pg->n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Traversed %D local bins, not the same as expected: %D", b, pg->n);
1781: /* Broadcast from the last rank the largest global bin size. */
1782: MPI_Bcast(&NN,1,MPI_INT, rank-1,((PetscObject)A)->comm);
1783: N = NN;
1784: /* Now construct the new pseudograph. */
1785: /* The number of edges and the source vertices are the same as in the old pseudograph. */
1786: MatIJGetEdges(A,&len,&ixidx,PETSC_NULL);
1787: /* But instead of terminating at bins, the edges terminate at the local numbers within the bin. */
1788: #if defined PETSC_USE_DEBUG
1789: {
1790: PetscInt k,blen = 0;
1791: for(k = 0; k < pg->n; ++k) blen += bsizes[k];
1792: if(len != blen)
1793: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of edges in the original pseudograph %D and the renumbering pseudograph %D do not match", len, blen);
1794: }
1795: #endif
1796: MatCreate(((PetscObject)A)->comm, B);
1797: MatSetSizes(*B, A->rmap->n, PETSC_DETERMINE, PETSC_DETERMINE, N);
1798: MatIJSetEdges(*B, len, ixidx, iyidx);
1799: /* All ixidx indices are within the local ownership range, so no parallel assembly is required. */
1800: MatIJSetEdgesLocal_Private(*B, len, ixidx, iyidx);
1801: return(0);
1802: }
1804: #undef __FUNCT__
1806: PetscErrorCode MatTranspose_IJ(Mat A, MatReuse reuse, Mat *B)
1807: {
1809: PetscBool multivalued;
1810: IS ix, iy;
1812: MatCreate(((PetscObject)A)->comm, B);
1813: MatSetSizes(*B, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N);
1814: MatSetType(*B, MATIJ);
1815: MatIJGetMultivalued(A,&multivalued);
1816: MatIJGetEdgesIS(A, &ix, &iy);
1817: MatIJSetEdgesIS(*B, iy,ix);
1818: ISDestroy(&(ix));
1819: ISDestroy(&(iy));
1820: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
1821: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
1822: return(0);
1823: }
1825: #undef __FUNCT__
1827: PetscErrorCode MatTransposeMatMult_IJ_IJ(Mat A, Mat B, MatReuse reuse, PetscReal fill, Mat *CC)
1828: {
1830: Mat C;
1831: PetscInt nsupp1, nsupp2, nsupp3, *supp1 = PETSC_NULL, *supp2 = PETSC_NULL, *supp3, imgsize1, *imgsizes1 = PETSC_NULL,
1832: imgsize2, *imgsizes2 = PETSC_NULL, *image1 = PETSC_NULL, *image2 = PETSC_NULL,
1833: *ixidx, *iyidx, count, i1,i2,i1low,i1high,i2low,i2high,k;
1836: /*
1837: C _
1838: ... |
1839: |-----| |-----| ------> |
1840: ^ ^ |
1841: | ======> | -
1842: A | A |
1843: ... | B _ ... |
1844: | ... | |
1845: |-----| ------> | |-----|
1846: |
1847: -
1848: */
1849: MatIJGetSupport(A, &nsupp1, &supp1);
1850: MatIJGetSupport(B, &nsupp2, &supp2);
1851: /* Avoid computing the intersection, which may be unscalable in storage. */
1852: /*
1853: Count the number of images of the intersection of supports under the "upward" (A) and "rightward" (B) maps.
1854: It is done this way: supp1 is mapped by B obtaining offsets2, and supp2 is mapped by A obtaining offsets1.
1855: */
1856: MatIJMap(A,MATIJ_GLOBAL,nsupp2,supp2,PETSC_NULL,PETSC_NULL,MATIJ_GLOBAL,&imgsize1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&imgsizes1);
1857: MatIJMap(B,MATIJ_GLOBAL,nsupp1,supp1,PETSC_NULL,PETSC_NULL,MATIJ_GLOBAL,&imgsize2,PETSC_NULL,PETSC_NULL,PETSC_NULL,&imgsizes2);
1858: /* Count the number of supp1 indices with nonzero images in B -- that's the size of the intersection. */
1859: nsupp3 = 0;
1860: for(k = 0; k < nsupp1; ++k) nsupp3 += (imgsizes1[k]>0);
1861: #if defined(PETSC_USE_DEBUG)
1862: /* Now count the number of supp2 indices with nonzero images in map1: should be the same. */
1863: {
1864: PetscInt nsupp3_2 = 0;
1865: for(k = 0; k < nsupp2; ++k) nsupp3_2 += (imgsizes2[k]>0);
1866: if(nsupp3 != nsupp3_2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Intersections supports different in map1: %D and map2: %D", nsupp3, nsupp3_2);
1867: }
1868: #endif
1869: /* Allocate indices for the intersection. */
1870: PetscMalloc(sizeof(PetscInt)*nsupp3, &supp3);
1871: nsupp3 = 0;
1872: for(k = 0; k < nsupp2; ++k) {
1873: if(imgsizes1[k]) {
1874: supp3[nsupp3] = supp2[k];
1875: ++nsupp3;
1876: }
1877: }
1878: PetscFree(supp1);
1879: PetscFree(supp2);
1880:
1881: /*
1882: Now obtain the "up" (A) and "right" (B) images of supp3.
1883: Recall that imgsizes1 are allocated for lsupp2, and imgsizes2 for lsupp1.
1884: */
1885: /* Reuse imgsizes1 and imgsizes2, even though they are a bit bigger, than necessary now and full of junk from the previous calls. Saves malloc/free.*/
1886: MatIJMap(A,MATIJ_GLOBAL,nsupp3,supp3,PETSC_NULL,PETSC_NULL,MATIJ_GLOBAL,&imgsize1,&image1,PETSC_NULL,PETSC_NULL,&imgsizes1);
1887: MatIJMap(B,MATIJ_GLOBAL,nsupp3,supp3,PETSC_NULL,PETSC_NULL,MATIJ_GLOBAL,&imgsize2,&image2,PETSC_NULL,PETSC_NULL,&imgsizes2);
1888: PetscFree(supp3);
1890: /* Count the total number of arrows to add to the pushed forward MatIJ. */
1891: count = 0;
1892: for(k = 0; k < nsupp3; ++k) {
1893: count += (imgsizes1[k])*(imgsizes2[k]);
1894: }
1895: /* Allocate storage for the composed indices. */
1896: PetscMalloc2(count, PetscInt, &ixidx, count, PetscInt, &iyidx);
1897: count= 0;
1898: i1low = 0;
1899: i2low = 0;
1900: for(k = 0; k < nsupp3; ++k) {
1901: i1high = i1low + imgsizes1[k];
1902: i2high = i2low + imgsizes2[k];
1903: for(i1 = i1low; i1 < i1high; ++i1) {
1904: for(i2 = i2low; i1 < i2high; ++i2) {
1905: ixidx[count] = image1[i1];
1906: iyidx[count] = image2[i2];
1907: ++count;
1908: }
1909: }
1910: i1low = i1high;
1911: i2low = i2high;
1912: }
1913: PetscFree(image1);
1914: PetscFree(image2);
1915: PetscFree(imgsizes1);
1916: PetscFree(imgsizes2);
1917: /* Now construct the new MatIJ. */
1918: MatCreate(((PetscObject)A)->comm, &C);
1919: MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N);
1920: MatSetType(C, MATIJ);
1921: MatIJSetEdges(C,count,ixidx,iyidx);
1922: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
1923: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
1924: PetscFree2(ixidx,iyidx);
1926: *CC = C;
1927: return(0);
1928: }
1934: #undef __FUNCT__
1936: PetscErrorCode MatMatMult_IJ_IJ(Mat A, Mat B, MatReuse reuse, PetscReal fill, Mat *CC)
1937: {
1939: Mat At,C;
1943: /*
1945: B _ B _
1946: ... | ... |
1947: |-----| ------> | |-----| ------> |
1948: ^ | ^ |
1949: | - | -
1950: A | ======> A |
1951: ... | ... | C _
1952: | | ... |
1953: |-----| |-----| ------> |
1954: |
1955: -
1956: Convert this to a pushforward by transposing A to At and then pushing B forward along At
1957: (reflect the second diagram with respect to a horizontal axis and then compare with Pushforward,
1958: of just push forward "downward".)
1960: B _ B _ B _
1961: ... | ... | ... |
1962: |-----| ------> | |-----| ------> | |-----| ------> |
1963: ^ | | | | |
1964: | - | - | -
1965: A | ======> At | ======> At |
1966: ... | ... | ... | C _
1967: | V V ... |
1968: |-----| |-----| |-----| ------> |
1969: |
1970: -
1971: */
1972: MatTranspose(A, MAT_INITIAL_MATRIX, &At);
1973: MatTransposeMatMult(At, B, MAT_INITIAL_MATRIX, 1.0, &C);
1974: MatDestroy(&At);
1975: C->ops->matmult = MatMatMult_IJ_IJ;
1976: *CC = C;
1977: return(0);
1978: }
1981: #undef __FUNCT__
1983: PetscErrorCode MatZeroEntries_IJ(Mat A)
1984: {
1985: Mat_IJ *pg = (Mat_IJ*) A->data;
1988: MatIJClear_Private(A);
1989: MatStashMPIIJClear_Private(pg->stash);
1990: A->was_assembled = A->assembled;
1991: A->assembled = PETSC_FALSE;
1992: return(0);
1993: }
1996: EXTERN_C_BEGIN
1997: #undef __FUNCT__
1999: PetscErrorCode MatView_IJ(Mat A, PetscViewer v)
2000: {
2001: Mat_IJ *pg = (Mat_IJ*) A->data;
2002: PetscBool isij, isascii;
2003: PetscInt indi, indj,i=-1,j;
2004: PetscHashIIter it=0;
2008: PetscObjectTypeCompare((PetscObject)A,MATIJ,&isij);
2009: if(!isij) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix not of type MATIJ: %s", ((PetscObject)A)->type);
2010: MatIJCheckAssembled(A,PETSC_TRUE,1);
2011: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&isascii);
2012: if (!isascii) {
2013: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)v)->type_name);
2014: }
2015: if(!pg->hsupp) {
2016: i = A->rmap->rstart;
2017: }
2018: else {
2019: PetscHashIIterBegin(pg->hsupp,it);
2020: }
2021: while(1) {
2022: if(pg->hsupp) {
2023: if(PetscHashIIterAtEnd(pg->hsupp,it)) break;
2024: else {
2025: PetscHashIIterGetKeyVal(pg->hsupp,it,i,indi);
2026: }
2027: }
2028: else {
2029: if(i == A->rmap->rend)
2030: break;
2031: else
2032: indi = i - A->rmap->rstart;
2033: }
2034: PetscViewerASCIISynchronizedPrintf(v, "%D --> ", i);
2035: for(indj = pg->ijlen[indi]; indj < pg->ijlen[indi+1]; ++indj) {
2036: MatIJGetIndexImage_Private(A,MATIJ_GLOBAL,indj,j);
2037: PetscViewerASCIISynchronizedPrintf(v, "%D ", j);
2038: }
2039: }
2040: PetscViewerFlush(v);
2041: return(0);
2042: }
2045: #undef __FUNCT__
2047: PetscErrorCode MatDestroy_IJ(Mat A) {
2049: Mat_IJ *pg = (Mat_IJ *)(A->data);
2050:
2052: MatIJClear_Private(A);
2053: MatStashMPIIJDestroy_Private(&(pg->stash));
2054: A->data = PETSC_NULL;
2055:
2056: PetscObjectComposeFunction((PetscObject)A,"MatMatMult_ij_ij_C", "",PETSC_NULL);
2057: PetscObjectComposeFunction((PetscObject)A,"MatTransposeMatMult_ij_ij_C", "",PETSC_NULL);
2058: return(0);
2059: }
2061: #undef __FUNCT__
2063: PetscErrorCode MatCreate_IJ(Mat A)
2064: {
2066: Mat_IJ *pg;
2069: PetscNewLog(A, Mat_IJ, &pg);
2070: A->data = (void*)pg;
2072: PetscLayoutSetUp(A->rmap);
2073: PetscLayoutSetUp(A->cmap);
2075: MatStashMPIIJCreate_Private(A->rmap, &(pg->stash));
2077: PetscMemzero(A->ops,sizeof(*(A->ops)));
2078: A->ops->assemblybegin = MatAssemblyBegin_IJ;
2079: A->ops->assemblyend = MatAssemblyEnd_IJ;
2080: A->ops->zeroentries = MatZeroEntries_IJ;
2081: A->ops->getrow = MatGetRow_IJ;
2082: A->ops->restorerow = MatRestoreRow_IJ;
2083: A->ops->duplicate = MatDuplicate_IJ;
2084: A->ops->destroy = MatDestroy_IJ;
2085: A->ops->view = MatView_IJ;
2087: MatRegisterOp(((PetscObject)A)->comm, PETSC_NULL, (PetscVoidFunction)MatMatMult_IJ_IJ, "MatMatMult", 2, MATIJ,MATIJ);
2088: MatRegisterOp(((PetscObject)A)->comm, PETSC_NULL, (PetscVoidFunction)MatTransposeMatMult_IJ_IJ, "MatTransposeMatMult", 2, MATIJ,MATIJ);
2089: PetscObjectChangeTypeName((PetscObject)A, MATIJ);
2090: return(0);
2091: }
2092: EXTERN_C_END