Actual source code: blockmat.c
petsc-3.12.5 2020-03-29
2: /*
3: This provides a matrix that consists of Mats
4: */
6: #include <petsc/private/matimpl.h>
7: #include <../src/mat/impls/baij/seq/baij.h> /* use the common AIJ data-structure */
9: typedef struct {
10: SEQAIJHEADER(Mat);
11: SEQBAIJHEADER;
12: Mat *diags;
14: Vec left,right,middle,workb; /* dummy vectors to perform local parts of product */
15: } Mat_BlockMat;
17: static PetscErrorCode MatSOR_BlockMat_Symmetric(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
18: {
19: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
20: PetscScalar *x;
21: const Mat *v;
22: const PetscScalar *b;
23: PetscErrorCode ierr;
24: PetscInt n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs;
25: const PetscInt *idx;
26: IS row,col;
27: MatFactorInfo info;
28: Vec left = a->left,right = a->right, middle = a->middle;
29: Mat *diag;
32: its = its*lits;
33: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
34: if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
35: if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0");
36: if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift");
37: if ((flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) && !(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) {
38: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot do backward sweep without forward sweep");
39: }
41: if (!a->diags) {
42: PetscMalloc1(mbs,&a->diags);
43: MatFactorInfoInitialize(&info);
44: for (i=0; i<mbs; i++) {
45: MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);
46: MatCholeskyFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,&info);
47: MatCholeskyFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);
48: ISDestroy(&row);
49: ISDestroy(&col);
50: }
51: VecDuplicate(bb,&a->workb);
52: }
53: diag = a->diags;
55: VecSet(xx,0.0);
56: VecGetArray(xx,&x);
57: /* copy right hand side because it must be modified during iteration */
58: VecCopy(bb,a->workb);
59: VecGetArrayRead(a->workb,&b);
61: /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */
62: while (its--) {
63: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
65: for (i=0; i<mbs; i++) {
66: n = a->i[i+1] - a->i[i] - 1;
67: idx = a->j + a->i[i] + 1;
68: v = a->a + a->i[i] + 1;
70: VecSet(left,0.0);
71: for (j=0; j<n; j++) {
72: VecPlaceArray(right,x + idx[j]*bs);
73: MatMultAdd(v[j],right,left,left);
74: VecResetArray(right);
75: }
76: VecPlaceArray(right,b + i*bs);
77: VecAYPX(left,-1.0,right);
78: VecResetArray(right);
80: VecPlaceArray(right,x + i*bs);
81: MatSolve(diag[i],left,right);
83: /* now adjust right hand side, see MatSOR_SeqSBAIJ */
84: for (j=0; j<n; j++) {
85: MatMultTranspose(v[j],right,left);
86: VecPlaceArray(middle,b + idx[j]*bs);
87: VecAXPY(middle,-1.0,left);
88: VecResetArray(middle);
89: }
90: VecResetArray(right);
92: }
93: }
94: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
96: for (i=mbs-1; i>=0; i--) {
97: n = a->i[i+1] - a->i[i] - 1;
98: idx = a->j + a->i[i] + 1;
99: v = a->a + a->i[i] + 1;
101: VecSet(left,0.0);
102: for (j=0; j<n; j++) {
103: VecPlaceArray(right,x + idx[j]*bs);
104: MatMultAdd(v[j],right,left,left);
105: VecResetArray(right);
106: }
107: VecPlaceArray(right,b + i*bs);
108: VecAYPX(left,-1.0,right);
109: VecResetArray(right);
111: VecPlaceArray(right,x + i*bs);
112: MatSolve(diag[i],left,right);
113: VecResetArray(right);
115: }
116: }
117: }
118: VecRestoreArray(xx,&x);
119: VecRestoreArrayRead(a->workb,&b);
120: return(0);
121: }
123: static PetscErrorCode MatSOR_BlockMat(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
124: {
125: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
126: PetscScalar *x;
127: const Mat *v;
128: const PetscScalar *b;
129: PetscErrorCode ierr;
130: PetscInt n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs;
131: const PetscInt *idx;
132: IS row,col;
133: MatFactorInfo info;
134: Vec left = a->left,right = a->right;
135: Mat *diag;
138: its = its*lits;
139: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
140: if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
141: if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0");
142: if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift");
144: if (!a->diags) {
145: PetscMalloc1(mbs,&a->diags);
146: MatFactorInfoInitialize(&info);
147: for (i=0; i<mbs; i++) {
148: MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);
149: MatLUFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,col,&info);
150: MatLUFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);
151: ISDestroy(&row);
152: ISDestroy(&col);
153: }
154: }
155: diag = a->diags;
157: VecSet(xx,0.0);
158: VecGetArray(xx,&x);
159: VecGetArrayRead(bb,&b);
161: /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */
162: while (its--) {
163: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
165: for (i=0; i<mbs; i++) {
166: n = a->i[i+1] - a->i[i];
167: idx = a->j + a->i[i];
168: v = a->a + a->i[i];
170: VecSet(left,0.0);
171: for (j=0; j<n; j++) {
172: if (idx[j] != i) {
173: VecPlaceArray(right,x + idx[j]*bs);
174: MatMultAdd(v[j],right,left,left);
175: VecResetArray(right);
176: }
177: }
178: VecPlaceArray(right,b + i*bs);
179: VecAYPX(left,-1.0,right);
180: VecResetArray(right);
182: VecPlaceArray(right,x + i*bs);
183: MatSolve(diag[i],left,right);
184: VecResetArray(right);
185: }
186: }
187: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
189: for (i=mbs-1; i>=0; i--) {
190: n = a->i[i+1] - a->i[i];
191: idx = a->j + a->i[i];
192: v = a->a + a->i[i];
194: VecSet(left,0.0);
195: for (j=0; j<n; j++) {
196: if (idx[j] != i) {
197: VecPlaceArray(right,x + idx[j]*bs);
198: MatMultAdd(v[j],right,left,left);
199: VecResetArray(right);
200: }
201: }
202: VecPlaceArray(right,b + i*bs);
203: VecAYPX(left,-1.0,right);
204: VecResetArray(right);
206: VecPlaceArray(right,x + i*bs);
207: MatSolve(diag[i],left,right);
208: VecResetArray(right);
210: }
211: }
212: }
213: VecRestoreArray(xx,&x);
214: VecRestoreArrayRead(bb,&b);
215: return(0);
216: }
218: static PetscErrorCode MatSetValues_BlockMat(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
219: {
220: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
221: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
222: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen;
223: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
225: PetscInt ridx,cidx;
226: PetscBool roworiented=a->roworiented;
227: MatScalar value;
228: Mat *ap,*aa = a->a;
231: for (k=0; k<m; k++) { /* loop over added rows */
232: row = im[k];
233: brow = row/bs;
234: if (row < 0) continue;
235: #if defined(PETSC_USE_DEBUG)
236: if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1);
237: #endif
238: rp = aj + ai[brow];
239: ap = aa + ai[brow];
240: rmax = imax[brow];
241: nrow = ailen[brow];
242: low = 0;
243: high = nrow;
244: for (l=0; l<n; l++) { /* loop over added columns */
245: if (in[l] < 0) continue;
246: #if defined(PETSC_USE_DEBUG)
247: if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1);
248: #endif
249: col = in[l]; bcol = col/bs;
250: if (A->symmetric && brow > bcol) continue;
251: ridx = row % bs; cidx = col % bs;
252: if (roworiented) value = v[l + k*n];
253: else value = v[k + l*m];
255: if (col <= lastcol) low = 0;
256: else high = nrow;
257: lastcol = col;
258: while (high-low > 7) {
259: t = (low+high)/2;
260: if (rp[t] > bcol) high = t;
261: else low = t;
262: }
263: for (i=low; i<high; i++) {
264: if (rp[i] > bcol) break;
265: if (rp[i] == bcol) goto noinsert1;
266: }
267: if (nonew == 1) goto noinsert1;
268: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
269: MatSeqXAIJReallocateAIJ(A,a->mbs,1,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,Mat);
270: N = nrow++ - 1; high++;
271: /* shift up all the later entries in this row */
272: for (ii=N; ii>=i; ii--) {
273: rp[ii+1] = rp[ii];
274: ap[ii+1] = ap[ii];
275: }
276: if (N>=i) ap[i] = 0;
277: rp[i] = bcol;
278: a->nz++;
279: A->nonzerostate++;
280: noinsert1:;
281: if (!*(ap+i)) {
282: MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,0,ap+i);
283: }
284: MatSetValues(ap[i],1,&ridx,1,&cidx,&value,is);
285: low = i;
286: }
287: ailen[brow] = nrow;
288: }
289: return(0);
290: }
292: static PetscErrorCode MatLoad_BlockMat(Mat newmat, PetscViewer viewer)
293: {
294: PetscErrorCode ierr;
295: Mat tmpA;
296: PetscInt i,j,m,n,bs = 1,ncols,*lens,currentcol,mbs,**ii,*ilens,nextcol,*llens,cnt = 0;
297: const PetscInt *cols;
298: const PetscScalar *values;
299: PetscBool flg = PETSC_FALSE,notdone;
300: Mat_SeqAIJ *a;
301: Mat_BlockMat *amat;
304: /* force binary viewer to load .info file if it has not yet done so */
305: PetscViewerSetUp(viewer);
306: MatCreate(PETSC_COMM_SELF,&tmpA);
307: MatSetType(tmpA,MATSEQAIJ);
308: MatLoad_SeqAIJ(tmpA,viewer);
310: MatGetLocalSize(tmpA,&m,&n);
311: PetscOptionsBegin(PETSC_COMM_SELF,NULL,"Options for loading BlockMat matrix 1","Mat");
312: PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);
313: PetscOptionsBool("-matload_symmetric","Store the matrix as symmetric","MatLoad",flg,&flg,NULL);
314: PetscOptionsEnd();
316: /* Determine number of nonzero blocks for each block row */
317: a = (Mat_SeqAIJ*) tmpA->data;
318: mbs = m/bs;
319: PetscMalloc3(mbs,&lens,bs,&ii,bs,&ilens);
320: PetscArrayzero(lens,mbs);
322: for (i=0; i<mbs; i++) {
323: for (j=0; j<bs; j++) {
324: ii[j] = a->j + a->i[i*bs + j];
325: ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j];
326: }
328: currentcol = -1;
329: while (PETSC_TRUE) {
330: notdone = PETSC_FALSE;
331: nextcol = 1000000000;
332: for (j=0; j<bs; j++) {
333: while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) {
334: ii[j]++;
335: ilens[j]--;
336: }
337: if (ilens[j] > 0) {
338: notdone = PETSC_TRUE;
339: nextcol = PetscMin(nextcol,ii[j][0]/bs);
340: }
341: }
342: if (!notdone) break;
343: if (!flg || (nextcol >= i)) lens[i]++;
344: currentcol = nextcol;
345: }
346: }
348: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
349: MatSetSizes(newmat,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
350: }
351: MatBlockMatSetPreallocation(newmat,bs,0,lens);
352: if (flg) {
353: MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);
354: }
355: amat = (Mat_BlockMat*)(newmat)->data;
357: /* preallocate the submatrices */
358: PetscMalloc1(bs,&llens);
359: for (i=0; i<mbs; i++) { /* loops for block rows */
360: for (j=0; j<bs; j++) {
361: ii[j] = a->j + a->i[i*bs + j];
362: ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j];
363: }
365: currentcol = 1000000000;
366: for (j=0; j<bs; j++) { /* loop over rows in block finding first nonzero block */
367: if (ilens[j] > 0) {
368: currentcol = PetscMin(currentcol,ii[j][0]/bs);
369: }
370: }
372: while (PETSC_TRUE) { /* loops over blocks in block row */
373: notdone = PETSC_FALSE;
374: nextcol = 1000000000;
375: PetscArrayzero(llens,bs);
376: for (j=0; j<bs; j++) { /* loop over rows in block */
377: while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { /* loop over columns in row */
378: ii[j]++;
379: ilens[j]--;
380: llens[j]++;
381: }
382: if (ilens[j] > 0) {
383: notdone = PETSC_TRUE;
384: nextcol = PetscMin(nextcol,ii[j][0]/bs);
385: }
386: }
387: if (cnt >= amat->maxnz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of blocks found greater than expected %D",cnt);
388: if (!flg || currentcol >= i) {
389: amat->j[cnt] = currentcol;
390: MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,llens,amat->a+cnt++);
391: }
393: if (!notdone) break;
394: currentcol = nextcol;
395: }
396: amat->ilen[i] = lens[i];
397: }
399: PetscFree3(lens,ii,ilens);
400: PetscFree(llens);
402: /* copy over the matrix, one row at a time */
403: for (i=0; i<m; i++) {
404: MatGetRow(tmpA,i,&ncols,&cols,&values);
405: MatSetValues(newmat,1,&i,ncols,cols,values,INSERT_VALUES);
406: MatRestoreRow(tmpA,i,&ncols,&cols,&values);
407: }
408: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
409: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
410: return(0);
411: }
413: static PetscErrorCode MatView_BlockMat(Mat A,PetscViewer viewer)
414: {
415: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
416: PetscErrorCode ierr;
417: const char *name;
418: PetscViewerFormat format;
421: PetscObjectGetName((PetscObject)A,&name);
422: PetscViewerGetFormat(viewer,&format);
423: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
424: PetscViewerASCIIPrintf(viewer,"Nonzero block matrices = %D \n",a->nz);
425: if (A->symmetric) {
426: PetscViewerASCIIPrintf(viewer,"Only upper triangular part of symmetric matrix is stored\n");
427: }
428: }
429: return(0);
430: }
432: static PetscErrorCode MatDestroy_BlockMat(Mat mat)
433: {
435: Mat_BlockMat *bmat = (Mat_BlockMat*)mat->data;
436: PetscInt i;
439: VecDestroy(&bmat->right);
440: VecDestroy(&bmat->left);
441: VecDestroy(&bmat->middle);
442: VecDestroy(&bmat->workb);
443: if (bmat->diags) {
444: for (i=0; i<mat->rmap->n/mat->rmap->bs; i++) {
445: MatDestroy(&bmat->diags[i]);
446: }
447: }
448: if (bmat->a) {
449: for (i=0; i<bmat->nz; i++) {
450: MatDestroy(&bmat->a[i]);
451: }
452: }
453: MatSeqXAIJFreeAIJ(mat,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);
454: PetscFree(mat->data);
455: return(0);
456: }
458: static PetscErrorCode MatMult_BlockMat(Mat A,Vec x,Vec y)
459: {
460: Mat_BlockMat *bmat = (Mat_BlockMat*)A->data;
462: PetscScalar *xx,*yy;
463: PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j;
464: Mat *aa;
467: /*
468: Standard CSR multiply except each entry is a Mat
469: */
470: VecGetArray(x,&xx);
472: VecSet(y,0.0);
473: VecGetArray(y,&yy);
474: aj = bmat->j;
475: aa = bmat->a;
476: ii = bmat->i;
477: for (i=0; i<m; i++) {
478: jrow = ii[i];
479: VecPlaceArray(bmat->left,yy + bs*i);
480: n = ii[i+1] - jrow;
481: for (j=0; j<n; j++) {
482: VecPlaceArray(bmat->right,xx + bs*aj[jrow]);
483: MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
484: VecResetArray(bmat->right);
485: jrow++;
486: }
487: VecResetArray(bmat->left);
488: }
489: VecRestoreArray(x,&xx);
490: VecRestoreArray(y,&yy);
491: return(0);
492: }
494: PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y)
495: {
496: Mat_BlockMat *bmat = (Mat_BlockMat*)A->data;
498: PetscScalar *xx,*yy;
499: PetscInt *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j;
500: Mat *aa;
503: /*
504: Standard CSR multiply except each entry is a Mat
505: */
506: VecGetArray(x,&xx);
508: VecSet(y,0.0);
509: VecGetArray(y,&yy);
510: aj = bmat->j;
511: aa = bmat->a;
512: ii = bmat->i;
513: for (i=0; i<m; i++) {
514: jrow = ii[i];
515: n = ii[i+1] - jrow;
516: VecPlaceArray(bmat->left,yy + bs*i);
517: VecPlaceArray(bmat->middle,xx + bs*i);
518: /* if we ALWAYS required a diagonal entry then could remove this if test */
519: if (aj[jrow] == i) {
520: VecPlaceArray(bmat->right,xx + bs*aj[jrow]);
521: MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
522: VecResetArray(bmat->right);
523: jrow++;
524: n--;
525: }
526: for (j=0; j<n; j++) {
527: VecPlaceArray(bmat->right,xx + bs*aj[jrow]); /* upper triangular part */
528: MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
529: VecResetArray(bmat->right);
531: VecPlaceArray(bmat->right,yy + bs*aj[jrow]); /* lower triangular part */
532: MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);
533: VecResetArray(bmat->right);
534: jrow++;
535: }
536: VecResetArray(bmat->left);
537: VecResetArray(bmat->middle);
538: }
539: VecRestoreArray(x,&xx);
540: VecRestoreArray(y,&yy);
541: return(0);
542: }
544: static PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z)
545: {
547: return(0);
548: }
550: static PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y)
551: {
553: return(0);
554: }
556: static PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z)
557: {
559: return(0);
560: }
562: /*
563: Adds diagonal pointers to sparse matrix structure.
564: */
565: static PetscErrorCode MatMarkDiagonal_BlockMat(Mat A)
566: {
567: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
569: PetscInt i,j,mbs = A->rmap->n/A->rmap->bs;
572: if (!a->diag) {
573: PetscMalloc1(mbs,&a->diag);
574: }
575: for (i=0; i<mbs; i++) {
576: a->diag[i] = a->i[i+1];
577: for (j=a->i[i]; j<a->i[i+1]; j++) {
578: if (a->j[j] == i) {
579: a->diag[i] = j;
580: break;
581: }
582: }
583: }
584: return(0);
585: }
587: static PetscErrorCode MatCreateSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
588: {
589: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
590: Mat_SeqAIJ *c;
592: PetscInt i,k,first,step,lensi,nrows,ncols;
593: PetscInt *j_new,*i_new,*aj = a->j,*ailen = a->ilen;
594: PetscScalar *a_new;
595: Mat C,*aa = a->a;
596: PetscBool stride,equal;
599: ISEqual(isrow,iscol,&equal);
600: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices");
601: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
602: if (!stride) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for stride indices");
603: ISStrideGetInfo(iscol,&first,&step);
604: if (step != A->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only select one entry from each block");
606: ISGetLocalSize(isrow,&nrows);
607: ncols = nrows;
609: /* create submatrix */
610: if (scall == MAT_REUSE_MATRIX) {
611: PetscInt n_cols,n_rows;
612: C = *B;
613: MatGetSize(C,&n_rows,&n_cols);
614: if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
615: MatZeroEntries(C);
616: } else {
617: MatCreate(PetscObjectComm((PetscObject)A),&C);
618: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
619: if (A->symmetric) {
620: MatSetType(C,MATSEQSBAIJ);
621: } else {
622: MatSetType(C,MATSEQAIJ);
623: }
624: MatSeqAIJSetPreallocation(C,0,ailen);
625: MatSeqSBAIJSetPreallocation(C,1,0,ailen);
626: }
627: c = (Mat_SeqAIJ*)C->data;
629: /* loop over rows inserting into submatrix */
630: a_new = c->a;
631: j_new = c->j;
632: i_new = c->i;
634: for (i=0; i<nrows; i++) {
635: lensi = ailen[i];
636: for (k=0; k<lensi; k++) {
637: *j_new++ = *aj++;
638: MatGetValue(*aa++,first,first,a_new++);
639: }
640: i_new[i+1] = i_new[i] + lensi;
641: c->ilen[i] = lensi;
642: }
644: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
645: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
646: *B = C;
647: return(0);
648: }
650: static PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode)
651: {
652: Mat_BlockMat *a = (Mat_BlockMat*)A->data;
654: PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
655: PetscInt m = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0;
656: Mat *aa = a->a,*ap;
659: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
661: if (m) rmax = ailen[0]; /* determine row with most nonzeros */
662: for (i=1; i<m; i++) {
663: /* move each row back by the amount of empty slots (fshift) before it*/
664: fshift += imax[i-1] - ailen[i-1];
665: rmax = PetscMax(rmax,ailen[i]);
666: if (fshift) {
667: ip = aj + ai[i];
668: ap = aa + ai[i];
669: N = ailen[i];
670: for (j=0; j<N; j++) {
671: ip[j-fshift] = ip[j];
672: ap[j-fshift] = ap[j];
673: }
674: }
675: ai[i] = ai[i-1] + ailen[i-1];
676: }
677: if (m) {
678: fshift += imax[m-1] - ailen[m-1];
679: ai[m] = ai[m-1] + ailen[m-1];
680: }
681: /* reset ilen and imax for each row */
682: for (i=0; i<m; i++) {
683: ailen[i] = imax[i] = ai[i+1] - ai[i];
684: }
685: a->nz = ai[m];
686: for (i=0; i<a->nz; i++) {
687: #if defined(PETSC_USE_DEBUG)
688: if (!aa[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz);
689: #endif
690: MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);
691: MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);
692: }
693: PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n/A->cmap->bs,fshift,a->nz);
694: PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);
695: PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);
697: A->info.mallocs += a->reallocs;
698: a->reallocs = 0;
699: A->info.nz_unneeded = (double)fshift;
700: a->rmax = rmax;
701: MatMarkDiagonal_BlockMat(A);
702: return(0);
703: }
705: static PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt,PetscBool flg)
706: {
709: if (opt == MAT_SYMMETRIC && flg) {
710: A->ops->sor = MatSOR_BlockMat_Symmetric;
711: A->ops->mult = MatMult_BlockMat_Symmetric;
712: } else {
713: PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]);
714: }
715: return(0);
716: }
719: static struct _MatOps MatOps_Values = {MatSetValues_BlockMat,
720: 0,
721: 0,
722: MatMult_BlockMat,
723: /* 4*/ MatMultAdd_BlockMat,
724: MatMultTranspose_BlockMat,
725: MatMultTransposeAdd_BlockMat,
726: 0,
727: 0,
728: 0,
729: /* 10*/ 0,
730: 0,
731: 0,
732: MatSOR_BlockMat,
733: 0,
734: /* 15*/ 0,
735: 0,
736: 0,
737: 0,
738: 0,
739: /* 20*/ 0,
740: MatAssemblyEnd_BlockMat,
741: MatSetOption_BlockMat,
742: 0,
743: /* 24*/ 0,
744: 0,
745: 0,
746: 0,
747: 0,
748: /* 29*/ 0,
749: 0,
750: 0,
751: 0,
752: 0,
753: /* 34*/ 0,
754: 0,
755: 0,
756: 0,
757: 0,
758: /* 39*/ 0,
759: 0,
760: 0,
761: 0,
762: 0,
763: /* 44*/ 0,
764: 0,
765: MatShift_Basic,
766: 0,
767: 0,
768: /* 49*/ 0,
769: 0,
770: 0,
771: 0,
772: 0,
773: /* 54*/ 0,
774: 0,
775: 0,
776: 0,
777: 0,
778: /* 59*/ MatCreateSubMatrix_BlockMat,
779: MatDestroy_BlockMat,
780: MatView_BlockMat,
781: 0,
782: 0,
783: /* 64*/ 0,
784: 0,
785: 0,
786: 0,
787: 0,
788: /* 69*/ 0,
789: 0,
790: 0,
791: 0,
792: 0,
793: /* 74*/ 0,
794: 0,
795: 0,
796: 0,
797: 0,
798: /* 79*/ 0,
799: 0,
800: 0,
801: 0,
802: MatLoad_BlockMat,
803: /* 84*/ 0,
804: 0,
805: 0,
806: 0,
807: 0,
808: /* 89*/ 0,
809: 0,
810: 0,
811: 0,
812: 0,
813: /* 94*/ 0,
814: 0,
815: 0,
816: 0,
817: 0,
818: /* 99*/ 0,
819: 0,
820: 0,
821: 0,
822: 0,
823: /*104*/ 0,
824: 0,
825: 0,
826: 0,
827: 0,
828: /*109*/ 0,
829: 0,
830: 0,
831: 0,
832: 0,
833: /*114*/ 0,
834: 0,
835: 0,
836: 0,
837: 0,
838: /*119*/ 0,
839: 0,
840: 0,
841: 0,
842: 0,
843: /*124*/ 0,
844: 0,
845: 0,
846: 0,
847: 0,
848: /*129*/ 0,
849: 0,
850: 0,
851: 0,
852: 0,
853: /*134*/ 0,
854: 0,
855: 0,
856: 0,
857: 0,
858: /*139*/ 0,
859: 0,
860: 0
861: };
863: /*@C
864: MatBlockMatSetPreallocation - For good matrix assembly performance
865: the user should preallocate the matrix storage by setting the parameter nz
866: (or the array nnz). By setting these parameters accurately, performance
867: during matrix assembly can be increased by more than a factor of 50.
869: Collective
871: Input Parameters:
872: + B - The matrix
873: . bs - size of each block in matrix
874: . nz - number of nonzeros per block row (same for all rows)
875: - nnz - array containing the number of nonzeros in the various block rows
876: (possibly different for each row) or NULL
878: Notes:
879: If nnz is given then nz is ignored
881: Specify the preallocated storage with either nz or nnz (not both).
882: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
883: allocation. For large problems you MUST preallocate memory or you
884: will get TERRIBLE performance, see the users' manual chapter on matrices.
886: Level: intermediate
888: .seealso: MatCreate(), MatCreateBlockMat(), MatSetValues()
890: @*/
891: PetscErrorCode MatBlockMatSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
892: {
896: PetscTryMethod(B,"MatBlockMatSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
897: return(0);
898: }
900: static PetscErrorCode MatBlockMatSetPreallocation_BlockMat(Mat A,PetscInt bs,PetscInt nz,PetscInt *nnz)
901: {
902: Mat_BlockMat *bmat = (Mat_BlockMat*)A->data;
904: PetscInt i;
907: PetscLayoutSetBlockSize(A->rmap,bs);
908: PetscLayoutSetBlockSize(A->cmap,bs);
909: PetscLayoutSetUp(A->rmap);
910: PetscLayoutSetUp(A->cmap);
911: PetscLayoutGetBlockSize(A->rmap,&bs);
913: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
914: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz);
915: if (nnz) {
916: for (i=0; i<A->rmap->n/bs; i++) {
917: if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]);
918: if (nnz[i] > A->cmap->n/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],A->cmap->n/bs);
919: }
920: }
921: bmat->mbs = A->rmap->n/bs;
923: VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->right);
924: VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->middle);
925: VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);
927: if (!bmat->imax) {
928: PetscMalloc2(A->rmap->n,&bmat->imax,A->rmap->n,&bmat->ilen);
929: PetscLogObjectMemory((PetscObject)A,2*A->rmap->n*sizeof(PetscInt));
930: }
931: if (nnz) {
932: nz = 0;
933: for (i=0; i<A->rmap->n/A->rmap->bs; i++) {
934: bmat->imax[i] = nnz[i];
935: nz += nnz[i];
936: }
937: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently requires block row by row preallocation");
939: /* bmat->ilen will count nonzeros in each row so far. */
940: for (i=0; i<bmat->mbs; i++) bmat->ilen[i] = 0;
942: /* allocate the matrix space */
943: MatSeqXAIJFreeAIJ(A,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);
944: PetscMalloc3(nz,&bmat->a,nz,&bmat->j,A->rmap->n+1,&bmat->i);
945: PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
946: bmat->i[0] = 0;
947: for (i=1; i<bmat->mbs+1; i++) {
948: bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1];
949: }
950: bmat->singlemalloc = PETSC_TRUE;
951: bmat->free_a = PETSC_TRUE;
952: bmat->free_ij = PETSC_TRUE;
954: bmat->nz = 0;
955: bmat->maxnz = nz;
956: A->info.nz_unneeded = (double)bmat->maxnz;
957: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
958: return(0);
959: }
961: /*MC
962: MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix
963: consisting of (usually) sparse blocks.
965: Level: advanced
967: .seealso: MatCreateBlockMat()
969: M*/
971: PETSC_EXTERN PetscErrorCode MatCreate_BlockMat(Mat A)
972: {
973: Mat_BlockMat *b;
977: PetscNewLog(A,&b);
978: A->data = (void*)b;
979: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
981: A->assembled = PETSC_TRUE;
982: A->preallocated = PETSC_FALSE;
983: PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);
985: PetscObjectComposeFunction((PetscObject)A,"MatBlockMatSetPreallocation_C",MatBlockMatSetPreallocation_BlockMat);
986: return(0);
987: }
989: /*@C
990: MatCreateBlockMat - Creates a new matrix in which each block contains a uniform-size sequential Mat object
992: Collective
994: Input Parameters:
995: + comm - MPI communicator
996: . m - number of rows
997: . n - number of columns
998: . bs - size of each submatrix
999: . nz - expected maximum number of nonzero blocks in row (use PETSC_DEFAULT if not known)
1000: - nnz - expected number of nonzers per block row if known (use NULL otherwise)
1003: Output Parameter:
1004: . A - the matrix
1006: Level: intermediate
1008: Notes:
1009: Matrices of this type are nominally-sparse matrices in which each "entry" is a Mat object. Each Mat must
1010: have the same size and be sequential. The local and global sizes must be compatible with this decomposition.
1012: For matrices containing parallel submatrices and variable block sizes, see MATNEST.
1014: .seealso: MATBLOCKMAT, MatCreateNest()
1015: @*/
1016: PetscErrorCode MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,PetscInt nz,PetscInt *nnz, Mat *A)
1017: {
1021: MatCreate(comm,A);
1022: MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
1023: MatSetType(*A,MATBLOCKMAT);
1024: MatBlockMatSetPreallocation(*A,bs,nz,nnz);
1025: return(0);
1026: }