Actual source code: sbaij.c
petsc-3.12.5 2020-03-29
2: /*
3: Defines the basic matrix operations for the SBAIJ (compressed row)
4: matrix storage format.
5: */
6: #include <../src/mat/impls/baij/seq/baij.h>
7: #include <../src/mat/impls/sbaij/seq/sbaij.h>
8: #include <petscblaslapack.h>
10: #include <../src/mat/impls/sbaij/seq/relax.h>
11: #define USESHORT
12: #include <../src/mat/impls/sbaij/seq/relax.h>
14: #if defined(PETSC_HAVE_ELEMENTAL)
15: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
16: #endif
17: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat,MatType,MatReuse,Mat*);
19: /*
20: Checks for missing diagonals
21: */
22: PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd)
23: {
24: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
26: PetscInt *diag,*ii = a->i,i;
29: MatMarkDiagonal_SeqSBAIJ(A);
30: *missing = PETSC_FALSE;
31: if (A->rmap->n > 0 && !ii) {
32: *missing = PETSC_TRUE;
33: if (dd) *dd = 0;
34: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
35: } else {
36: diag = a->diag;
37: for (i=0; i<a->mbs; i++) {
38: if (diag[i] >= ii[i+1]) {
39: *missing = PETSC_TRUE;
40: if (dd) *dd = i;
41: break;
42: }
43: }
44: }
45: return(0);
46: }
48: PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
49: {
50: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
52: PetscInt i,j;
55: if (!a->diag) {
56: PetscMalloc1(a->mbs,&a->diag);
57: PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt));
58: a->free_diag = PETSC_TRUE;
59: }
60: for (i=0; i<a->mbs; i++) {
61: a->diag[i] = a->i[i+1];
62: for (j=a->i[i]; j<a->i[i+1]; j++) {
63: if (a->j[j] == i) {
64: a->diag[i] = j;
65: break;
66: }
67: }
68: }
69: return(0);
70: }
72: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done)
73: {
74: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
76: PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt;
77: PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja;
80: *nn = n;
81: if (!ia) return(0);
82: if (symmetric) {
83: MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja);
84: nz = tia[n];
85: } else {
86: tia = a->i; tja = a->j;
87: }
89: if (!blockcompressed && bs > 1) {
90: (*nn) *= bs;
91: /* malloc & create the natural set of indices */
92: PetscMalloc1((n+1)*bs,ia);
93: if (n) {
94: (*ia)[0] = oshift;
95: for (j=1; j<bs; j++) {
96: (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1];
97: }
98: }
100: for (i=1; i<n; i++) {
101: (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1];
102: for (j=1; j<bs; j++) {
103: (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1];
104: }
105: }
106: if (n) {
107: (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1];
108: }
110: if (inja) {
111: PetscMalloc1(nz*bs*bs,ja);
112: cnt = 0;
113: for (i=0; i<n; i++) {
114: for (j=0; j<bs; j++) {
115: for (k=tia[i]; k<tia[i+1]; k++) {
116: for (l=0; l<bs; l++) {
117: (*ja)[cnt++] = bs*tja[k] + l;
118: }
119: }
120: }
121: }
122: }
124: if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
125: PetscFree(tia);
126: PetscFree(tja);
127: }
128: } else if (oshift == 1) {
129: if (symmetric) {
130: nz = tia[A->rmap->n/bs];
131: /* add 1 to i and j indices */
132: for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1;
133: *ia = tia;
134: if (ja) {
135: for (i=0; i<nz; i++) tja[i] = tja[i] + 1;
136: *ja = tja;
137: }
138: } else {
139: nz = a->i[A->rmap->n/bs];
140: /* malloc space and add 1 to i and j indices */
141: PetscMalloc1(A->rmap->n/bs+1,ia);
142: for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1;
143: if (ja) {
144: PetscMalloc1(nz,ja);
145: for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1;
146: }
147: }
148: } else {
149: *ia = tia;
150: if (ja) *ja = tja;
151: }
152: return(0);
153: }
155: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
156: {
160: if (!ia) return(0);
161: if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
162: PetscFree(*ia);
163: if (ja) {PetscFree(*ja);}
164: }
165: return(0);
166: }
168: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
169: {
170: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
174: #if defined(PETSC_USE_LOG)
175: PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz);
176: #endif
177: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
178: if (a->free_diag) {PetscFree(a->diag);}
179: ISDestroy(&a->row);
180: ISDestroy(&a->col);
181: ISDestroy(&a->icol);
182: PetscFree(a->idiag);
183: PetscFree(a->inode.size);
184: if (a->free_imax_ilen) {PetscFree2(a->imax,a->ilen);}
185: PetscFree(a->solve_work);
186: PetscFree(a->sor_work);
187: PetscFree(a->solves_work);
188: PetscFree(a->mult_work);
189: PetscFree(a->saved_values);
190: if (a->free_jshort) {PetscFree(a->jshort);}
191: PetscFree(a->inew);
192: MatDestroy(&a->parent);
193: PetscFree(A->data);
195: PetscObjectChangeTypeName((PetscObject)A,0);
196: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
197: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
198: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL);
199: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL);
200: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL);
201: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL);
202: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL);
203: #if defined(PETSC_HAVE_ELEMENTAL)
204: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL);
205: #endif
206: return(0);
207: }
209: PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg)
210: {
211: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
215: switch (op) {
216: case MAT_ROW_ORIENTED:
217: a->roworiented = flg;
218: break;
219: case MAT_KEEP_NONZERO_PATTERN:
220: a->keepnonzeropattern = flg;
221: break;
222: case MAT_NEW_NONZERO_LOCATIONS:
223: a->nonew = (flg ? 0 : 1);
224: break;
225: case MAT_NEW_NONZERO_LOCATION_ERR:
226: a->nonew = (flg ? -1 : 0);
227: break;
228: case MAT_NEW_NONZERO_ALLOCATION_ERR:
229: a->nonew = (flg ? -2 : 0);
230: break;
231: case MAT_UNUSED_NONZERO_LOCATION_ERR:
232: a->nounused = (flg ? -1 : 0);
233: break;
234: case MAT_NEW_DIAGONALS:
235: case MAT_IGNORE_OFF_PROC_ENTRIES:
236: case MAT_USE_HASH_TABLE:
237: case MAT_SORTED_FULL:
238: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
239: break;
240: case MAT_HERMITIAN:
241: #if defined(PETSC_USE_COMPLEX) /* MAT_HERMITIAN is a synonym for MAT_SYMMETRIC with reals */
242: if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatAssemblyEnd() first");
243: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
244: A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort;
245: } else if (A->cmap->bs == 1) {
246: A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian;
247: } else if (!A->symmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1");
248: #endif
249: break;
250: case MAT_SPD:
251: /* These options are handled directly by MatSetOption() */
252: break;
253: case MAT_SYMMETRIC:
254: case MAT_STRUCTURALLY_SYMMETRIC:
255: case MAT_SYMMETRY_ETERNAL:
256: case MAT_STRUCTURE_ONLY:
257: /* These options are handled directly by MatSetOption() */
258: break;
259: case MAT_IGNORE_LOWER_TRIANGULAR:
260: a->ignore_ltriangular = flg;
261: break;
262: case MAT_ERROR_LOWER_TRIANGULAR:
263: a->ignore_ltriangular = flg;
264: break;
265: case MAT_GETROW_UPPERTRIANGULAR:
266: a->getrow_utriangular = flg;
267: break;
268: case MAT_SUBMAT_SINGLEIS:
269: break;
270: default:
271: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
272: }
273: return(0);
274: }
276: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
277: {
278: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
282: if (A && !a->getrow_utriangular) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()");
284: /* Get the upper triangular part of the row */
285: MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);
286: return(0);
287: }
289: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
290: {
294: if (idx) {PetscFree(*idx);}
295: if (v) {PetscFree(*v);}
296: return(0);
297: }
299: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
300: {
301: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
304: a->getrow_utriangular = PETSC_TRUE;
305: return(0);
306: }
308: PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
309: {
310: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
313: a->getrow_utriangular = PETSC_FALSE;
314: return(0);
315: }
317: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B)
318: {
322: if (reuse == MAT_INITIAL_MATRIX) {
323: MatDuplicate(A,MAT_COPY_VALUES,B);
324: } else if (reuse == MAT_REUSE_MATRIX) {
325: MatCopy(A,*B,SAME_NONZERO_PATTERN);
326: }
327: return(0);
328: }
330: PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
331: {
332: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
333: PetscErrorCode ierr;
334: PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
335: PetscViewerFormat format;
336: PetscInt *diag;
339: PetscViewerGetFormat(viewer,&format);
340: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
341: PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);
342: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
343: Mat aij;
344: const char *matname;
346: if (A->factortype && bs>1) {
347: PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");
348: return(0);
349: }
350: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);
351: PetscObjectGetName((PetscObject)A,&matname);
352: PetscObjectSetName((PetscObject)aij,matname);
353: MatView(aij,viewer);
354: MatDestroy(&aij);
355: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
356: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
357: for (i=0; i<a->mbs; i++) {
358: for (j=0; j<bs; j++) {
359: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
360: for (k=a->i[i]; k<a->i[i+1]; k++) {
361: for (l=0; l<bs; l++) {
362: #if defined(PETSC_USE_COMPLEX)
363: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
364: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
365: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
366: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
367: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
368: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
369: } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
370: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
371: }
372: #else
373: if (a->a[bs2*k + l*bs + j] != 0.0) {
374: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
375: }
376: #endif
377: }
378: }
379: PetscViewerASCIIPrintf(viewer,"\n");
380: }
381: }
382: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
383: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
384: return(0);
385: } else {
386: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
387: if (A->factortype) { /* for factored matrix */
388: if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet");
390: diag=a->diag;
391: for (i=0; i<a->mbs; i++) { /* for row block i */
392: PetscViewerASCIIPrintf(viewer,"row %D:",i);
393: /* diagonal entry */
394: #if defined(PETSC_USE_COMPLEX)
395: if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
396: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),(double)PetscImaginaryPart(1.0/a->a[diag[i]]));
397: } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
398: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),-(double)PetscImaginaryPart(1.0/a->a[diag[i]]));
399: } else {
400: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));
401: }
402: #else
403: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));
404: #endif
405: /* off-diagonal entries */
406: for (k=a->i[i]; k<a->i[i+1]-1; k++) {
407: #if defined(PETSC_USE_COMPLEX)
408: if (PetscImaginaryPart(a->a[k]) > 0.0) {
409: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));
410: } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
411: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));
412: } else {
413: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));
414: }
415: #else
416: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);
417: #endif
418: }
419: PetscViewerASCIIPrintf(viewer,"\n");
420: }
422: } else { /* for non-factored matrix */
423: for (i=0; i<a->mbs; i++) { /* for row block i */
424: for (j=0; j<bs; j++) { /* for row bs*i + j */
425: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
426: for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */
427: for (l=0; l<bs; l++) { /* for column */
428: #if defined(PETSC_USE_COMPLEX)
429: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
430: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
431: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
432: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
433: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
434: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
435: } else {
436: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
437: }
438: #else
439: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
440: #endif
441: }
442: }
443: PetscViewerASCIIPrintf(viewer,"\n");
444: }
445: }
446: }
447: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
448: }
449: PetscViewerFlush(viewer);
450: return(0);
451: }
453: #include <petscdraw.h>
454: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
455: {
456: Mat A = (Mat) Aa;
457: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
459: PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
460: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
461: MatScalar *aa;
462: PetscViewer viewer;
465: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
466: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
468: /* loop over matrix elements drawing boxes */
470: PetscDrawCollectiveBegin(draw);
471: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
472: /* Blue for negative, Cyan for zero and Red for positive */
473: color = PETSC_DRAW_BLUE;
474: for (i=0,row=0; i<mbs; i++,row+=bs) {
475: for (j=a->i[i]; j<a->i[i+1]; j++) {
476: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
477: x_l = a->j[j]*bs; x_r = x_l + 1.0;
478: aa = a->a + j*bs2;
479: for (k=0; k<bs; k++) {
480: for (l=0; l<bs; l++) {
481: if (PetscRealPart(*aa++) >= 0.) continue;
482: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
483: }
484: }
485: }
486: }
487: color = PETSC_DRAW_CYAN;
488: for (i=0,row=0; i<mbs; i++,row+=bs) {
489: for (j=a->i[i]; j<a->i[i+1]; j++) {
490: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
491: x_l = a->j[j]*bs; x_r = x_l + 1.0;
492: aa = a->a + j*bs2;
493: for (k=0; k<bs; k++) {
494: for (l=0; l<bs; l++) {
495: if (PetscRealPart(*aa++) != 0.) continue;
496: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
497: }
498: }
499: }
500: }
501: color = PETSC_DRAW_RED;
502: for (i=0,row=0; i<mbs; i++,row+=bs) {
503: for (j=a->i[i]; j<a->i[i+1]; j++) {
504: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
505: x_l = a->j[j]*bs; x_r = x_l + 1.0;
506: aa = a->a + j*bs2;
507: for (k=0; k<bs; k++) {
508: for (l=0; l<bs; l++) {
509: if (PetscRealPart(*aa++) <= 0.) continue;
510: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
511: }
512: }
513: }
514: }
515: PetscDrawCollectiveEnd(draw);
516: return(0);
517: }
519: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
520: {
522: PetscReal xl,yl,xr,yr,w,h;
523: PetscDraw draw;
524: PetscBool isnull;
527: PetscViewerDrawGetDraw(viewer,0,&draw);
528: PetscDrawIsNull(draw,&isnull);
529: if (isnull) return(0);
531: xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
532: xr += w; yr += h; xl = -w; yl = -h;
533: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
534: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
535: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
536: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
537: PetscDrawSave(draw);
538: return(0);
539: }
541: PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
542: {
544: PetscBool iascii,isdraw;
545: FILE *file = 0;
548: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
549: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
550: if (iascii) {
551: MatView_SeqSBAIJ_ASCII(A,viewer);
552: } else if (isdraw) {
553: MatView_SeqSBAIJ_Draw(A,viewer);
554: } else {
555: Mat B;
556: const char *matname;
557: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
558: PetscObjectGetName((PetscObject)A,&matname);
559: PetscObjectSetName((PetscObject)B,matname);
560: MatView(B,viewer);
561: MatDestroy(&B);
562: PetscViewerBinaryGetInfoPointer(viewer,&file);
563: if (file) {
564: fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs);
565: }
566: }
567: return(0);
568: }
571: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
572: {
573: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
574: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
575: PetscInt *ai = a->i,*ailen = a->ilen;
576: PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
577: MatScalar *ap,*aa = a->a;
580: for (k=0; k<m; k++) { /* loop over rows */
581: row = im[k]; brow = row/bs;
582: if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
583: 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);
584: rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
585: nrow = ailen[brow];
586: for (l=0; l<n; l++) { /* loop over columns */
587: if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
588: 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);
589: col = in[l];
590: bcol = col/bs;
591: cidx = col%bs;
592: ridx = row%bs;
593: high = nrow;
594: low = 0; /* assume unsorted */
595: while (high-low > 5) {
596: t = (low+high)/2;
597: if (rp[t] > bcol) high = t;
598: else low = t;
599: }
600: for (i=low; i<high; i++) {
601: if (rp[i] > bcol) break;
602: if (rp[i] == bcol) {
603: *v++ = ap[bs2*i+bs*cidx+ridx];
604: goto finished;
605: }
606: }
607: *v++ = 0.0;
608: finished:;
609: }
610: }
611: return(0);
612: }
615: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
616: {
617: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
618: PetscErrorCode ierr;
619: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
620: PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
621: PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
622: PetscBool roworiented=a->roworiented;
623: const PetscScalar *value = v;
624: MatScalar *ap,*aa = a->a,*bap;
627: if (roworiented) stepval = (n-1)*bs;
628: else stepval = (m-1)*bs;
630: for (k=0; k<m; k++) { /* loop over added rows */
631: row = im[k];
632: if (row < 0) continue;
633: #if defined(PETSC_USE_DEBUG)
634: if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index row too large %D max %D",row,a->mbs-1);
635: #endif
636: rp = aj + ai[row];
637: ap = aa + bs2*ai[row];
638: rmax = imax[row];
639: nrow = ailen[row];
640: low = 0;
641: high = nrow;
642: for (l=0; l<n; l++) { /* loop over added columns */
643: if (in[l] < 0) continue;
644: col = in[l];
645: #if defined(PETSC_USE_DEBUG)
646: if (col >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index column too large %D max %D",col,a->nbs-1);
647: #endif
648: if (col < row) {
649: if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
650: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
651: }
652: if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
653: else value = v + l*(stepval+bs)*bs + k*bs;
655: if (col <= lastcol) low = 0;
656: else high = nrow;
658: lastcol = col;
659: while (high-low > 7) {
660: t = (low+high)/2;
661: if (rp[t] > col) high = t;
662: else low = t;
663: }
664: for (i=low; i<high; i++) {
665: if (rp[i] > col) break;
666: if (rp[i] == col) {
667: bap = ap + bs2*i;
668: if (roworiented) {
669: if (is == ADD_VALUES) {
670: for (ii=0; ii<bs; ii++,value+=stepval) {
671: for (jj=ii; jj<bs2; jj+=bs) {
672: bap[jj] += *value++;
673: }
674: }
675: } else {
676: for (ii=0; ii<bs; ii++,value+=stepval) {
677: for (jj=ii; jj<bs2; jj+=bs) {
678: bap[jj] = *value++;
679: }
680: }
681: }
682: } else {
683: if (is == ADD_VALUES) {
684: for (ii=0; ii<bs; ii++,value+=stepval) {
685: for (jj=0; jj<bs; jj++) {
686: *bap++ += *value++;
687: }
688: }
689: } else {
690: for (ii=0; ii<bs; ii++,value+=stepval) {
691: for (jj=0; jj<bs; jj++) {
692: *bap++ = *value++;
693: }
694: }
695: }
696: }
697: goto noinsert2;
698: }
699: }
700: if (nonew == 1) goto noinsert2;
701: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new block index nonzero block (%D, %D) in the matrix", row, col);
702: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
703: N = nrow++ - 1; high++;
704: /* shift up all the later entries in this row */
705: PetscArraymove(rp+i+1,rp+i,N-i+1);
706: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
707: PetscArrayzero(ap+bs2*i,bs2);
708: rp[i] = col;
709: bap = ap + bs2*i;
710: if (roworiented) {
711: for (ii=0; ii<bs; ii++,value+=stepval) {
712: for (jj=ii; jj<bs2; jj+=bs) {
713: bap[jj] = *value++;
714: }
715: }
716: } else {
717: for (ii=0; ii<bs; ii++,value+=stepval) {
718: for (jj=0; jj<bs; jj++) {
719: *bap++ = *value++;
720: }
721: }
722: }
723: noinsert2:;
724: low = i;
725: }
726: ailen[row] = nrow;
727: }
728: return(0);
729: }
731: /*
732: This is not yet used
733: */
734: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
735: {
736: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
738: const PetscInt *ai = a->i, *aj = a->j,*cols;
739: PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
740: PetscBool flag;
743: PetscMalloc1(m,&ns);
744: while (i < m) {
745: nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
746: /* Limits the number of elements in a node to 'a->inode.limit' */
747: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
748: nzy = ai[j+1] - ai[j];
749: if (nzy != (nzx - j + i)) break;
750: PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);
751: if (!flag) break;
752: }
753: ns[node_count++] = blk_size;
755: i = j;
756: }
757: if (!a->inode.size && m && node_count > .9*m) {
758: PetscFree(ns);
759: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
760: } else {
761: a->inode.node_count = node_count;
763: PetscMalloc1(node_count,&a->inode.size);
764: PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));
765: PetscArraycpy(a->inode.size,ns,node_count);
766: PetscFree(ns);
767: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
769: /* count collections of adjacent columns in each inode */
770: row = 0;
771: cnt = 0;
772: for (i=0; i<node_count; i++) {
773: cols = aj + ai[row] + a->inode.size[i];
774: nz = ai[row+1] - ai[row] - a->inode.size[i];
775: for (j=1; j<nz; j++) {
776: if (cols[j] != cols[j-1]+1) cnt++;
777: }
778: cnt++;
779: row += a->inode.size[i];
780: }
781: PetscMalloc1(2*cnt,&counts);
782: cnt = 0;
783: row = 0;
784: for (i=0; i<node_count; i++) {
785: cols = aj + ai[row] + a->inode.size[i];
786: counts[2*cnt] = cols[0];
787: nz = ai[row+1] - ai[row] - a->inode.size[i];
788: cnt2 = 1;
789: for (j=1; j<nz; j++) {
790: if (cols[j] != cols[j-1]+1) {
791: counts[2*(cnt++)+1] = cnt2;
792: counts[2*cnt] = cols[j];
793: cnt2 = 1;
794: } else cnt2++;
795: }
796: counts[2*(cnt++)+1] = cnt2;
797: row += a->inode.size[i];
798: }
799: PetscIntView(2*cnt,counts,0);
800: }
801: return(0);
802: }
804: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
805: {
806: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
808: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
809: PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
810: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
811: MatScalar *aa = a->a,*ap;
814: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
816: if (m) rmax = ailen[0];
817: for (i=1; i<mbs; i++) {
818: /* move each row back by the amount of empty slots (fshift) before it*/
819: fshift += imax[i-1] - ailen[i-1];
820: rmax = PetscMax(rmax,ailen[i]);
821: if (fshift) {
822: ip = aj + ai[i];
823: ap = aa + bs2*ai[i];
824: N = ailen[i];
825: PetscArraymove(ip-fshift,ip,N);
826: PetscArraymove(ap-bs2*fshift,ap,bs2*N);
827: }
828: ai[i] = ai[i-1] + ailen[i-1];
829: }
830: if (mbs) {
831: fshift += imax[mbs-1] - ailen[mbs-1];
832: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
833: }
834: /* reset ilen and imax for each row */
835: for (i=0; i<mbs; i++) {
836: ailen[i] = imax[i] = ai[i+1] - ai[i];
837: }
838: a->nz = ai[mbs];
840: /* diagonals may have moved, reset it */
841: if (a->diag) {
842: PetscArraycpy(a->diag,ai,mbs);
843: }
844: if (fshift && a->nounused == -1) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2);
846: PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);
847: PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
848: PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);
850: A->info.mallocs += a->reallocs;
851: a->reallocs = 0;
852: A->info.nz_unneeded = (PetscReal)fshift*bs2;
853: a->idiagvalid = PETSC_FALSE;
854: a->rmax = rmax;
856: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
857: if (a->jshort && a->free_jshort) {
858: /* when matrix data structure is changed, previous jshort must be replaced */
859: PetscFree(a->jshort);
860: }
861: PetscMalloc1(a->i[A->rmap->n],&a->jshort);
862: PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));
863: for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
864: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
865: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
866: a->free_jshort = PETSC_TRUE;
867: }
868: return(0);
869: }
871: /*
872: This function returns an array of flags which indicate the locations of contiguous
873: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
874: then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
875: Assume: sizes should be long enough to hold all the values.
876: */
877: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
878: {
879: PetscInt i,j,k,row;
880: PetscBool flg;
883: for (i=0,j=0; i<n; j++) {
884: row = idx[i];
885: if (row%bs!=0) { /* Not the begining of a block */
886: sizes[j] = 1;
887: i++;
888: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
889: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
890: i++;
891: } else { /* Begining of the block, so check if the complete block exists */
892: flg = PETSC_TRUE;
893: for (k=1; k<bs; k++) {
894: if (row+k != idx[i+k]) { /* break in the block */
895: flg = PETSC_FALSE;
896: break;
897: }
898: }
899: if (flg) { /* No break in the bs */
900: sizes[j] = bs;
901: i += bs;
902: } else {
903: sizes[j] = 1;
904: i++;
905: }
906: }
907: }
908: *bs_max = j;
909: return(0);
910: }
913: /* Only add/insert a(i,j) with i<=j (blocks).
914: Any a(i,j) with i>j input by user is ingored.
915: */
917: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
918: {
919: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
921: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
922: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
923: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
924: PetscInt ridx,cidx,bs2=a->bs2;
925: MatScalar *ap,value,*aa=a->a,*bap;
928: for (k=0; k<m; k++) { /* loop over added rows */
929: row = im[k]; /* row number */
930: brow = row/bs; /* block row number */
931: if (row < 0) continue;
932: #if defined(PETSC_USE_DEBUG)
933: 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);
934: #endif
935: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
936: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
937: rmax = imax[brow]; /* maximum space allocated for this row */
938: nrow = ailen[brow]; /* actual length of this row */
939: low = 0;
940: high = nrow;
941: for (l=0; l<n; l++) { /* loop over added columns */
942: if (in[l] < 0) continue;
943: #if defined(PETSC_USE_DEBUG)
944: 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);
945: #endif
946: col = in[l];
947: bcol = col/bs; /* block col number */
949: if (brow > bcol) {
950: if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
951: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
952: }
954: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
955: if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
956: /* element value a(k,l) */
957: if (roworiented) value = v[l + k*n];
958: else value = v[k + l*m];
960: /* move pointer bap to a(k,l) quickly and add/insert value */
961: if (col <= lastcol) low = 0;
962: else high = nrow;
964: lastcol = col;
965: while (high-low > 7) {
966: t = (low+high)/2;
967: if (rp[t] > bcol) high = t;
968: else low = t;
969: }
970: for (i=low; i<high; i++) {
971: if (rp[i] > bcol) break;
972: if (rp[i] == bcol) {
973: bap = ap + bs2*i + bs*cidx + ridx;
974: if (is == ADD_VALUES) *bap += value;
975: else *bap = value;
976: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
977: if (brow == bcol && ridx < cidx) {
978: bap = ap + bs2*i + bs*ridx + cidx;
979: if (is == ADD_VALUES) *bap += value;
980: else *bap = value;
981: }
982: goto noinsert1;
983: }
984: }
986: if (nonew == 1) goto noinsert1;
987: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
988: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
990: N = nrow++ - 1; high++;
991: /* shift up all the later entries in this row */
992: PetscArraymove(rp+i+1,rp+i,N-i+1);
993: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
994: PetscArrayzero(ap+bs2*i,bs2);
995: rp[i] = bcol;
996: ap[bs2*i + bs*cidx + ridx] = value;
997: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
998: if (brow == bcol && ridx < cidx) {
999: ap[bs2*i + bs*ridx + cidx] = value;
1000: }
1001: A->nonzerostate++;
1002: noinsert1:;
1003: low = i;
1004: }
1005: } /* end of loop over added columns */
1006: ailen[brow] = nrow;
1007: } /* end of loop over added rows */
1008: return(0);
1009: }
1011: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
1012: {
1013: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1014: Mat outA;
1016: PetscBool row_identity;
1019: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
1020: ISIdentity(row,&row_identity);
1021: if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1022: if (inA->rmap->bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */
1024: outA = inA;
1025: inA->factortype = MAT_FACTOR_ICC;
1026: PetscFree(inA->solvertype);
1027: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
1029: MatMarkDiagonal_SeqSBAIJ(inA);
1030: MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);
1032: PetscObjectReference((PetscObject)row);
1033: ISDestroy(&a->row);
1034: a->row = row;
1035: PetscObjectReference((PetscObject)row);
1036: ISDestroy(&a->col);
1037: a->col = row;
1039: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1040: if (a->icol) {ISInvertPermutation(row,PETSC_DECIDE, &a->icol);}
1041: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
1043: if (!a->solve_work) {
1044: PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
1045: PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
1046: }
1048: MatCholeskyFactorNumeric(outA,inA,info);
1049: return(0);
1050: }
1052: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1053: {
1054: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1055: PetscInt i,nz,n;
1059: nz = baij->maxnz;
1060: n = mat->cmap->n;
1061: for (i=0; i<nz; i++) baij->j[i] = indices[i];
1063: baij->nz = nz;
1064: for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1066: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1067: return(0);
1068: }
1070: /*@
1071: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1072: in the matrix.
1074: Input Parameters:
1075: + mat - the SeqSBAIJ matrix
1076: - indices - the column indices
1078: Level: advanced
1080: Notes:
1081: This can be called if you have precomputed the nonzero structure of the
1082: matrix and want to provide it to the matrix object to improve the performance
1083: of the MatSetValues() operation.
1085: You MUST have set the correct numbers of nonzeros per row in the call to
1086: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1088: MUST be called before any calls to MatSetValues()
1090: .seealso: MatCreateSeqSBAIJ
1091: @*/
1092: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1093: {
1099: PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
1100: return(0);
1101: }
1103: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1104: {
1106: PetscBool isbaij;
1109: PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");
1110: if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name);
1111: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1112: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1113: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1114: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1116: if (a->i[a->mbs] != b->i[b->mbs]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different");
1117: if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different");
1118: if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size");
1119: PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);
1120: PetscObjectStateIncrease((PetscObject)B);
1121: } else {
1122: MatGetRowUpperTriangular(A);
1123: MatCopy_Basic(A,B,str);
1124: MatRestoreRowUpperTriangular(A);
1125: }
1126: return(0);
1127: }
1129: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1130: {
1134: MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,0);
1135: return(0);
1136: }
1138: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1139: {
1140: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1143: *array = a->a;
1144: return(0);
1145: }
1147: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1148: {
1150: *array = NULL;
1151: return(0);
1152: }
1154: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1155: {
1156: PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1157: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1158: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1162: /* Set the number of nonzeros in the new matrix */
1163: MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
1164: return(0);
1165: }
1167: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1168: {
1169: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1171: PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1172: PetscBLASInt one = 1;
1175: if (str == SAME_NONZERO_PATTERN) {
1176: PetscScalar alpha = a;
1177: PetscBLASInt bnz;
1178: PetscBLASIntCast(x->nz*bs2,&bnz);
1179: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1180: PetscObjectStateIncrease((PetscObject)Y);
1181: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1182: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);
1183: MatAXPY_Basic(Y,a,X,str);
1184: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);
1185: } else {
1186: Mat B;
1187: PetscInt *nnz;
1188: if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
1189: MatGetRowUpperTriangular(X);
1190: MatGetRowUpperTriangular(Y);
1191: PetscMalloc1(Y->rmap->N,&nnz);
1192: MatCreate(PetscObjectComm((PetscObject)Y),&B);
1193: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
1194: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
1195: MatSetBlockSizesFromMats(B,Y,Y);
1196: MatSetType(B,((PetscObject)Y)->type_name);
1197: MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);
1198: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1200: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1202: MatHeaderReplace(Y,&B);
1203: PetscFree(nnz);
1204: MatRestoreRowUpperTriangular(X);
1205: MatRestoreRowUpperTriangular(Y);
1206: }
1207: return(0);
1208: }
1210: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1211: {
1213: *flg = PETSC_TRUE;
1214: return(0);
1215: }
1217: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1218: {
1220: *flg = PETSC_TRUE;
1221: return(0);
1222: }
1224: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1225: {
1227: *flg = PETSC_FALSE;
1228: return(0);
1229: }
1231: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1232: {
1233: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1234: PetscInt i,nz = a->bs2*a->i[a->mbs];
1235: MatScalar *aa = a->a;
1238: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1239: return(0);
1240: }
1242: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1243: {
1244: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1245: PetscInt i,nz = a->bs2*a->i[a->mbs];
1246: MatScalar *aa = a->a;
1249: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1250: return(0);
1251: }
1253: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1254: {
1255: Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1256: PetscErrorCode ierr;
1257: PetscInt i,j,k,count;
1258: PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1259: PetscScalar zero = 0.0;
1260: MatScalar *aa;
1261: const PetscScalar *xx;
1262: PetscScalar *bb;
1263: PetscBool *zeroed,vecs = PETSC_FALSE;
1266: /* fix right hand side if needed */
1267: if (x && b) {
1268: VecGetArrayRead(x,&xx);
1269: VecGetArray(b,&bb);
1270: vecs = PETSC_TRUE;
1271: }
1273: /* zero the columns */
1274: PetscCalloc1(A->rmap->n,&zeroed);
1275: for (i=0; i<is_n; i++) {
1276: if (is_idx[i] < 0 || is_idx[i] >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",is_idx[i]);
1277: zeroed[is_idx[i]] = PETSC_TRUE;
1278: }
1279: if (vecs) {
1280: for (i=0; i<A->rmap->N; i++) {
1281: row = i/bs;
1282: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1283: for (k=0; k<bs; k++) {
1284: col = bs*baij->j[j] + k;
1285: if (col <= i) continue;
1286: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1287: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1288: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1289: }
1290: }
1291: }
1292: for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1293: }
1295: for (i=0; i<A->rmap->N; i++) {
1296: if (!zeroed[i]) {
1297: row = i/bs;
1298: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1299: for (k=0; k<bs; k++) {
1300: col = bs*baij->j[j] + k;
1301: if (zeroed[col]) {
1302: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1303: aa[0] = 0.0;
1304: }
1305: }
1306: }
1307: }
1308: }
1309: PetscFree(zeroed);
1310: if (vecs) {
1311: VecRestoreArrayRead(x,&xx);
1312: VecRestoreArray(b,&bb);
1313: }
1315: /* zero the rows */
1316: for (i=0; i<is_n; i++) {
1317: row = is_idx[i];
1318: count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1319: aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1320: for (k=0; k<count; k++) {
1321: aa[0] = zero;
1322: aa += bs;
1323: }
1324: if (diag != 0.0) {
1325: (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
1326: }
1327: }
1328: MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);
1329: return(0);
1330: }
1332: PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1333: {
1335: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1338: if (!Y->preallocated || !aij->nz) {
1339: MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
1340: }
1341: MatShift_Basic(Y,a);
1342: return(0);
1343: }
1345: /* -------------------------------------------------------------------*/
1346: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1347: MatGetRow_SeqSBAIJ,
1348: MatRestoreRow_SeqSBAIJ,
1349: MatMult_SeqSBAIJ_N,
1350: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1351: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1352: MatMultAdd_SeqSBAIJ_N,
1353: 0,
1354: 0,
1355: 0,
1356: /* 10*/ 0,
1357: 0,
1358: MatCholeskyFactor_SeqSBAIJ,
1359: MatSOR_SeqSBAIJ,
1360: MatTranspose_SeqSBAIJ,
1361: /* 15*/ MatGetInfo_SeqSBAIJ,
1362: MatEqual_SeqSBAIJ,
1363: MatGetDiagonal_SeqSBAIJ,
1364: MatDiagonalScale_SeqSBAIJ,
1365: MatNorm_SeqSBAIJ,
1366: /* 20*/ 0,
1367: MatAssemblyEnd_SeqSBAIJ,
1368: MatSetOption_SeqSBAIJ,
1369: MatZeroEntries_SeqSBAIJ,
1370: /* 24*/ 0,
1371: 0,
1372: 0,
1373: 0,
1374: 0,
1375: /* 29*/ MatSetUp_SeqSBAIJ,
1376: 0,
1377: 0,
1378: 0,
1379: 0,
1380: /* 34*/ MatDuplicate_SeqSBAIJ,
1381: 0,
1382: 0,
1383: 0,
1384: MatICCFactor_SeqSBAIJ,
1385: /* 39*/ MatAXPY_SeqSBAIJ,
1386: MatCreateSubMatrices_SeqSBAIJ,
1387: MatIncreaseOverlap_SeqSBAIJ,
1388: MatGetValues_SeqSBAIJ,
1389: MatCopy_SeqSBAIJ,
1390: /* 44*/ 0,
1391: MatScale_SeqSBAIJ,
1392: MatShift_SeqSBAIJ,
1393: 0,
1394: MatZeroRowsColumns_SeqSBAIJ,
1395: /* 49*/ 0,
1396: MatGetRowIJ_SeqSBAIJ,
1397: MatRestoreRowIJ_SeqSBAIJ,
1398: 0,
1399: 0,
1400: /* 54*/ 0,
1401: 0,
1402: 0,
1403: 0,
1404: MatSetValuesBlocked_SeqSBAIJ,
1405: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1406: 0,
1407: 0,
1408: 0,
1409: 0,
1410: /* 64*/ 0,
1411: 0,
1412: 0,
1413: 0,
1414: 0,
1415: /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1416: 0,
1417: MatConvert_MPISBAIJ_Basic,
1418: 0,
1419: 0,
1420: /* 74*/ 0,
1421: 0,
1422: 0,
1423: 0,
1424: 0,
1425: /* 79*/ 0,
1426: 0,
1427: 0,
1428: MatGetInertia_SeqSBAIJ,
1429: MatLoad_SeqSBAIJ,
1430: /* 84*/ MatIsSymmetric_SeqSBAIJ,
1431: MatIsHermitian_SeqSBAIJ,
1432: MatIsStructurallySymmetric_SeqSBAIJ,
1433: 0,
1434: 0,
1435: /* 89*/ 0,
1436: 0,
1437: 0,
1438: 0,
1439: 0,
1440: /* 94*/ 0,
1441: 0,
1442: 0,
1443: 0,
1444: 0,
1445: /* 99*/ 0,
1446: 0,
1447: 0,
1448: 0,
1449: 0,
1450: /*104*/ 0,
1451: MatRealPart_SeqSBAIJ,
1452: MatImaginaryPart_SeqSBAIJ,
1453: MatGetRowUpperTriangular_SeqSBAIJ,
1454: MatRestoreRowUpperTriangular_SeqSBAIJ,
1455: /*109*/ 0,
1456: 0,
1457: 0,
1458: 0,
1459: MatMissingDiagonal_SeqSBAIJ,
1460: /*114*/ 0,
1461: 0,
1462: 0,
1463: 0,
1464: 0,
1465: /*119*/ 0,
1466: 0,
1467: 0,
1468: 0,
1469: 0,
1470: /*124*/ 0,
1471: 0,
1472: 0,
1473: 0,
1474: 0,
1475: /*129*/ 0,
1476: 0,
1477: 0,
1478: 0,
1479: 0,
1480: /*134*/ 0,
1481: 0,
1482: 0,
1483: 0,
1484: 0,
1485: /*139*/ MatSetBlockSizes_Default,
1486: 0,
1487: 0,
1488: 0,
1489: 0,
1490: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ
1491: };
1493: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1494: {
1495: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1496: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1500: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1502: /* allocate space for values if not already there */
1503: if (!aij->saved_values) {
1504: PetscMalloc1(nz+1,&aij->saved_values);
1505: }
1507: /* copy values over */
1508: PetscArraycpy(aij->saved_values,aij->a,nz);
1509: return(0);
1510: }
1512: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1513: {
1514: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1516: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1519: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1520: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1522: /* copy values over */
1523: PetscArraycpy(aij->a,aij->saved_values,nz);
1524: return(0);
1525: }
1527: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1528: {
1529: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1531: PetscInt i,mbs,nbs,bs2;
1532: PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1535: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1537: MatSetBlockSize(B,PetscAbs(bs));
1538: PetscLayoutSetUp(B->rmap);
1539: PetscLayoutSetUp(B->cmap);
1540: if (B->rmap->N > B->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"SEQSBAIJ matrix cannot have more rows %D than columns %D",B->rmap->N,B->cmap->N);
1541: PetscLayoutGetBlockSize(B->rmap,&bs);
1543: B->preallocated = PETSC_TRUE;
1545: mbs = B->rmap->N/bs;
1546: nbs = B->cmap->n/bs;
1547: bs2 = bs*bs;
1549: if (mbs*bs != B->rmap->N || nbs*bs!=B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize");
1551: if (nz == MAT_SKIP_ALLOCATION) {
1552: skipallocation = PETSC_TRUE;
1553: nz = 0;
1554: }
1556: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1557: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1558: if (nnz) {
1559: for (i=0; i<mbs; i++) {
1560: 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]);
1561: if (nnz[i] > nbs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D block rowlength %D",i,nnz[i],nbs);
1562: }
1563: }
1565: B->ops->mult = MatMult_SeqSBAIJ_N;
1566: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1567: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1568: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1570: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1571: if (!flg) {
1572: switch (bs) {
1573: case 1:
1574: B->ops->mult = MatMult_SeqSBAIJ_1;
1575: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1576: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1577: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1578: break;
1579: case 2:
1580: B->ops->mult = MatMult_SeqSBAIJ_2;
1581: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1582: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1583: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1584: break;
1585: case 3:
1586: B->ops->mult = MatMult_SeqSBAIJ_3;
1587: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1588: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1589: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1590: break;
1591: case 4:
1592: B->ops->mult = MatMult_SeqSBAIJ_4;
1593: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1594: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1595: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1596: break;
1597: case 5:
1598: B->ops->mult = MatMult_SeqSBAIJ_5;
1599: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1600: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1601: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1602: break;
1603: case 6:
1604: B->ops->mult = MatMult_SeqSBAIJ_6;
1605: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1606: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1607: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1608: break;
1609: case 7:
1610: B->ops->mult = MatMult_SeqSBAIJ_7;
1611: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1612: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1613: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1614: break;
1615: }
1616: }
1618: b->mbs = mbs;
1619: b->nbs = nbs;
1620: if (!skipallocation) {
1621: if (!b->imax) {
1622: PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
1624: b->free_imax_ilen = PETSC_TRUE;
1626: PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));
1627: }
1628: if (!nnz) {
1629: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1630: else if (nz <= 0) nz = 1;
1631: nz = PetscMin(nbs,nz);
1632: for (i=0; i<mbs; i++) b->imax[i] = nz;
1633: nz = nz*mbs; /* total nz */
1634: } else {
1635: PetscInt64 nz64 = 0;
1636: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1637: PetscIntCast(nz64,&nz);
1638: }
1639: /* b->ilen will count nonzeros in each block row so far. */
1640: for (i=0; i<mbs; i++) b->ilen[i] = 0;
1641: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1643: /* allocate the matrix space */
1644: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
1645: PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
1646: PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
1647: PetscArrayzero(b->a,nz*bs2);
1648: PetscArrayzero(b->j,nz);
1650: b->singlemalloc = PETSC_TRUE;
1652: /* pointer to beginning of each row */
1653: b->i[0] = 0;
1654: for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1656: b->free_a = PETSC_TRUE;
1657: b->free_ij = PETSC_TRUE;
1658: } else {
1659: b->free_a = PETSC_FALSE;
1660: b->free_ij = PETSC_FALSE;
1661: }
1663: b->bs2 = bs2;
1664: b->nz = 0;
1665: b->maxnz = nz;
1666: b->inew = 0;
1667: b->jnew = 0;
1668: b->anew = 0;
1669: b->a2anew = 0;
1670: b->permute = PETSC_FALSE;
1672: B->was_assembled = PETSC_FALSE;
1673: B->assembled = PETSC_FALSE;
1674: if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
1675: return(0);
1676: }
1678: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1679: {
1680: PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1681: PetscScalar *values=0;
1682: PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1686: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
1687: PetscLayoutSetBlockSize(B->rmap,bs);
1688: PetscLayoutSetBlockSize(B->cmap,bs);
1689: PetscLayoutSetUp(B->rmap);
1690: PetscLayoutSetUp(B->cmap);
1691: PetscLayoutGetBlockSize(B->rmap,&bs);
1692: m = B->rmap->n/bs;
1694: if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]);
1695: PetscMalloc1(m+1,&nnz);
1696: for (i=0; i<m; i++) {
1697: nz = ii[i+1] - ii[i];
1698: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz);
1699: anz = 0;
1700: for (j=0; j<nz; j++) {
1701: /* count only values on the diagonal or above */
1702: if (jj[ii[i] + j] >= i) {
1703: anz = nz - j;
1704: break;
1705: }
1706: }
1707: nz_max = PetscMax(nz_max,anz);
1708: nnz[i] = anz;
1709: }
1710: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1711: PetscFree(nnz);
1713: values = (PetscScalar*)V;
1714: if (!values) {
1715: PetscCalloc1(bs*bs*nz_max,&values);
1716: }
1717: for (i=0; i<m; i++) {
1718: PetscInt ncols = ii[i+1] - ii[i];
1719: const PetscInt *icols = jj + ii[i];
1720: if (!roworiented || bs == 1) {
1721: const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1722: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
1723: } else {
1724: for (j=0; j<ncols; j++) {
1725: const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1726: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
1727: }
1728: }
1729: }
1730: if (!V) { PetscFree(values); }
1731: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1732: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1733: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1734: return(0);
1735: }
1737: /*
1738: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1739: */
1740: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1741: {
1743: PetscBool flg = PETSC_FALSE;
1744: PetscInt bs = B->rmap->bs;
1747: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1748: if (flg) bs = 8;
1750: if (!natural) {
1751: switch (bs) {
1752: case 1:
1753: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1754: break;
1755: case 2:
1756: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1757: break;
1758: case 3:
1759: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1760: break;
1761: case 4:
1762: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1763: break;
1764: case 5:
1765: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1766: break;
1767: case 6:
1768: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1769: break;
1770: case 7:
1771: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1772: break;
1773: default:
1774: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1775: break;
1776: }
1777: } else {
1778: switch (bs) {
1779: case 1:
1780: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1781: break;
1782: case 2:
1783: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1784: break;
1785: case 3:
1786: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1787: break;
1788: case 4:
1789: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1790: break;
1791: case 5:
1792: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1793: break;
1794: case 6:
1795: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1796: break;
1797: case 7:
1798: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1799: break;
1800: default:
1801: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1802: break;
1803: }
1804: }
1805: return(0);
1806: }
1808: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
1809: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*);
1811: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1812: {
1813: PetscInt n = A->rmap->n;
1817: #if defined(PETSC_USE_COMPLEX)
1818: if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported");
1819: #endif
1820: MatCreate(PetscObjectComm((PetscObject)A),B);
1821: MatSetSizes(*B,n,n,n,n);
1822: if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1823: MatSetType(*B,MATSEQSBAIJ);
1824: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
1826: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1827: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1828: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1830: (*B)->factortype = ftype;
1831: PetscFree((*B)->solvertype);
1832: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
1833: return(0);
1834: }
1836: /*@C
1837: MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1839: Not Collective
1841: Input Parameter:
1842: . mat - a MATSEQSBAIJ matrix
1844: Output Parameter:
1845: . array - pointer to the data
1847: Level: intermediate
1849: .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1850: @*/
1851: PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1852: {
1856: PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
1857: return(0);
1858: }
1860: /*@C
1861: MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1863: Not Collective
1865: Input Parameters:
1866: + mat - a MATSEQSBAIJ matrix
1867: - array - pointer to the data
1869: Level: intermediate
1871: .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1872: @*/
1873: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1874: {
1878: PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
1879: return(0);
1880: }
1882: /*MC
1883: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1884: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1886: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1887: can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd()
1889: Options Database Keys:
1890: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1892: Notes:
1893: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1894: stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1895: the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion.
1897: The number of rows in the matrix must be less than or equal to the number of columns
1899: Level: beginner
1901: .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1902: M*/
1903: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1904: {
1905: Mat_SeqSBAIJ *b;
1907: PetscMPIInt size;
1908: PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1911: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1912: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1914: PetscNewLog(B,&b);
1915: B->data = (void*)b;
1916: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1918: B->ops->destroy = MatDestroy_SeqSBAIJ;
1919: B->ops->view = MatView_SeqSBAIJ;
1920: b->row = 0;
1921: b->icol = 0;
1922: b->reallocs = 0;
1923: b->saved_values = 0;
1924: b->inode.limit = 5;
1925: b->inode.max_limit = 5;
1927: b->roworiented = PETSC_TRUE;
1928: b->nonew = 0;
1929: b->diag = 0;
1930: b->solve_work = 0;
1931: b->mult_work = 0;
1932: B->spptr = 0;
1933: B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1934: b->keepnonzeropattern = PETSC_FALSE;
1936: b->inew = 0;
1937: b->jnew = 0;
1938: b->anew = 0;
1939: b->a2anew = 0;
1940: b->permute = PETSC_FALSE;
1942: b->ignore_ltriangular = PETSC_TRUE;
1944: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);
1946: b->getrow_utriangular = PETSC_FALSE;
1948: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);
1950: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);
1951: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);
1952: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);
1953: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);
1954: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1955: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);
1956: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);
1957: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);
1958: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
1959: #if defined(PETSC_HAVE_ELEMENTAL)
1960: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);
1961: #endif
1963: B->symmetric = PETSC_TRUE;
1964: B->structurally_symmetric = PETSC_TRUE;
1965: B->symmetric_set = PETSC_TRUE;
1966: B->structurally_symmetric_set = PETSC_TRUE;
1967: B->symmetric_eternal = PETSC_TRUE;
1969: B->hermitian = PETSC_FALSE;
1970: B->hermitian_set = PETSC_FALSE;
1972: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
1974: PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");
1975: PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);
1976: if (no_unroll) {
1977: PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");
1978: }
1979: PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);
1980: if (no_inode) {
1981: PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");
1982: }
1983: PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);
1984: PetscOptionsEnd();
1985: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
1986: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
1987: return(0);
1988: }
1990: /*@C
1991: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
1992: compressed row) format. For good matrix assembly performance the
1993: user should preallocate the matrix storage by setting the parameter nz
1994: (or the array nnz). By setting these parameters accurately, performance
1995: during matrix assembly can be increased by more than a factor of 50.
1997: Collective on Mat
1999: Input Parameters:
2000: + B - the symmetric matrix
2001: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2002: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2003: . nz - number of block nonzeros per block row (same for all rows)
2004: - nnz - array containing the number of block nonzeros in the upper triangular plus
2005: diagonal portion of each block (possibly different for each block row) or NULL
2007: Options Database Keys:
2008: + -mat_no_unroll - uses code that does not unroll the loops in the
2009: block calculations (much slower)
2010: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2012: Level: intermediate
2014: Notes:
2015: Specify the preallocated storage with either nz or nnz (not both).
2016: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2017: allocation. See Users-Manual: ch_mat for details.
2019: You can call MatGetInfo() to get information on how effective the preallocation was;
2020: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2021: You can also run with the option -info and look for messages with the string
2022: malloc in them to see if additional memory allocation was needed.
2024: If the nnz parameter is given then the nz parameter is ignored
2027: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2028: @*/
2029: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2030: {
2037: PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
2038: return(0);
2039: }
2041: /*@C
2042: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2044: Input Parameters:
2045: + B - the matrix
2046: . bs - size of block, the blocks are ALWAYS square.
2047: . i - the indices into j for the start of each local row (starts with zero)
2048: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2049: - v - optional values in the matrix
2051: Level: advanced
2053: Notes:
2054: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2055: may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2056: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2057: MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2058: block column and the second index is over columns within a block.
2060: Any entries below the diagonal are ignored
2062: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2063: and usually the numerical values as well
2065: .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2066: @*/
2067: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2068: {
2075: PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
2076: return(0);
2077: }
2079: /*@C
2080: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2081: compressed row) format. For good matrix assembly performance the
2082: user should preallocate the matrix storage by setting the parameter nz
2083: (or the array nnz). By setting these parameters accurately, performance
2084: during matrix assembly can be increased by more than a factor of 50.
2086: Collective
2088: Input Parameters:
2089: + comm - MPI communicator, set to PETSC_COMM_SELF
2090: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2091: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2092: . m - number of rows, or number of columns
2093: . nz - number of block nonzeros per block row (same for all rows)
2094: - nnz - array containing the number of block nonzeros in the upper triangular plus
2095: diagonal portion of each block (possibly different for each block row) or NULL
2097: Output Parameter:
2098: . A - the symmetric matrix
2100: Options Database Keys:
2101: + -mat_no_unroll - uses code that does not unroll the loops in the
2102: block calculations (much slower)
2103: - -mat_block_size - size of the blocks to use
2105: Level: intermediate
2107: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2108: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2109: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2111: Notes:
2112: The number of rows and columns must be divisible by blocksize.
2113: This matrix type does not support complex Hermitian operation.
2115: Specify the preallocated storage with either nz or nnz (not both).
2116: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2117: allocation. See Users-Manual: ch_mat for details.
2119: If the nnz parameter is given then the nz parameter is ignored
2121: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2122: @*/
2123: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2124: {
2128: MatCreate(comm,A);
2129: MatSetSizes(*A,m,n,m,n);
2130: MatSetType(*A,MATSEQSBAIJ);
2131: MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
2132: return(0);
2133: }
2135: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2136: {
2137: Mat C;
2138: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2140: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2143: if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");
2145: *B = 0;
2146: MatCreate(PetscObjectComm((PetscObject)A),&C);
2147: MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
2148: MatSetBlockSizesFromMats(C,A,A);
2149: MatSetType(C,MATSEQSBAIJ);
2150: c = (Mat_SeqSBAIJ*)C->data;
2152: C->preallocated = PETSC_TRUE;
2153: C->factortype = A->factortype;
2154: c->row = 0;
2155: c->icol = 0;
2156: c->saved_values = 0;
2157: c->keepnonzeropattern = a->keepnonzeropattern;
2158: C->assembled = PETSC_TRUE;
2160: PetscLayoutReference(A->rmap,&C->rmap);
2161: PetscLayoutReference(A->cmap,&C->cmap);
2162: c->bs2 = a->bs2;
2163: c->mbs = a->mbs;
2164: c->nbs = a->nbs;
2166: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2167: c->imax = a->imax;
2168: c->ilen = a->ilen;
2169: c->free_imax_ilen = PETSC_FALSE;
2170: } else {
2171: PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);
2172: PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));
2173: for (i=0; i<mbs; i++) {
2174: c->imax[i] = a->imax[i];
2175: c->ilen[i] = a->ilen[i];
2176: }
2177: c->free_imax_ilen = PETSC_TRUE;
2178: }
2180: /* allocate the matrix space */
2181: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2182: PetscMalloc1(bs2*nz,&c->a);
2183: PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));
2184: c->i = a->i;
2185: c->j = a->j;
2186: c->singlemalloc = PETSC_FALSE;
2187: c->free_a = PETSC_TRUE;
2188: c->free_ij = PETSC_FALSE;
2189: c->parent = A;
2190: PetscObjectReference((PetscObject)A);
2191: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2192: MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2193: } else {
2194: PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
2195: PetscArraycpy(c->i,a->i,mbs+1);
2196: PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
2197: c->singlemalloc = PETSC_TRUE;
2198: c->free_a = PETSC_TRUE;
2199: c->free_ij = PETSC_TRUE;
2200: }
2201: if (mbs > 0) {
2202: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2203: PetscArraycpy(c->j,a->j,nz);
2204: }
2205: if (cpvalues == MAT_COPY_VALUES) {
2206: PetscArraycpy(c->a,a->a,bs2*nz);
2207: } else {
2208: PetscArrayzero(c->a,bs2*nz);
2209: }
2210: if (a->jshort) {
2211: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2212: /* if the parent matrix is reassembled, this child matrix will never notice */
2213: PetscMalloc1(nz,&c->jshort);
2214: PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));
2215: PetscArraycpy(c->jshort,a->jshort,nz);
2217: c->free_jshort = PETSC_TRUE;
2218: }
2219: }
2221: c->roworiented = a->roworiented;
2222: c->nonew = a->nonew;
2224: if (a->diag) {
2225: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2226: c->diag = a->diag;
2227: c->free_diag = PETSC_FALSE;
2228: } else {
2229: PetscMalloc1(mbs,&c->diag);
2230: PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));
2231: for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2232: c->free_diag = PETSC_TRUE;
2233: }
2234: }
2235: c->nz = a->nz;
2236: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2237: c->solve_work = 0;
2238: c->mult_work = 0;
2240: *B = C;
2241: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
2242: return(0);
2243: }
2245: PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer)
2246: {
2247: Mat_SeqSBAIJ *a;
2249: int fd;
2250: PetscMPIInt size;
2251: PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs;
2252: PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount;
2253: PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols;
2254: PetscInt *masked,nmask,tmp,bs2,ishift;
2255: PetscScalar *aa;
2256: MPI_Comm comm;
2257: PetscBool isbinary;
2260: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2261: if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)newmat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newmat)->type_name);
2263: /* force binary viewer to load .info file if it has not yet done so */
2264: PetscViewerSetUp(viewer);
2265: PetscObjectGetComm((PetscObject)viewer,&comm);
2266: PetscOptionsGetInt(((PetscObject)newmat)->options,((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);
2267: if (bs < 0) bs = 1;
2268: bs2 = bs*bs;
2270: MPI_Comm_size(comm,&size);
2271: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
2272: PetscViewerBinaryGetDescriptor(viewer,&fd);
2273: PetscBinaryRead(fd,header,4,NULL,PETSC_INT);
2274: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object");
2275: M = header[1]; N = header[2]; nz = header[3];
2277: if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ");
2279: if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices");
2281: /*
2282: This code adds extra rows to make sure the number of rows is
2283: divisible by the blocksize
2284: */
2285: mbs = M/bs;
2286: extra_rows = bs - M + bs*(mbs);
2287: if (extra_rows == bs) extra_rows = 0;
2288: else mbs++;
2289: if (extra_rows) {
2290: PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");
2291: }
2293: /* Set global sizes if not already set */
2294: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
2295: MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);
2296: } else { /* Check if the matrix global sizes are correct */
2297: MatGetSize(newmat,&rows,&cols);
2298: if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);
2299: }
2301: /* read in row lengths */
2302: PetscMalloc1(M+extra_rows,&rowlengths);
2303: PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);
2304: for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;
2306: /* read in column indices */
2307: PetscMalloc1(nz+extra_rows,&jj);
2308: PetscBinaryRead(fd,jj,nz,NULL,PETSC_INT);
2309: for (i=0; i<extra_rows; i++) jj[nz+i] = M+i;
2311: /* loop over row lengths determining block row lengths */
2312: PetscCalloc2(mbs,&s_browlengths,mbs,&mask);
2313: PetscMalloc1(mbs,&masked);
2314: rowcount = 0;
2315: nzcount = 0;
2316: for (i=0; i<mbs; i++) {
2317: nmask = 0;
2318: for (j=0; j<bs; j++) {
2319: kmax = rowlengths[rowcount];
2320: for (k=0; k<kmax; k++) {
2321: tmp = jj[nzcount++]/bs; /* block col. index */
2322: if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;}
2323: }
2324: rowcount++;
2325: }
2326: s_browlengths[i] += nmask;
2328: /* zero out the mask elements we set */
2329: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
2330: }
2332: /* Do preallocation */
2333: MatSeqSBAIJSetPreallocation(newmat,bs,0,s_browlengths);
2334: a = (Mat_SeqSBAIJ*)newmat->data;
2336: /* set matrix "i" values */
2337: a->i[0] = 0;
2338: for (i=1; i<= mbs; i++) {
2339: a->i[i] = a->i[i-1] + s_browlengths[i-1];
2340: a->ilen[i-1] = s_browlengths[i-1];
2341: }
2342: a->nz = a->i[mbs];
2344: /* read in nonzero values */
2345: PetscMalloc1(nz+extra_rows,&aa);
2346: PetscBinaryRead(fd,aa,nz,NULL,PETSC_SCALAR);
2347: for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0;
2349: /* set "a" and "j" values into matrix */
2350: nzcount = 0; jcount = 0;
2351: for (i=0; i<mbs; i++) {
2352: nzcountb = nzcount;
2353: nmask = 0;
2354: for (j=0; j<bs; j++) {
2355: kmax = rowlengths[i*bs+j];
2356: for (k=0; k<kmax; k++) {
2357: tmp = jj[nzcount++]/bs; /* block col. index */
2358: if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;}
2359: }
2360: }
2361: /* sort the masked values */
2362: PetscSortInt(nmask,masked);
2364: /* set "j" values into matrix */
2365: maskcount = 1;
2366: for (j=0; j<nmask; j++) {
2367: a->j[jcount++] = masked[j];
2368: mask[masked[j]] = maskcount++;
2369: }
2371: /* set "a" values into matrix */
2372: ishift = bs2*a->i[i];
2373: for (j=0; j<bs; j++) {
2374: kmax = rowlengths[i*bs+j];
2375: for (k=0; k<kmax; k++) {
2376: tmp = jj[nzcountb]/bs; /* block col. index */
2377: if (tmp >= i) {
2378: block = mask[tmp] - 1;
2379: point = jj[nzcountb] - bs*tmp;
2380: idx = ishift + bs2*block + j + bs*point;
2381: a->a[idx] = aa[nzcountb];
2382: }
2383: nzcountb++;
2384: }
2385: }
2386: /* zero out the mask elements we set */
2387: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
2388: }
2389: if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix");
2391: PetscFree(rowlengths);
2392: PetscFree2(s_browlengths,mask);
2393: PetscFree(aa);
2394: PetscFree(jj);
2395: PetscFree(masked);
2397: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2398: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
2399: return(0);
2400: }
2402: /*@
2403: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2404: (upper triangular entries in CSR format) provided by the user.
2406: Collective
2408: Input Parameters:
2409: + comm - must be an MPI communicator of size 1
2410: . bs - size of block
2411: . m - number of rows
2412: . n - number of columns
2413: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix
2414: . j - column indices
2415: - a - matrix values
2417: Output Parameter:
2418: . mat - the matrix
2420: Level: advanced
2422: Notes:
2423: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2424: once the matrix is destroyed
2426: You cannot set new nonzero locations into this matrix, that will generate an error.
2428: The i and j indices are 0 based
2430: When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1
2431: it is the regular CSR format excluding the lower triangular elements.
2433: .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2435: @*/
2436: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2437: {
2439: PetscInt ii;
2440: Mat_SeqSBAIJ *sbaij;
2443: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2444: if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2446: MatCreate(comm,mat);
2447: MatSetSizes(*mat,m,n,m,n);
2448: MatSetType(*mat,MATSEQSBAIJ);
2449: MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,0);
2450: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2451: PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);
2452: PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));
2454: sbaij->i = i;
2455: sbaij->j = j;
2456: sbaij->a = a;
2458: sbaij->singlemalloc = PETSC_FALSE;
2459: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2460: sbaij->free_a = PETSC_FALSE;
2461: sbaij->free_ij = PETSC_FALSE;
2462: sbaij->free_imax_ilen = PETSC_TRUE;
2464: for (ii=0; ii<m; ii++) {
2465: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2466: #if defined(PETSC_USE_DEBUG)
2467: if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]);
2468: #endif
2469: }
2470: #if defined(PETSC_USE_DEBUG)
2471: for (ii=0; ii<sbaij->i[m]; ii++) {
2472: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2473: if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]);
2474: }
2475: #endif
2477: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2478: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2479: return(0);
2480: }
2482: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2483: {
2485: PetscMPIInt size;
2488: MPI_Comm_size(comm,&size);
2489: if (size == 1 && scall == MAT_REUSE_MATRIX) {
2490: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2491: } else {
2492: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);
2493: }
2494: return(0);
2495: }