Actual source code: mpibaij.c
1: #define PETSCMAT_DLL
3: #include src/mat/impls/baij/mpi/mpibaij.h
5: EXTERN PetscErrorCode MatSetUpMultiply_MPIBAIJ(Mat);
6: EXTERN PetscErrorCode DisAssemble_MPIBAIJ(Mat);
7: EXTERN PetscErrorCode MatIncreaseOverlap_MPIBAIJ(Mat,PetscInt,IS[],PetscInt);
8: EXTERN PetscErrorCode MatGetSubMatrices_MPIBAIJ(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *[]);
9: EXTERN PetscErrorCode MatGetValues_SeqBAIJ(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt [],PetscScalar []);
10: EXTERN PetscErrorCode MatSetValues_SeqBAIJ(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt [],const PetscScalar [],InsertMode);
11: EXTERN PetscErrorCode MatSetValuesBlocked_SeqBAIJ(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
12: EXTERN PetscErrorCode MatGetRow_SeqBAIJ(Mat,PetscInt,PetscInt*,PetscInt*[],PetscScalar*[]);
13: EXTERN PetscErrorCode MatRestoreRow_SeqBAIJ(Mat,PetscInt,PetscInt*,PetscInt*[],PetscScalar*[]);
14: EXTERN PetscErrorCode MatZeroRows_SeqBAIJ(Mat,PetscInt,const PetscInt[],PetscScalar);
16: /* UGLY, ugly, ugly
17: When MatScalar == PetscScalar the function MatSetValuesBlocked_MPIBAIJ_MatScalar() does
18: not exist. Otherwise ..._MatScalar() takes matrix elements in single precision and
19: inserts them into the single precision data structure. The function MatSetValuesBlocked_MPIBAIJ()
20: converts the entries into single precision and then calls ..._MatScalar() to put them
21: into the single precision data structures.
22: */
23: #if defined(PETSC_USE_MAT_SINGLE)
24: EXTERN PetscErrorCode MatSetValuesBlocked_SeqBAIJ_MatScalar(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const MatScalar*,InsertMode);
25: EXTERN PetscErrorCode MatSetValues_MPIBAIJ_MatScalar(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const MatScalar*,InsertMode);
26: EXTERN PetscErrorCode MatSetValuesBlocked_MPIBAIJ_MatScalar(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const MatScalar*,InsertMode);
27: EXTERN PetscErrorCode MatSetValues_MPIBAIJ_HT_MatScalar(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const MatScalar*,InsertMode);
28: EXTERN PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT_MatScalar(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const MatScalar*,InsertMode);
29: #else
30: #define MatSetValuesBlocked_SeqBAIJ_MatScalar MatSetValuesBlocked_SeqBAIJ
31: #define MatSetValues_MPIBAIJ_MatScalar MatSetValues_MPIBAIJ
32: #define MatSetValuesBlocked_MPIBAIJ_MatScalar MatSetValuesBlocked_MPIBAIJ
33: #define MatSetValues_MPIBAIJ_HT_MatScalar MatSetValues_MPIBAIJ_HT
34: #define MatSetValuesBlocked_MPIBAIJ_HT_MatScalar MatSetValuesBlocked_MPIBAIJ_HT
35: #endif
39: PetscErrorCode MatGetRowMaxAbs_MPIBAIJ(Mat A,Vec v,PetscInt idx[])
40: {
41: Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data;
43: PetscInt i,*idxb = 0;
44: PetscScalar *va,*vb;
45: Vec vtmp;
48:
49: MatGetRowMaxAbs(a->A,v,idx);
50: VecGetArray(v,&va);
51: if (idx) {
52: for (i=0; i<A->cmap.n; i++) {if (PetscAbsScalar(va[i])) idx[i] += A->cmap.rstart;}
53: }
55: VecCreateSeq(PETSC_COMM_SELF,A->rmap.n,&vtmp);
56: if (idx) {PetscMalloc(A->rmap.n*sizeof(PetscInt),&idxb);}
57: MatGetRowMaxAbs(a->B,vtmp,idxb);
58: VecGetArray(vtmp,&vb);
60: for (i=0; i<A->rmap.n; i++){
61: if (PetscAbsScalar(va[i]) < PetscAbsScalar(vb[i])) {va[i] = vb[i]; if (idx) idx[i] = A->cmap.bs*a->garray[idxb[i]/A->cmap.bs] + (idxb[i] % A->cmap.bs);}
62: }
64: VecRestoreArray(v,&va);
65: VecRestoreArray(vtmp,&vb);
66: if (idxb) {PetscFree(idxb);}
67: VecDestroy(vtmp);
68: return(0);
69: }
74: PetscErrorCode MatStoreValues_MPIBAIJ(Mat mat)
75: {
76: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;
80: MatStoreValues(aij->A);
81: MatStoreValues(aij->B);
82: return(0);
83: }
89: PetscErrorCode MatRetrieveValues_MPIBAIJ(Mat mat)
90: {
91: Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;
95: MatRetrieveValues(aij->A);
96: MatRetrieveValues(aij->B);
97: return(0);
98: }
101: /*
102: Local utility routine that creates a mapping from the global column
103: number to the local number in the off-diagonal part of the local
104: storage of the matrix. This is done in a non scable way since the
105: length of colmap equals the global matrix length.
106: */
109: PetscErrorCode CreateColmap_MPIBAIJ_Private(Mat mat)
110: {
111: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
112: Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)baij->B->data;
114: PetscInt nbs = B->nbs,i,bs=mat->rmap.bs;
117: #if defined (PETSC_USE_CTABLE)
118: PetscTableCreate(baij->nbs,&baij->colmap);
119: for (i=0; i<nbs; i++){
120: PetscTableAdd(baij->colmap,baij->garray[i]+1,i*bs+1);
121: }
122: #else
123: PetscMalloc((baij->Nbs+1)*sizeof(PetscInt),&baij->colmap);
124: PetscLogObjectMemory(mat,baij->Nbs*sizeof(PetscInt));
125: PetscMemzero(baij->colmap,baij->Nbs*sizeof(PetscInt));
126: for (i=0; i<nbs; i++) baij->colmap[baij->garray[i]] = i*bs+1;
127: #endif
128: return(0);
129: }
131: #define CHUNKSIZE 10
133: #define MatSetValues_SeqBAIJ_A_Private(row,col,value,addv) \
134: { \
135: \
136: brow = row/bs; \
137: rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; \
138: rmax = aimax[brow]; nrow = ailen[brow]; \
139: bcol = col/bs; \
140: ridx = row % bs; cidx = col % bs; \
141: low = 0; high = nrow; \
142: while (high-low > 3) { \
143: t = (low+high)/2; \
144: if (rp[t] > bcol) high = t; \
145: else low = t; \
146: } \
147: for (_i=low; _i<high; _i++) { \
148: if (rp[_i] > bcol) break; \
149: if (rp[_i] == bcol) { \
150: bap = ap + bs2*_i + bs*cidx + ridx; \
151: if (addv == ADD_VALUES) *bap += value; \
152: else *bap = value; \
153: goto a_noinsert; \
154: } \
155: } \
156: if (a->nonew == 1) goto a_noinsert; \
157: if (a->nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) into matrix", row, col); \
158: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,aimax,a->nonew,MatScalar); \
159: N = nrow++ - 1; \
160: /* shift up all the later entries in this row */ \
161: for (ii=N; ii>=_i; ii--) { \
162: rp[ii+1] = rp[ii]; \
163: PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar)); \
164: } \
165: if (N>=_i) { PetscMemzero(ap+bs2*_i,bs2*sizeof(MatScalar)); } \
166: rp[_i] = bcol; \
167: ap[bs2*_i + bs*cidx + ridx] = value; \
168: a_noinsert:; \
169: ailen[brow] = nrow; \
170: }
172: #define MatSetValues_SeqBAIJ_B_Private(row,col,value,addv) \
173: { \
174: brow = row/bs; \
175: rp = bj + bi[brow]; ap = ba + bs2*bi[brow]; \
176: rmax = bimax[brow]; nrow = bilen[brow]; \
177: bcol = col/bs; \
178: ridx = row % bs; cidx = col % bs; \
179: low = 0; high = nrow; \
180: while (high-low > 3) { \
181: t = (low+high)/2; \
182: if (rp[t] > bcol) high = t; \
183: else low = t; \
184: } \
185: for (_i=low; _i<high; _i++) { \
186: if (rp[_i] > bcol) break; \
187: if (rp[_i] == bcol) { \
188: bap = ap + bs2*_i + bs*cidx + ridx; \
189: if (addv == ADD_VALUES) *bap += value; \
190: else *bap = value; \
191: goto b_noinsert; \
192: } \
193: } \
194: if (b->nonew == 1) goto b_noinsert; \
195: if (b->nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) into matrix", row, col); \
196: MatSeqXAIJReallocateAIJ(B,b->mbs,bs2,nrow,brow,bcol,rmax,ba,bi,bj,rp,ap,bimax,b->nonew,MatScalar); \
197: CHKMEMQ;\
198: N = nrow++ - 1; \
199: /* shift up all the later entries in this row */ \
200: for (ii=N; ii>=_i; ii--) { \
201: rp[ii+1] = rp[ii]; \
202: PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar)); \
203: } \
204: if (N>=_i) { PetscMemzero(ap+bs2*_i,bs2*sizeof(MatScalar));} \
205: rp[_i] = bcol; \
206: ap[bs2*_i + bs*cidx + ridx] = value; \
207: b_noinsert:; \
208: bilen[brow] = nrow; \
209: }
211: #if defined(PETSC_USE_MAT_SINGLE)
214: PetscErrorCode MatSetValues_MPIBAIJ(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
215: {
216: Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)mat->data;
218: PetscInt i,N = m*n;
219: MatScalar *vsingle;
222: if (N > b->setvalueslen) {
223: PetscFree(b->setvaluescopy);
224: PetscMalloc(N*sizeof(MatScalar),&b->setvaluescopy);
225: b->setvalueslen = N;
226: }
227: vsingle = b->setvaluescopy;
229: for (i=0; i<N; i++) {
230: vsingle[i] = v[i];
231: }
232: MatSetValues_MPIBAIJ_MatScalar(mat,m,im,n,in,vsingle,addv);
233: return(0);
234: }
238: PetscErrorCode MatSetValuesBlocked_MPIBAIJ(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
239: {
240: Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)mat->data;
242: PetscInt i,N = m*n*b->bs2;
243: MatScalar *vsingle;
246: if (N > b->setvalueslen) {
247: PetscFree(b->setvaluescopy);
248: PetscMalloc(N*sizeof(MatScalar),&b->setvaluescopy);
249: b->setvalueslen = N;
250: }
251: vsingle = b->setvaluescopy;
252: for (i=0; i<N; i++) {
253: vsingle[i] = v[i];
254: }
255: MatSetValuesBlocked_MPIBAIJ_MatScalar(mat,m,im,n,in,vsingle,addv);
256: return(0);
257: }
261: PetscErrorCode MatSetValues_MPIBAIJ_HT(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
262: {
263: Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)mat->data;
265: PetscInt i,N = m*n;
266: MatScalar *vsingle;
269: if (N > b->setvalueslen) {
270: PetscFree(b->setvaluescopy);
271: PetscMalloc(N*sizeof(MatScalar),&b->setvaluescopy);
272: b->setvalueslen = N;
273: }
274: vsingle = b->setvaluescopy;
275: for (i=0; i<N; i++) {
276: vsingle[i] = v[i];
277: }
278: MatSetValues_MPIBAIJ_HT_MatScalar(mat,m,im,n,in,vsingle,addv);
279: return(0);
280: }
284: PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv)
285: {
286: Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)mat->data;
288: PetscInt i,N = m*n*b->bs2;
289: MatScalar *vsingle;
292: if (N > b->setvalueslen) {
293: PetscFree(b->setvaluescopy);
294: PetscMalloc(N*sizeof(MatScalar),&b->setvaluescopy);
295: b->setvalueslen = N;
296: }
297: vsingle = b->setvaluescopy;
298: for (i=0; i<N; i++) {
299: vsingle[i] = v[i];
300: }
301: MatSetValuesBlocked_MPIBAIJ_HT_MatScalar(mat,m,im,n,in,vsingle,addv);
302: return(0);
303: }
304: #endif
308: PetscErrorCode MatSetValues_MPIBAIJ_MatScalar(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const MatScalar v[],InsertMode addv)
309: {
310: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
311: MatScalar value;
312: PetscTruth roworiented = baij->roworiented;
314: PetscInt i,j,row,col;
315: PetscInt rstart_orig=mat->rmap.rstart;
316: PetscInt rend_orig=mat->rmap.rend,cstart_orig=mat->cmap.rstart;
317: PetscInt cend_orig=mat->cmap.rend,bs=mat->rmap.bs;
319: /* Some Variables required in the macro */
320: Mat A = baij->A;
321: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)(A)->data;
322: PetscInt *aimax=a->imax,*ai=a->i,*ailen=a->ilen,*aj=a->j;
323: MatScalar *aa=a->a;
325: Mat B = baij->B;
326: Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(B)->data;
327: PetscInt *bimax=b->imax,*bi=b->i,*bilen=b->ilen,*bj=b->j;
328: MatScalar *ba=b->a;
330: PetscInt *rp,ii,nrow,_i,rmax,N,brow,bcol;
331: PetscInt low,high,t,ridx,cidx,bs2=a->bs2;
332: MatScalar *ap,*bap;
335: for (i=0; i<m; i++) {
336: if (im[i] < 0) continue;
337: #if defined(PETSC_USE_DEBUG)
338: if (im[i] >= mat->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap.N-1);
339: #endif
340: if (im[i] >= rstart_orig && im[i] < rend_orig) {
341: row = im[i] - rstart_orig;
342: for (j=0; j<n; j++) {
343: if (in[j] >= cstart_orig && in[j] < cend_orig){
344: col = in[j] - cstart_orig;
345: if (roworiented) value = v[i*n+j]; else value = v[i+j*m];
346: MatSetValues_SeqBAIJ_A_Private(row,col,value,addv);
347: /* MatSetValues_SeqBAIJ(baij->A,1,&row,1,&col,&value,addv); */
348: } else if (in[j] < 0) continue;
349: #if defined(PETSC_USE_DEBUG)
350: else if (in[j] >= mat->cmap.N) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[i],mat->cmap.N-1);}
351: #endif
352: else {
353: if (mat->was_assembled) {
354: if (!baij->colmap) {
355: CreateColmap_MPIBAIJ_Private(mat);
356: }
357: #if defined (PETSC_USE_CTABLE)
358: PetscTableFind(baij->colmap,in[j]/bs + 1,&col);
359: col = col - 1;
360: #else
361: col = baij->colmap[in[j]/bs] - 1;
362: #endif
363: if (col < 0 && !((Mat_SeqBAIJ*)(baij->A->data))->nonew) {
364: DisAssemble_MPIBAIJ(mat);
365: col = in[j];
366: /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */
367: B = baij->B;
368: b = (Mat_SeqBAIJ*)(B)->data;
369: bimax=b->imax;bi=b->i;bilen=b->ilen;bj=b->j;
370: ba=b->a;
371: } else col += in[j]%bs;
372: } else col = in[j];
373: if (roworiented) value = v[i*n+j]; else value = v[i+j*m];
374: MatSetValues_SeqBAIJ_B_Private(row,col,value,addv);
375: /* MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv); */
376: }
377: }
378: } else {
379: if (!baij->donotstash) {
380: if (roworiented) {
381: MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n);
382: } else {
383: MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m);
384: }
385: }
386: }
387: }
388: return(0);
389: }
393: PetscErrorCode MatSetValuesBlocked_MPIBAIJ_MatScalar(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const MatScalar v[],InsertMode addv)
394: {
395: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
396: const MatScalar *value;
397: MatScalar *barray=baij->barray;
398: PetscTruth roworiented = baij->roworiented;
399: PetscErrorCode ierr;
400: PetscInt i,j,ii,jj,row,col,rstart=baij->rstartbs;
401: PetscInt rend=baij->rendbs,cstart=baij->cstartbs,stepval;
402: PetscInt cend=baij->cendbs,bs=mat->rmap.bs,bs2=baij->bs2;
403:
405: if(!barray) {
406: PetscMalloc(bs2*sizeof(MatScalar),&barray);
407: baij->barray = barray;
408: }
410: if (roworiented) {
411: stepval = (n-1)*bs;
412: } else {
413: stepval = (m-1)*bs;
414: }
415: for (i=0; i<m; i++) {
416: if (im[i] < 0) continue;
417: #if defined(PETSC_USE_DEBUG)
418: if (im[i] >= baij->Mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large, row %D max %D",im[i],baij->Mbs-1);
419: #endif
420: if (im[i] >= rstart && im[i] < rend) {
421: row = im[i] - rstart;
422: for (j=0; j<n; j++) {
423: /* If NumCol = 1 then a copy is not required */
424: if ((roworiented) && (n == 1)) {
425: barray = (MatScalar*)v + i*bs2;
426: } else if((!roworiented) && (m == 1)) {
427: barray = (MatScalar*)v + j*bs2;
428: } else { /* Here a copy is required */
429: if (roworiented) {
430: value = v + i*(stepval+bs)*bs + j*bs;
431: } else {
432: value = v + j*(stepval+bs)*bs + i*bs;
433: }
434: for (ii=0; ii<bs; ii++,value+=stepval) {
435: for (jj=0; jj<bs; jj++) {
436: *barray++ = *value++;
437: }
438: }
439: barray -=bs2;
440: }
441:
442: if (in[j] >= cstart && in[j] < cend){
443: col = in[j] - cstart;
444: MatSetValuesBlocked_SeqBAIJ_MatScalar(baij->A,1,&row,1,&col,barray,addv);
445: }
446: else if (in[j] < 0) continue;
447: #if defined(PETSC_USE_DEBUG)
448: else if (in[j] >= baij->Nbs) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large, col %D max %D",in[j],baij->Nbs-1);}
449: #endif
450: else {
451: if (mat->was_assembled) {
452: if (!baij->colmap) {
453: CreateColmap_MPIBAIJ_Private(mat);
454: }
456: #if defined(PETSC_USE_DEBUG)
457: #if defined (PETSC_USE_CTABLE)
458: { PetscInt data;
459: PetscTableFind(baij->colmap,in[j]+1,&data);
460: if ((data - 1) % bs) SETERRQ(PETSC_ERR_PLIB,"Incorrect colmap");
461: }
462: #else
463: if ((baij->colmap[in[j]] - 1) % bs) SETERRQ(PETSC_ERR_PLIB,"Incorrect colmap");
464: #endif
465: #endif
466: #if defined (PETSC_USE_CTABLE)
467: PetscTableFind(baij->colmap,in[j]+1,&col);
468: col = (col - 1)/bs;
469: #else
470: col = (baij->colmap[in[j]] - 1)/bs;
471: #endif
472: if (col < 0 && !((Mat_SeqBAIJ*)(baij->A->data))->nonew) {
473: DisAssemble_MPIBAIJ(mat);
474: col = in[j];
475: }
476: }
477: else col = in[j];
478: MatSetValuesBlocked_SeqBAIJ_MatScalar(baij->B,1,&row,1,&col,barray,addv);
479: }
480: }
481: } else {
482: if (!baij->donotstash) {
483: if (roworiented) {
484: MatStashValuesRowBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);
485: } else {
486: MatStashValuesColBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);
487: }
488: }
489: }
490: }
491: return(0);
492: }
494: #define HASH_KEY 0.6180339887
495: #define HASH(size,key,tmp) (tmp = (key)*HASH_KEY,(PetscInt)((size)*(tmp-(PetscInt)tmp)))
496: /* #define HASH(size,key) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
497: /* #define HASH(size,key,tmp) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
500: PetscErrorCode MatSetValues_MPIBAIJ_HT_MatScalar(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const MatScalar v[],InsertMode addv)
501: {
502: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
503: PetscTruth roworiented = baij->roworiented;
505: PetscInt i,j,row,col;
506: PetscInt rstart_orig=mat->rmap.rstart;
507: PetscInt rend_orig=mat->rmap.rend,Nbs=baij->Nbs;
508: PetscInt h1,key,size=baij->ht_size,bs=mat->rmap.bs,*HT=baij->ht,idx;
509: PetscReal tmp;
510: MatScalar **HD = baij->hd,value;
511: #if defined(PETSC_USE_DEBUG)
512: PetscInt total_ct=baij->ht_total_ct,insert_ct=baij->ht_insert_ct;
513: #endif
517: for (i=0; i<m; i++) {
518: #if defined(PETSC_USE_DEBUG)
519: if (im[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row");
520: if (im[i] >= mat->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap.N-1);
521: #endif
522: row = im[i];
523: if (row >= rstart_orig && row < rend_orig) {
524: for (j=0; j<n; j++) {
525: col = in[j];
526: if (roworiented) value = v[i*n+j]; else value = v[i+j*m];
527: /* Look up PetscInto the Hash Table */
528: key = (row/bs)*Nbs+(col/bs)+1;
529: h1 = HASH(size,key,tmp);
531:
532: idx = h1;
533: #if defined(PETSC_USE_DEBUG)
534: insert_ct++;
535: total_ct++;
536: if (HT[idx] != key) {
537: for (idx=h1; (idx<size) && (HT[idx]!=key); idx++,total_ct++);
538: if (idx == size) {
539: for (idx=0; (idx<h1) && (HT[idx]!=key); idx++,total_ct++);
540: if (idx == h1) {
541: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col);
542: }
543: }
544: }
545: #else
546: if (HT[idx] != key) {
547: for (idx=h1; (idx<size) && (HT[idx]!=key); idx++);
548: if (idx == size) {
549: for (idx=0; (idx<h1) && (HT[idx]!=key); idx++);
550: if (idx == h1) {
551: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col);
552: }
553: }
554: }
555: #endif
556: /* A HASH table entry is found, so insert the values at the correct address */
557: if (addv == ADD_VALUES) *(HD[idx]+ (col % bs)*bs + (row % bs)) += value;
558: else *(HD[idx]+ (col % bs)*bs + (row % bs)) = value;
559: }
560: } else {
561: if (!baij->donotstash) {
562: if (roworiented) {
563: MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n);
564: } else {
565: MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m);
566: }
567: }
568: }
569: }
570: #if defined(PETSC_USE_DEBUG)
571: baij->ht_total_ct = total_ct;
572: baij->ht_insert_ct = insert_ct;
573: #endif
574: return(0);
575: }
579: PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT_MatScalar(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const MatScalar v[],InsertMode addv)
580: {
581: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
582: PetscTruth roworiented = baij->roworiented;
583: PetscErrorCode ierr;
584: PetscInt i,j,ii,jj,row,col;
585: PetscInt rstart=baij->rstartbs;
586: PetscInt rend=mat->rmap.rend,stepval,bs=mat->rmap.bs,bs2=baij->bs2,nbs2=n*bs2;
587: PetscInt h1,key,size=baij->ht_size,idx,*HT=baij->ht,Nbs=baij->Nbs;
588: PetscReal tmp;
589: MatScalar **HD = baij->hd,*baij_a;
590: const MatScalar *v_t,*value;
591: #if defined(PETSC_USE_DEBUG)
592: PetscInt total_ct=baij->ht_total_ct,insert_ct=baij->ht_insert_ct;
593: #endif
594:
597: if (roworiented) {
598: stepval = (n-1)*bs;
599: } else {
600: stepval = (m-1)*bs;
601: }
602: for (i=0; i<m; i++) {
603: #if defined(PETSC_USE_DEBUG)
604: if (im[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",im[i]);
605: if (im[i] >= baij->Mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],baij->Mbs-1);
606: #endif
607: row = im[i];
608: v_t = v + i*nbs2;
609: if (row >= rstart && row < rend) {
610: for (j=0; j<n; j++) {
611: col = in[j];
613: /* Look up into the Hash Table */
614: key = row*Nbs+col+1;
615: h1 = HASH(size,key,tmp);
616:
617: idx = h1;
618: #if defined(PETSC_USE_DEBUG)
619: total_ct++;
620: insert_ct++;
621: if (HT[idx] != key) {
622: for (idx=h1; (idx<size) && (HT[idx]!=key); idx++,total_ct++);
623: if (idx == size) {
624: for (idx=0; (idx<h1) && (HT[idx]!=key); idx++,total_ct++);
625: if (idx == h1) {
626: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col);
627: }
628: }
629: }
630: #else
631: if (HT[idx] != key) {
632: for (idx=h1; (idx<size) && (HT[idx]!=key); idx++);
633: if (idx == size) {
634: for (idx=0; (idx<h1) && (HT[idx]!=key); idx++);
635: if (idx == h1) {
636: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col);
637: }
638: }
639: }
640: #endif
641: baij_a = HD[idx];
642: if (roworiented) {
643: /*value = v + i*(stepval+bs)*bs + j*bs;*/
644: /* value = v + (i*(stepval+bs)+j)*bs; */
645: value = v_t;
646: v_t += bs;
647: if (addv == ADD_VALUES) {
648: for (ii=0; ii<bs; ii++,value+=stepval) {
649: for (jj=ii; jj<bs2; jj+=bs) {
650: baij_a[jj] += *value++;
651: }
652: }
653: } else {
654: for (ii=0; ii<bs; ii++,value+=stepval) {
655: for (jj=ii; jj<bs2; jj+=bs) {
656: baij_a[jj] = *value++;
657: }
658: }
659: }
660: } else {
661: value = v + j*(stepval+bs)*bs + i*bs;
662: if (addv == ADD_VALUES) {
663: for (ii=0; ii<bs; ii++,value+=stepval,baij_a+=bs) {
664: for (jj=0; jj<bs; jj++) {
665: baij_a[jj] += *value++;
666: }
667: }
668: } else {
669: for (ii=0; ii<bs; ii++,value+=stepval,baij_a+=bs) {
670: for (jj=0; jj<bs; jj++) {
671: baij_a[jj] = *value++;
672: }
673: }
674: }
675: }
676: }
677: } else {
678: if (!baij->donotstash) {
679: if (roworiented) {
680: MatStashValuesRowBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);
681: } else {
682: MatStashValuesColBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);
683: }
684: }
685: }
686: }
687: #if defined(PETSC_USE_DEBUG)
688: baij->ht_total_ct = total_ct;
689: baij->ht_insert_ct = insert_ct;
690: #endif
691: return(0);
692: }
696: PetscErrorCode MatGetValues_MPIBAIJ(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
697: {
698: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
700: PetscInt bs=mat->rmap.bs,i,j,bsrstart = mat->rmap.rstart,bsrend = mat->rmap.rend;
701: PetscInt bscstart = mat->cmap.rstart,bscend = mat->cmap.rend,row,col,data;
704: for (i=0; i<m; i++) {
705: if (idxm[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);
706: if (idxm[i] >= mat->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm[i],mat->rmap.N-1);
707: if (idxm[i] >= bsrstart && idxm[i] < bsrend) {
708: row = idxm[i] - bsrstart;
709: for (j=0; j<n; j++) {
710: if (idxn[j] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]);
711: if (idxn[j] >= mat->cmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",idxn[j],mat->cmap.N-1);
712: if (idxn[j] >= bscstart && idxn[j] < bscend){
713: col = idxn[j] - bscstart;
714: MatGetValues_SeqBAIJ(baij->A,1,&row,1,&col,v+i*n+j);
715: } else {
716: if (!baij->colmap) {
717: CreateColmap_MPIBAIJ_Private(mat);
718: }
719: #if defined (PETSC_USE_CTABLE)
720: PetscTableFind(baij->colmap,idxn[j]/bs+1,&data);
721: data --;
722: #else
723: data = baij->colmap[idxn[j]/bs]-1;
724: #endif
725: if((data < 0) || (baij->garray[data/bs] != idxn[j]/bs)) *(v+i*n+j) = 0.0;
726: else {
727: col = data + idxn[j]%bs;
728: MatGetValues_SeqBAIJ(baij->B,1,&row,1,&col,v+i*n+j);
729: }
730: }
731: }
732: } else {
733: SETERRQ(PETSC_ERR_SUP,"Only local values currently supported");
734: }
735: }
736: return(0);
737: }
741: PetscErrorCode MatNorm_MPIBAIJ(Mat mat,NormType type,PetscReal *nrm)
742: {
743: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
744: Mat_SeqBAIJ *amat = (Mat_SeqBAIJ*)baij->A->data,*bmat = (Mat_SeqBAIJ*)baij->B->data;
746: PetscInt i,j,bs2=baij->bs2,bs=baij->A->rmap.bs,nz,row,col;
747: PetscReal sum = 0.0;
748: MatScalar *v;
751: if (baij->size == 1) {
752: MatNorm(baij->A,type,nrm);
753: } else {
754: if (type == NORM_FROBENIUS) {
755: v = amat->a;
756: nz = amat->nz*bs2;
757: for (i=0; i<nz; i++) {
758: #if defined(PETSC_USE_COMPLEX)
759: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
760: #else
761: sum += (*v)*(*v); v++;
762: #endif
763: }
764: v = bmat->a;
765: nz = bmat->nz*bs2;
766: for (i=0; i<nz; i++) {
767: #if defined(PETSC_USE_COMPLEX)
768: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
769: #else
770: sum += (*v)*(*v); v++;
771: #endif
772: }
773: MPI_Allreduce(&sum,nrm,1,MPIU_REAL,MPI_SUM,mat->comm);
774: *nrm = sqrt(*nrm);
775: } else if (type == NORM_1) { /* max column sum */
776: PetscReal *tmp,*tmp2;
777: PetscInt *jj,*garray=baij->garray,cstart=baij->rstartbs;
778: PetscMalloc((2*mat->cmap.N+1)*sizeof(PetscReal),&tmp);
779: tmp2 = tmp + mat->cmap.N;
780: PetscMemzero(tmp,mat->cmap.N*sizeof(PetscReal));
781: v = amat->a; jj = amat->j;
782: for (i=0; i<amat->nz; i++) {
783: for (j=0; j<bs; j++){
784: col = bs*(cstart + *jj) + j; /* column index */
785: for (row=0; row<bs; row++){
786: tmp[col] += PetscAbsScalar(*v); v++;
787: }
788: }
789: jj++;
790: }
791: v = bmat->a; jj = bmat->j;
792: for (i=0; i<bmat->nz; i++) {
793: for (j=0; j<bs; j++){
794: col = bs*garray[*jj] + j;
795: for (row=0; row<bs; row++){
796: tmp[col] += PetscAbsScalar(*v); v++;
797: }
798: }
799: jj++;
800: }
801: MPI_Allreduce(tmp,tmp2,mat->cmap.N,MPIU_REAL,MPI_SUM,mat->comm);
802: *nrm = 0.0;
803: for (j=0; j<mat->cmap.N; j++) {
804: if (tmp2[j] > *nrm) *nrm = tmp2[j];
805: }
806: PetscFree(tmp);
807: } else if (type == NORM_INFINITY) { /* max row sum */
808: PetscReal *sums;
809: PetscMalloc(bs*sizeof(PetscReal),&sums);CHKERRQ(ierr)
810: sum = 0.0;
811: for (j=0; j<amat->mbs; j++) {
812: for (row=0; row<bs; row++) sums[row] = 0.0;
813: v = amat->a + bs2*amat->i[j];
814: nz = amat->i[j+1]-amat->i[j];
815: for (i=0; i<nz; i++) {
816: for (col=0; col<bs; col++){
817: for (row=0; row<bs; row++){
818: sums[row] += PetscAbsScalar(*v); v++;
819: }
820: }
821: }
822: v = bmat->a + bs2*bmat->i[j];
823: nz = bmat->i[j+1]-bmat->i[j];
824: for (i=0; i<nz; i++) {
825: for (col=0; col<bs; col++){
826: for (row=0; row<bs; row++){
827: sums[row] += PetscAbsScalar(*v); v++;
828: }
829: }
830: }
831: for (row=0; row<bs; row++){
832: if (sums[row] > sum) sum = sums[row];
833: }
834: }
835: MPI_Allreduce(&sum,nrm,1,MPIU_REAL,MPI_MAX,mat->comm);
836: PetscFree(sums);
837: } else {
838: SETERRQ(PETSC_ERR_SUP,"No support for this norm yet");
839: }
840: }
841: return(0);
842: }
844: /*
845: Creates the hash table, and sets the table
846: This table is created only once.
847: If new entried need to be added to the matrix
848: then the hash table has to be destroyed and
849: recreated.
850: */
853: PetscErrorCode MatCreateHashTable_MPIBAIJ_Private(Mat mat,PetscReal factor)
854: {
855: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
856: Mat A = baij->A,B=baij->B;
857: Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data,*b=(Mat_SeqBAIJ *)B->data;
858: PetscInt i,j,k,nz=a->nz+b->nz,h1,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
860: PetscInt size,bs2=baij->bs2,rstart=baij->rstartbs;
861: PetscInt cstart=baij->cstartbs,*garray=baij->garray,row,col,Nbs=baij->Nbs;
862: PetscInt *HT,key;
863: MatScalar **HD;
864: PetscReal tmp;
865: #if defined(PETSC_USE_INFO)
866: PetscInt ct=0,max=0;
867: #endif
870: baij->ht_size=(PetscInt)(factor*nz);
871: size = baij->ht_size;
873: if (baij->ht) {
874: return(0);
875: }
876:
877: /* Allocate Memory for Hash Table */
878: PetscMalloc((size)*(sizeof(PetscInt)+sizeof(MatScalar*))+1,&baij->hd);
879: baij->ht = (PetscInt*)(baij->hd + size);
880: HD = baij->hd;
881: HT = baij->ht;
884: PetscMemzero(HD,size*(sizeof(PetscInt)+sizeof(PetscScalar*)));
885:
887: /* Loop Over A */
888: for (i=0; i<a->mbs; i++) {
889: for (j=ai[i]; j<ai[i+1]; j++) {
890: row = i+rstart;
891: col = aj[j]+cstart;
892:
893: key = row*Nbs + col + 1;
894: h1 = HASH(size,key,tmp);
895: for (k=0; k<size; k++){
896: if (!HT[(h1+k)%size]) {
897: HT[(h1+k)%size] = key;
898: HD[(h1+k)%size] = a->a + j*bs2;
899: break;
900: #if defined(PETSC_USE_INFO)
901: } else {
902: ct++;
903: #endif
904: }
905: }
906: #if defined(PETSC_USE_INFO)
907: if (k> max) max = k;
908: #endif
909: }
910: }
911: /* Loop Over B */
912: for (i=0; i<b->mbs; i++) {
913: for (j=bi[i]; j<bi[i+1]; j++) {
914: row = i+rstart;
915: col = garray[bj[j]];
916: key = row*Nbs + col + 1;
917: h1 = HASH(size,key,tmp);
918: for (k=0; k<size; k++){
919: if (!HT[(h1+k)%size]) {
920: HT[(h1+k)%size] = key;
921: HD[(h1+k)%size] = b->a + j*bs2;
922: break;
923: #if defined(PETSC_USE_INFO)
924: } else {
925: ct++;
926: #endif
927: }
928: }
929: #if defined(PETSC_USE_INFO)
930: if (k> max) max = k;
931: #endif
932: }
933: }
934:
935: /* Print Summary */
936: #if defined(PETSC_USE_INFO)
937: for (i=0,j=0; i<size; i++) {
938: if (HT[i]) {j++;}
939: }
940: PetscInfo2(0,"Average Search = %5.2f,max search = %D\n",(!j)? 0.0:((PetscReal)(ct+j))/j,max);
941: #endif
942: return(0);
943: }
947: PetscErrorCode MatAssemblyBegin_MPIBAIJ(Mat mat,MatAssemblyType mode)
948: {
949: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
951: PetscInt nstash,reallocs;
952: InsertMode addv;
955: if (baij->donotstash) {
956: return(0);
957: }
959: /* make sure all processors are either in INSERTMODE or ADDMODE */
960: MPI_Allreduce(&mat->insertmode,&addv,1,MPI_INT,MPI_BOR,mat->comm);
961: if (addv == (ADD_VALUES|INSERT_VALUES)) {
962: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Some processors inserted others added");
963: }
964: mat->insertmode = addv; /* in case this processor had no cache */
966: MatStashScatterBegin_Private(&mat->stash,mat->rmap.range);
967: MatStashScatterBegin_Private(&mat->bstash,baij->rangebs);
968: MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);
969: PetscInfo2(0,"Stash has %D entries,uses %D mallocs.\n",nstash,reallocs);
970: MatStashGetInfo_Private(&mat->bstash,&nstash,&reallocs);
971: PetscInfo2(0,"Block-Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);
972: return(0);
973: }
977: PetscErrorCode MatAssemblyEnd_MPIBAIJ(Mat mat,MatAssemblyType mode)
978: {
979: Mat_MPIBAIJ *baij=(Mat_MPIBAIJ*)mat->data;
980: Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)baij->A->data;
982: PetscInt i,j,rstart,ncols,flg,bs2=baij->bs2;
983: PetscInt *row,*col,other_disassembled;
984: PetscTruth r1,r2,r3;
985: MatScalar *val;
986: InsertMode addv = mat->insertmode;
987: PetscMPIInt n;
989: /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */
991: if (!baij->donotstash) {
992: while (1) {
993: MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);
994: if (!flg) break;
996: for (i=0; i<n;) {
997: /* Now identify the consecutive vals belonging to the same row */
998: for (j=i,rstart=row[j]; j<n; j++) { if (row[j] != rstart) break; }
999: if (j < n) ncols = j-i;
1000: else ncols = n-i;
1001: /* Now assemble all these values with a single function call */
1002: MatSetValues_MPIBAIJ_MatScalar(mat,1,row+i,ncols,col+i,val+i,addv);
1003: i = j;
1004: }
1005: }
1006: MatStashScatterEnd_Private(&mat->stash);
1007: /* Now process the block-stash. Since the values are stashed column-oriented,
1008: set the roworiented flag to column oriented, and after MatSetValues()
1009: restore the original flags */
1010: r1 = baij->roworiented;
1011: r2 = a->roworiented;
1012: r3 = ((Mat_SeqBAIJ*)baij->B->data)->roworiented;
1013: baij->roworiented = PETSC_FALSE;
1014: a->roworiented = PETSC_FALSE;
1015: (((Mat_SeqBAIJ*)baij->B->data))->roworiented = PETSC_FALSE; /* b->roworiented */
1016: while (1) {
1017: MatStashScatterGetMesg_Private(&mat->bstash,&n,&row,&col,&val,&flg);
1018: if (!flg) break;
1019:
1020: for (i=0; i<n;) {
1021: /* Now identify the consecutive vals belonging to the same row */
1022: for (j=i,rstart=row[j]; j<n; j++) { if (row[j] != rstart) break; }
1023: if (j < n) ncols = j-i;
1024: else ncols = n-i;
1025: MatSetValuesBlocked_MPIBAIJ_MatScalar(mat,1,row+i,ncols,col+i,val+i*bs2,addv);
1026: i = j;
1027: }
1028: }
1029: MatStashScatterEnd_Private(&mat->bstash);
1030: baij->roworiented = r1;
1031: a->roworiented = r2;
1032: ((Mat_SeqBAIJ*)baij->B->data)->roworiented = r3; /* b->roworiented */
1033: }
1034:
1035: MatAssemblyBegin(baij->A,mode);
1036: MatAssemblyEnd(baij->A,mode);
1038: /* determine if any processor has disassembled, if so we must
1039: also disassemble ourselfs, in order that we may reassemble. */
1040: /*
1041: if nonzero structure of submatrix B cannot change then we know that
1042: no processor disassembled thus we can skip this stuff
1043: */
1044: if (!((Mat_SeqBAIJ*)baij->B->data)->nonew) {
1045: MPI_Allreduce(&mat->was_assembled,&other_disassembled,1,MPI_INT,MPI_PROD,mat->comm);
1046: if (mat->was_assembled && !other_disassembled) {
1047: DisAssemble_MPIBAIJ(mat);
1048: }
1049: }
1051: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
1052: MatSetUpMultiply_MPIBAIJ(mat);
1053: }
1054: ((Mat_SeqBAIJ*)baij->B->data)->compressedrow.use = PETSC_TRUE; /* b->compressedrow.use */
1055: MatAssemblyBegin(baij->B,mode);
1056: MatAssemblyEnd(baij->B,mode);
1057:
1058: #if defined(PETSC_USE_INFO)
1059: if (baij->ht && mode== MAT_FINAL_ASSEMBLY) {
1060: PetscInfo1(0,"Average Hash Table Search in MatSetValues = %5.2f\n",((PetscReal)baij->ht_total_ct)/baij->ht_insert_ct);
1061: baij->ht_total_ct = 0;
1062: baij->ht_insert_ct = 0;
1063: }
1064: #endif
1065: if (baij->ht_flag && !baij->ht && mode == MAT_FINAL_ASSEMBLY) {
1066: MatCreateHashTable_MPIBAIJ_Private(mat,baij->ht_fact);
1067: mat->ops->setvalues = MatSetValues_MPIBAIJ_HT;
1068: mat->ops->setvaluesblocked = MatSetValuesBlocked_MPIBAIJ_HT;
1069: }
1071: PetscFree(baij->rowvalues);
1072: baij->rowvalues = 0;
1073: return(0);
1074: }
1078: static PetscErrorCode MatView_MPIBAIJ_ASCIIorDraworSocket(Mat mat,PetscViewer viewer)
1079: {
1080: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
1081: PetscErrorCode ierr;
1082: PetscMPIInt size = baij->size,rank = baij->rank;
1083: PetscInt bs = mat->rmap.bs;
1084: PetscTruth iascii,isdraw;
1085: PetscViewer sviewer;
1086: PetscViewerFormat format;
1089: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1090: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
1091: if (iascii) {
1092: PetscViewerGetFormat(viewer,&format);
1093: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1094: MatInfo info;
1095: MPI_Comm_rank(mat->comm,&rank);
1096: MatGetInfo(mat,MAT_LOCAL,&info);
1097: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D bs %D mem %D\n",
1098: rank,mat->rmap.N,(PetscInt)info.nz_used*bs,(PetscInt)info.nz_allocated*bs,
1099: mat->rmap.bs,(PetscInt)info.memory);
1100: MatGetInfo(baij->A,MAT_LOCAL,&info);
1101: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used*bs);
1102: MatGetInfo(baij->B,MAT_LOCAL,&info);
1103: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used*bs);
1104: PetscViewerFlush(viewer);
1105: PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");
1106: VecScatterView(baij->Mvctx,viewer);
1107: return(0);
1108: } else if (format == PETSC_VIEWER_ASCII_INFO) {
1109: PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);
1110: return(0);
1111: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
1112: return(0);
1113: }
1114: }
1116: if (isdraw) {
1117: PetscDraw draw;
1118: PetscTruth isnull;
1119: PetscViewerDrawGetDraw(viewer,0,&draw);
1120: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
1121: }
1123: if (size == 1) {
1124: PetscObjectSetName((PetscObject)baij->A,mat->name);
1125: MatView(baij->A,viewer);
1126: } else {
1127: /* assemble the entire matrix onto first processor. */
1128: Mat A;
1129: Mat_SeqBAIJ *Aloc;
1130: PetscInt M = mat->rmap.N,N = mat->cmap.N,*ai,*aj,col,i,j,k,*rvals,mbs = baij->mbs;
1131: MatScalar *a;
1133: /* Here we are creating a temporary matrix, so will assume MPIBAIJ is acceptable */
1134: /* Perhaps this should be the type of mat? */
1135: MatCreate(mat->comm,&A);
1136: if (!rank) {
1137: MatSetSizes(A,M,N,M,N);
1138: } else {
1139: MatSetSizes(A,0,0,M,N);
1140: }
1141: MatSetType(A,MATMPIBAIJ);
1142: MatMPIBAIJSetPreallocation(A,mat->rmap.bs,0,PETSC_NULL,0,PETSC_NULL);
1143: PetscLogObjectParent(mat,A);
1145: /* copy over the A part */
1146: Aloc = (Mat_SeqBAIJ*)baij->A->data;
1147: ai = Aloc->i; aj = Aloc->j; a = Aloc->a;
1148: PetscMalloc(bs*sizeof(PetscInt),&rvals);
1150: for (i=0; i<mbs; i++) {
1151: rvals[0] = bs*(baij->rstartbs + i);
1152: for (j=1; j<bs; j++) { rvals[j] = rvals[j-1] + 1; }
1153: for (j=ai[i]; j<ai[i+1]; j++) {
1154: col = (baij->cstartbs+aj[j])*bs;
1155: for (k=0; k<bs; k++) {
1156: MatSetValues_MPIBAIJ_MatScalar(A,bs,rvals,1,&col,a,INSERT_VALUES);
1157: col++; a += bs;
1158: }
1159: }
1160: }
1161: /* copy over the B part */
1162: Aloc = (Mat_SeqBAIJ*)baij->B->data;
1163: ai = Aloc->i; aj = Aloc->j; a = Aloc->a;
1164: for (i=0; i<mbs; i++) {
1165: rvals[0] = bs*(baij->rstartbs + i);
1166: for (j=1; j<bs; j++) { rvals[j] = rvals[j-1] + 1; }
1167: for (j=ai[i]; j<ai[i+1]; j++) {
1168: col = baij->garray[aj[j]]*bs;
1169: for (k=0; k<bs; k++) {
1170: MatSetValues_MPIBAIJ_MatScalar(A,bs,rvals,1,&col,a,INSERT_VALUES);
1171: col++; a += bs;
1172: }
1173: }
1174: }
1175: PetscFree(rvals);
1176: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1177: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1178: /*
1179: Everyone has to call to draw the matrix since the graphics waits are
1180: synchronized across all processors that share the PetscDraw object
1181: */
1182: PetscViewerGetSingleton(viewer,&sviewer);
1183: if (!rank) {
1184: PetscObjectSetName((PetscObject)((Mat_MPIBAIJ*)(A->data))->A,mat->name);
1185: MatView(((Mat_MPIBAIJ*)(A->data))->A,sviewer);
1186: }
1187: PetscViewerRestoreSingleton(viewer,&sviewer);
1188: MatDestroy(A);
1189: }
1190: return(0);
1191: }
1195: PetscErrorCode MatView_MPIBAIJ(Mat mat,PetscViewer viewer)
1196: {
1198: PetscTruth iascii,isdraw,issocket,isbinary;
1201: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1202: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
1203: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);
1204: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
1205: if (iascii || isdraw || issocket || isbinary) {
1206: MatView_MPIBAIJ_ASCIIorDraworSocket(mat,viewer);
1207: } else {
1208: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by MPIBAIJ matrices",((PetscObject)viewer)->type_name);
1209: }
1210: return(0);
1211: }
1215: PetscErrorCode MatDestroy_MPIBAIJ(Mat mat)
1216: {
1217: Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data;
1221: #if defined(PETSC_USE_LOG)
1222: PetscLogObjectState((PetscObject)mat,"Rows=%D,Cols=%D",mat->rmap.N,mat->cmap.N);
1223: #endif
1224: MatStashDestroy_Private(&mat->stash);
1225: MatStashDestroy_Private(&mat->bstash);
1226: MatDestroy(baij->A);
1227: MatDestroy(baij->B);
1228: #if defined (PETSC_USE_CTABLE)
1229: if (baij->colmap) {PetscTableDestroy(baij->colmap);}
1230: #else
1231: PetscFree(baij->colmap);
1232: #endif
1233: PetscFree(baij->garray);
1234: if (baij->lvec) {VecDestroy(baij->lvec);}
1235: if (baij->Mvctx) {VecScatterDestroy(baij->Mvctx);}
1236: PetscFree(baij->rowvalues);
1237: PetscFree(baij->barray);
1238: PetscFree(baij->hd);
1239: #if defined(PETSC_USE_MAT_SINGLE)
1240: PetscFree(baij->setvaluescopy);
1241: #endif
1242: PetscFree(baij->rangebs);
1243: PetscFree(baij);
1245: PetscObjectChangeTypeName((PetscObject)mat,0);
1246: PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C","",PETSC_NULL);
1247: PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C","",PETSC_NULL);
1248: PetscObjectComposeFunction((