Actual source code: sbaij.c
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: #if defined(PETSC_HAVE_SCALAPACK)
18: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
19: #endif
20: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat,MatType,MatReuse,Mat*);
22: /*
23: Checks for missing diagonals
24: */
25: PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd)
26: {
27: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
28: PetscInt *diag,*ii = a->i,i;
30: MatMarkDiagonal_SeqSBAIJ(A);
31: *missing = PETSC_FALSE;
32: if (A->rmap->n > 0 && !ii) {
33: *missing = PETSC_TRUE;
34: if (dd) *dd = 0;
35: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
36: } else {
37: diag = a->diag;
38: for (i=0; i<a->mbs; i++) {
39: if (diag[i] >= ii[i+1]) {
40: *missing = PETSC_TRUE;
41: if (dd) *dd = i;
42: break;
43: }
44: }
45: }
46: return 0;
47: }
49: PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
50: {
51: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
52: PetscInt i,j;
54: if (!a->diag) {
55: PetscMalloc1(a->mbs,&a->diag);
56: PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt));
57: a->free_diag = PETSC_TRUE;
58: }
59: for (i=0; i<a->mbs; i++) {
60: a->diag[i] = a->i[i+1];
61: for (j=a->i[i]; j<a->i[i+1]; j++) {
62: if (a->j[j] == i) {
63: a->diag[i] = j;
64: break;
65: }
66: }
67: }
68: return 0;
69: }
71: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done)
72: {
73: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
74: PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt;
75: PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja;
77: *nn = n;
78: if (!ia) return 0;
79: if (symmetric) {
80: MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja);
81: nz = tia[n];
82: } else {
83: tia = a->i; tja = a->j;
84: }
86: if (!blockcompressed && bs > 1) {
87: (*nn) *= bs;
88: /* malloc & create the natural set of indices */
89: PetscMalloc1((n+1)*bs,ia);
90: if (n) {
91: (*ia)[0] = oshift;
92: for (j=1; j<bs; j++) {
93: (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1];
94: }
95: }
97: for (i=1; i<n; i++) {
98: (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1];
99: for (j=1; j<bs; j++) {
100: (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1];
101: }
102: }
103: if (n) {
104: (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1];
105: }
107: if (inja) {
108: PetscMalloc1(nz*bs*bs,ja);
109: cnt = 0;
110: for (i=0; i<n; i++) {
111: for (j=0; j<bs; j++) {
112: for (k=tia[i]; k<tia[i+1]; k++) {
113: for (l=0; l<bs; l++) {
114: (*ja)[cnt++] = bs*tja[k] + l;
115: }
116: }
117: }
118: }
119: }
121: if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
122: PetscFree(tia);
123: PetscFree(tja);
124: }
125: } else if (oshift == 1) {
126: if (symmetric) {
127: nz = tia[A->rmap->n/bs];
128: /* add 1 to i and j indices */
129: for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1;
130: *ia = tia;
131: if (ja) {
132: for (i=0; i<nz; i++) tja[i] = tja[i] + 1;
133: *ja = tja;
134: }
135: } else {
136: nz = a->i[A->rmap->n/bs];
137: /* malloc space and add 1 to i and j indices */
138: PetscMalloc1(A->rmap->n/bs+1,ia);
139: for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1;
140: if (ja) {
141: PetscMalloc1(nz,ja);
142: for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1;
143: }
144: }
145: } else {
146: *ia = tia;
147: if (ja) *ja = tja;
148: }
149: return 0;
150: }
152: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
153: {
154: if (!ia) return 0;
155: if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
156: PetscFree(*ia);
157: if (ja) PetscFree(*ja);
158: }
159: return 0;
160: }
162: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
163: {
164: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
166: #if defined(PETSC_USE_LOG)
167: PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->N,a->nz);
168: #endif
169: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
170: if (a->free_diag) PetscFree(a->diag);
171: ISDestroy(&a->row);
172: ISDestroy(&a->col);
173: ISDestroy(&a->icol);
174: PetscFree(a->idiag);
175: PetscFree(a->inode.size);
176: if (a->free_imax_ilen) PetscFree2(a->imax,a->ilen);
177: PetscFree(a->solve_work);
178: PetscFree(a->sor_work);
179: PetscFree(a->solves_work);
180: PetscFree(a->mult_work);
181: PetscFree(a->saved_values);
182: if (a->free_jshort) PetscFree(a->jshort);
183: PetscFree(a->inew);
184: MatDestroy(&a->parent);
185: PetscFree(A->data);
187: PetscObjectChangeTypeName((PetscObject)A,NULL);
188: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
189: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
190: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL);
191: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL);
192: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL);
193: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL);
194: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL);
195: #if defined(PETSC_HAVE_ELEMENTAL)
196: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL);
197: #endif
198: #if defined(PETSC_HAVE_SCALAPACK)
199: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_scalapack_C",NULL);
200: #endif
201: return 0;
202: }
204: PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg)
205: {
206: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
207: #if defined(PETSC_USE_COMPLEX)
208: PetscInt bs;
209: #endif
211: #if defined(PETSC_USE_COMPLEX)
212: MatGetBlockSize(A,&bs);
213: #endif
214: switch (op) {
215: case MAT_ROW_ORIENTED:
216: a->roworiented = flg;
217: break;
218: case MAT_KEEP_NONZERO_PATTERN:
219: a->keepnonzeropattern = flg;
220: break;
221: case MAT_NEW_NONZERO_LOCATIONS:
222: a->nonew = (flg ? 0 : 1);
223: break;
224: case MAT_NEW_NONZERO_LOCATION_ERR:
225: a->nonew = (flg ? -1 : 0);
226: break;
227: case MAT_NEW_NONZERO_ALLOCATION_ERR:
228: a->nonew = (flg ? -2 : 0);
229: break;
230: case MAT_UNUSED_NONZERO_LOCATION_ERR:
231: a->nounused = (flg ? -1 : 0);
232: break;
233: case MAT_FORCE_DIAGONAL_ENTRIES:
234: case MAT_IGNORE_OFF_PROC_ENTRIES:
235: case MAT_USE_HASH_TABLE:
236: case MAT_SORTED_FULL:
237: PetscInfo(A,"Option %s ignored\n",MatOptions[op]);
238: break;
239: case MAT_HERMITIAN:
240: #if defined(PETSC_USE_COMPLEX)
241: if (flg) { /* disable transpose ops */
243: A->ops->multtranspose = NULL;
244: A->ops->multtransposeadd = NULL;
245: A->symmetric = PETSC_FALSE;
246: }
247: #endif
248: break;
249: case MAT_SYMMETRIC:
250: case MAT_SPD:
251: #if defined(PETSC_USE_COMPLEX)
252: if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */
253: A->ops->multtranspose = A->ops->mult;
254: A->ops->multtransposeadd = A->ops->multadd;
255: }
256: #endif
257: break;
258: /* These options are handled directly by MatSetOption() */
259: case MAT_STRUCTURALLY_SYMMETRIC:
260: case MAT_SYMMETRY_ETERNAL:
261: case MAT_STRUCTURE_ONLY:
262: /* These options are handled directly by MatSetOption() */
263: break;
264: case MAT_IGNORE_LOWER_TRIANGULAR:
265: a->ignore_ltriangular = flg;
266: break;
267: case MAT_ERROR_LOWER_TRIANGULAR:
268: a->ignore_ltriangular = flg;
269: break;
270: case MAT_GETROW_UPPERTRIANGULAR:
271: a->getrow_utriangular = flg;
272: break;
273: case MAT_SUBMAT_SINGLEIS:
274: break;
275: default:
276: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
277: }
278: return 0;
279: }
281: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
282: {
283: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
287: /* Get the upper triangular part of the row */
288: MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);
289: return 0;
290: }
292: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
293: {
294: if (nz) *nz = 0;
295: if (idx) PetscFree(*idx);
296: if (v) PetscFree(*v);
297: return 0;
298: }
300: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
301: {
302: 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;
312: a->getrow_utriangular = PETSC_FALSE;
313: return 0;
314: }
316: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B)
317: {
318: if (reuse == MAT_INITIAL_MATRIX) {
319: MatDuplicate(A,MAT_COPY_VALUES,B);
320: } else if (reuse == MAT_REUSE_MATRIX) {
321: MatCopy(A,*B,SAME_NONZERO_PATTERN);
322: }
323: return 0;
324: }
326: PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
327: {
328: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
329: PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
330: PetscViewerFormat format;
331: PetscInt *diag;
333: PetscViewerGetFormat(viewer,&format);
334: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
335: PetscViewerASCIIPrintf(viewer," block size is %" PetscInt_FMT "\n",bs);
336: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
337: Mat aij;
338: const char *matname;
340: if (A->factortype && bs>1) {
341: PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");
342: return 0;
343: }
344: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);
345: PetscObjectGetName((PetscObject)A,&matname);
346: PetscObjectSetName((PetscObject)aij,matname);
347: MatView(aij,viewer);
348: MatDestroy(&aij);
349: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
350: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
351: for (i=0; i<a->mbs; i++) {
352: for (j=0; j<bs; j++) {
353: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j);
354: for (k=a->i[i]; k<a->i[i+1]; k++) {
355: for (l=0; l<bs; l++) {
356: #if defined(PETSC_USE_COMPLEX)
357: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
358: PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l,
359: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j])));
360: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
361: PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l,
362: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j])));
363: } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
364: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
365: }
366: #else
367: if (a->a[bs2*k + l*bs + j] != 0.0) {
368: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
369: }
370: #endif
371: }
372: }
373: PetscViewerASCIIPrintf(viewer,"\n");
374: }
375: }
376: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
377: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
378: return 0;
379: } else {
380: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
381: if (A->factortype) { /* for factored matrix */
384: diag=a->diag;
385: for (i=0; i<a->mbs; i++) { /* for row block i */
386: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
387: /* diagonal entry */
388: #if defined(PETSC_USE_COMPLEX)
389: if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
390: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),(double)PetscImaginaryPart(1.0/a->a[diag[i]]));
391: } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
392: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),-(double)PetscImaginaryPart(1.0/a->a[diag[i]]));
393: } else {
394: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));
395: }
396: #else
397: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));
398: #endif
399: /* off-diagonal entries */
400: for (k=a->i[i]; k<a->i[i+1]-1; k++) {
401: #if defined(PETSC_USE_COMPLEX)
402: if (PetscImaginaryPart(a->a[k]) > 0.0) {
403: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));
404: } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
405: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));
406: } else {
407: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));
408: }
409: #else
410: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[k],(double)a->a[k]);
411: #endif
412: }
413: PetscViewerASCIIPrintf(viewer,"\n");
414: }
416: } else { /* for non-factored matrix */
417: for (i=0; i<a->mbs; i++) { /* for row block i */
418: for (j=0; j<bs; j++) { /* for row bs*i + j */
419: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j);
420: for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */
421: for (l=0; l<bs; l++) { /* for column */
422: #if defined(PETSC_USE_COMPLEX)
423: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
424: PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l,
425: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j])));
426: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
427: PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l,
428: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j])));
429: } else {
430: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
431: }
432: #else
433: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
434: #endif
435: }
436: }
437: PetscViewerASCIIPrintf(viewer,"\n");
438: }
439: }
440: }
441: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
442: }
443: PetscViewerFlush(viewer);
444: return 0;
445: }
447: #include <petscdraw.h>
448: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
449: {
450: Mat A = (Mat) Aa;
451: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
452: PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
453: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
454: MatScalar *aa;
455: PetscViewer viewer;
458: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
459: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
461: /* loop over matrix elements drawing boxes */
463: PetscDrawCollectiveBegin(draw);
464: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
465: /* Blue for negative, Cyan for zero and Red for positive */
466: color = PETSC_DRAW_BLUE;
467: for (i=0,row=0; i<mbs; i++,row+=bs) {
468: for (j=a->i[i]; j<a->i[i+1]; j++) {
469: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
470: x_l = a->j[j]*bs; x_r = x_l + 1.0;
471: aa = a->a + j*bs2;
472: for (k=0; k<bs; k++) {
473: for (l=0; l<bs; l++) {
474: if (PetscRealPart(*aa++) >= 0.) continue;
475: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
476: }
477: }
478: }
479: }
480: color = PETSC_DRAW_CYAN;
481: for (i=0,row=0; i<mbs; i++,row+=bs) {
482: for (j=a->i[i]; j<a->i[i+1]; j++) {
483: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
484: x_l = a->j[j]*bs; x_r = x_l + 1.0;
485: aa = a->a + j*bs2;
486: for (k=0; k<bs; k++) {
487: for (l=0; l<bs; l++) {
488: if (PetscRealPart(*aa++) != 0.) continue;
489: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
490: }
491: }
492: }
493: }
494: color = PETSC_DRAW_RED;
495: for (i=0,row=0; i<mbs; i++,row+=bs) {
496: for (j=a->i[i]; j<a->i[i+1]; j++) {
497: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
498: x_l = a->j[j]*bs; x_r = x_l + 1.0;
499: aa = a->a + j*bs2;
500: for (k=0; k<bs; k++) {
501: for (l=0; l<bs; l++) {
502: if (PetscRealPart(*aa++) <= 0.) continue;
503: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
504: }
505: }
506: }
507: }
508: PetscDrawCollectiveEnd(draw);
509: return 0;
510: }
512: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
513: {
514: PetscReal xl,yl,xr,yr,w,h;
515: PetscDraw draw;
516: PetscBool isnull;
518: PetscViewerDrawGetDraw(viewer,0,&draw);
519: PetscDrawIsNull(draw,&isnull);
520: if (isnull) return 0;
522: xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
523: xr += w; yr += h; xl = -w; yl = -h;
524: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
525: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
526: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
527: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
528: PetscDrawSave(draw);
529: return 0;
530: }
532: /* Used for both MPIBAIJ and MPISBAIJ matrices */
533: #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary
535: PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
536: {
537: PetscBool iascii,isbinary,isdraw;
539: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
540: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
541: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
542: if (iascii) {
543: MatView_SeqSBAIJ_ASCII(A,viewer);
544: } else if (isbinary) {
545: MatView_SeqSBAIJ_Binary(A,viewer);
546: } else if (isdraw) {
547: MatView_SeqSBAIJ_Draw(A,viewer);
548: } else {
549: Mat B;
550: const char *matname;
551: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
552: PetscObjectGetName((PetscObject)A,&matname);
553: PetscObjectSetName((PetscObject)B,matname);
554: MatView(B,viewer);
555: MatDestroy(&B);
556: }
557: return 0;
558: }
560: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
561: {
562: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
563: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
564: PetscInt *ai = a->i,*ailen = a->ilen;
565: PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
566: MatScalar *ap,*aa = a->a;
568: for (k=0; k<m; k++) { /* loop over rows */
569: row = im[k]; brow = row/bs;
570: if (row < 0) {v += n; continue;} /* negative row */
572: rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
573: nrow = ailen[brow];
574: for (l=0; l<n; l++) { /* loop over columns */
575: if (in[l] < 0) {v++; continue;} /* negative column */
577: col = in[l];
578: bcol = col/bs;
579: cidx = col%bs;
580: ridx = row%bs;
581: high = nrow;
582: low = 0; /* assume unsorted */
583: while (high-low > 5) {
584: t = (low+high)/2;
585: if (rp[t] > bcol) high = t;
586: else low = t;
587: }
588: for (i=low; i<high; i++) {
589: if (rp[i] > bcol) break;
590: if (rp[i] == bcol) {
591: *v++ = ap[bs2*i+bs*cidx+ridx];
592: goto finished;
593: }
594: }
595: *v++ = 0.0;
596: finished:;
597: }
598: }
599: return 0;
600: }
602: PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B)
603: {
604: Mat C;
606: MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C);
607: MatPermute(C,rowp,colp,B);
608: MatDestroy(&C);
609: if (rowp == colp) {
610: MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B);
611: }
612: return 0;
613: }
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: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
619: PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
620: PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
621: PetscBool roworiented=a->roworiented;
622: const PetscScalar *value = v;
623: MatScalar *ap,*aa = a->a,*bap;
625: if (roworiented) stepval = (n-1)*bs;
626: else stepval = (m-1)*bs;
628: for (k=0; k<m; k++) { /* loop over added rows */
629: row = im[k];
630: if (row < 0) continue;
632: rp = aj + ai[row];
633: ap = aa + bs2*ai[row];
634: rmax = imax[row];
635: nrow = ailen[row];
636: low = 0;
637: high = nrow;
638: for (l=0; l<n; l++) { /* loop over added columns */
639: if (in[l] < 0) continue;
640: col = in[l];
642: if (col < row) {
643: if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
644: 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)");
645: }
646: if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
647: else value = v + l*(stepval+bs)*bs + k*bs;
649: if (col <= lastcol) low = 0;
650: else high = nrow;
652: lastcol = col;
653: while (high-low > 7) {
654: t = (low+high)/2;
655: if (rp[t] > col) high = t;
656: else low = t;
657: }
658: for (i=low; i<high; i++) {
659: if (rp[i] > col) break;
660: if (rp[i] == col) {
661: bap = ap + bs2*i;
662: if (roworiented) {
663: if (is == ADD_VALUES) {
664: for (ii=0; ii<bs; ii++,value+=stepval) {
665: for (jj=ii; jj<bs2; jj+=bs) {
666: bap[jj] += *value++;
667: }
668: }
669: } else {
670: for (ii=0; ii<bs; ii++,value+=stepval) {
671: for (jj=ii; jj<bs2; jj+=bs) {
672: bap[jj] = *value++;
673: }
674: }
675: }
676: } else {
677: if (is == ADD_VALUES) {
678: for (ii=0; ii<bs; ii++,value+=stepval) {
679: for (jj=0; jj<bs; jj++) {
680: *bap++ += *value++;
681: }
682: }
683: } else {
684: for (ii=0; ii<bs; ii++,value+=stepval) {
685: for (jj=0; jj<bs; jj++) {
686: *bap++ = *value++;
687: }
688: }
689: }
690: }
691: goto noinsert2;
692: }
693: }
694: if (nonew == 1) goto noinsert2;
696: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
697: N = nrow++ - 1; high++;
698: /* shift up all the later entries in this row */
699: PetscArraymove(rp+i+1,rp+i,N-i+1);
700: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
701: PetscArrayzero(ap+bs2*i,bs2);
702: rp[i] = col;
703: bap = ap + bs2*i;
704: if (roworiented) {
705: for (ii=0; ii<bs; ii++,value+=stepval) {
706: for (jj=ii; jj<bs2; jj+=bs) {
707: bap[jj] = *value++;
708: }
709: }
710: } else {
711: for (ii=0; ii<bs; ii++,value+=stepval) {
712: for (jj=0; jj<bs; jj++) {
713: *bap++ = *value++;
714: }
715: }
716: }
717: noinsert2:;
718: low = i;
719: }
720: ailen[row] = nrow;
721: }
722: return 0;
723: }
725: /*
726: This is not yet used
727: */
728: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
729: {
730: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
731: const PetscInt *ai = a->i, *aj = a->j,*cols;
732: PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
733: PetscBool flag;
735: PetscMalloc1(m,&ns);
736: while (i < m) {
737: nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
738: /* Limits the number of elements in a node to 'a->inode.limit' */
739: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
740: nzy = ai[j+1] - ai[j];
741: if (nzy != (nzx - j + i)) break;
742: PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);
743: if (!flag) break;
744: }
745: ns[node_count++] = blk_size;
747: i = j;
748: }
749: if (!a->inode.size && m && node_count > .9*m) {
750: PetscFree(ns);
751: PetscInfo(A,"Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n",node_count,m);
752: } else {
753: a->inode.node_count = node_count;
755: PetscMalloc1(node_count,&a->inode.size);
756: PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));
757: PetscArraycpy(a->inode.size,ns,node_count);
758: PetscFree(ns);
759: PetscInfo(A,"Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n",node_count,m,a->inode.limit);
761: /* count collections of adjacent columns in each inode */
762: row = 0;
763: cnt = 0;
764: for (i=0; i<node_count; i++) {
765: cols = aj + ai[row] + a->inode.size[i];
766: nz = ai[row+1] - ai[row] - a->inode.size[i];
767: for (j=1; j<nz; j++) {
768: if (cols[j] != cols[j-1]+1) cnt++;
769: }
770: cnt++;
771: row += a->inode.size[i];
772: }
773: PetscMalloc1(2*cnt,&counts);
774: cnt = 0;
775: row = 0;
776: for (i=0; i<node_count; i++) {
777: cols = aj + ai[row] + a->inode.size[i];
778: counts[2*cnt] = cols[0];
779: nz = ai[row+1] - ai[row] - a->inode.size[i];
780: cnt2 = 1;
781: for (j=1; j<nz; j++) {
782: if (cols[j] != cols[j-1]+1) {
783: counts[2*(cnt++)+1] = cnt2;
784: counts[2*cnt] = cols[j];
785: cnt2 = 1;
786: } else cnt2++;
787: }
788: counts[2*(cnt++)+1] = cnt2;
789: row += a->inode.size[i];
790: }
791: PetscIntView(2*cnt,counts,NULL);
792: }
793: return 0;
794: }
796: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
797: {
798: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
799: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
800: PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
801: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
802: MatScalar *aa = a->a,*ap;
804: if (mode == MAT_FLUSH_ASSEMBLY) return 0;
806: if (m) rmax = ailen[0];
807: for (i=1; i<mbs; i++) {
808: /* move each row back by the amount of empty slots (fshift) before it*/
809: fshift += imax[i-1] - ailen[i-1];
810: rmax = PetscMax(rmax,ailen[i]);
811: if (fshift) {
812: ip = aj + ai[i];
813: ap = aa + bs2*ai[i];
814: N = ailen[i];
815: PetscArraymove(ip-fshift,ip,N);
816: PetscArraymove(ap-bs2*fshift,ap,bs2*N);
817: }
818: ai[i] = ai[i-1] + ailen[i-1];
819: }
820: if (mbs) {
821: fshift += imax[mbs-1] - ailen[mbs-1];
822: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
823: }
824: /* reset ilen and imax for each row */
825: for (i=0; i<mbs; i++) {
826: ailen[i] = imax[i] = ai[i+1] - ai[i];
827: }
828: a->nz = ai[mbs];
830: /* diagonals may have moved, reset it */
831: if (a->diag) {
832: PetscArraycpy(a->diag,ai,mbs);
833: }
836: PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);
837: PetscInfo(A,"Number of mallocs during MatSetValues is %" PetscInt_FMT "\n",a->reallocs);
838: PetscInfo(A,"Most nonzeros blocks in any row is %" PetscInt_FMT "\n",rmax);
840: A->info.mallocs += a->reallocs;
841: a->reallocs = 0;
842: A->info.nz_unneeded = (PetscReal)fshift*bs2;
843: a->idiagvalid = PETSC_FALSE;
844: a->rmax = rmax;
846: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
847: if (a->jshort && a->free_jshort) {
848: /* when matrix data structure is changed, previous jshort must be replaced */
849: PetscFree(a->jshort);
850: }
851: PetscMalloc1(a->i[A->rmap->n],&a->jshort);
852: PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));
853: for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
854: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
855: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
856: a->free_jshort = PETSC_TRUE;
857: }
858: return 0;
859: }
861: /*
862: This function returns an array of flags which indicate the locations of contiguous
863: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
864: then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)]
865: Assume: sizes should be long enough to hold all the values.
866: */
867: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
868: {
869: PetscInt i,j,k,row;
870: PetscBool flg;
872: for (i=0,j=0; i<n; j++) {
873: row = idx[i];
874: if (row%bs!=0) { /* Not the beginning of a block */
875: sizes[j] = 1;
876: i++;
877: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
878: sizes[j] = 1; /* Also makes sure at least 'bs' values exist for next else */
879: i++;
880: } else { /* Beginning of the block, so check if the complete block exists */
881: flg = PETSC_TRUE;
882: for (k=1; k<bs; k++) {
883: if (row+k != idx[i+k]) { /* break in the block */
884: flg = PETSC_FALSE;
885: break;
886: }
887: }
888: if (flg) { /* No break in the bs */
889: sizes[j] = bs;
890: i += bs;
891: } else {
892: sizes[j] = 1;
893: i++;
894: }
895: }
896: }
897: *bs_max = j;
898: return 0;
899: }
901: /* Only add/insert a(i,j) with i<=j (blocks).
902: Any a(i,j) with i>j input by user is ingored.
903: */
905: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
906: {
907: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
908: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
909: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
910: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
911: PetscInt ridx,cidx,bs2=a->bs2;
912: MatScalar *ap,value,*aa=a->a,*bap;
914: for (k=0; k<m; k++) { /* loop over added rows */
915: row = im[k]; /* row number */
916: brow = row/bs; /* block row number */
917: if (row < 0) continue;
919: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
920: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
921: rmax = imax[brow]; /* maximum space allocated for this row */
922: nrow = ailen[brow]; /* actual length of this row */
923: low = 0;
924: high = nrow;
925: for (l=0; l<n; l++) { /* loop over added columns */
926: if (in[l] < 0) continue;
928: col = in[l];
929: bcol = col/bs; /* block col number */
931: if (brow > bcol) {
932: if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
933: 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)");
934: }
936: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
937: if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
938: /* element value a(k,l) */
939: if (roworiented) value = v[l + k*n];
940: else value = v[k + l*m];
942: /* move pointer bap to a(k,l) quickly and add/insert value */
943: if (col <= lastcol) low = 0;
944: else high = nrow;
946: lastcol = col;
947: while (high-low > 7) {
948: t = (low+high)/2;
949: if (rp[t] > bcol) high = t;
950: else low = t;
951: }
952: for (i=low; i<high; i++) {
953: if (rp[i] > bcol) break;
954: if (rp[i] == bcol) {
955: bap = ap + bs2*i + bs*cidx + ridx;
956: if (is == ADD_VALUES) *bap += value;
957: else *bap = value;
958: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
959: if (brow == bcol && ridx < cidx) {
960: bap = ap + bs2*i + bs*ridx + cidx;
961: if (is == ADD_VALUES) *bap += value;
962: else *bap = value;
963: }
964: goto noinsert1;
965: }
966: }
968: if (nonew == 1) goto noinsert1;
970: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
972: N = nrow++ - 1; high++;
973: /* shift up all the later entries in this row */
974: PetscArraymove(rp+i+1,rp+i,N-i+1);
975: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
976: PetscArrayzero(ap+bs2*i,bs2);
977: rp[i] = bcol;
978: ap[bs2*i + bs*cidx + ridx] = value;
979: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
980: if (brow == bcol && ridx < cidx) {
981: ap[bs2*i + bs*ridx + cidx] = value;
982: }
983: A->nonzerostate++;
984: noinsert1:;
985: low = i;
986: }
987: } /* end of loop over added columns */
988: ailen[brow] = nrow;
989: } /* end of loop over added rows */
990: return 0;
991: }
993: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
994: {
995: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
996: Mat outA;
997: PetscBool row_identity;
1000: ISIdentity(row,&row_identity);
1004: outA = inA;
1005: inA->factortype = MAT_FACTOR_ICC;
1006: PetscFree(inA->solvertype);
1007: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
1009: MatMarkDiagonal_SeqSBAIJ(inA);
1010: MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);
1012: PetscObjectReference((PetscObject)row);
1013: ISDestroy(&a->row);
1014: a->row = row;
1015: PetscObjectReference((PetscObject)row);
1016: ISDestroy(&a->col);
1017: a->col = row;
1019: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1020: if (a->icol) ISInvertPermutation(row,PETSC_DECIDE, &a->icol);
1021: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
1023: if (!a->solve_work) {
1024: PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
1025: PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
1026: }
1028: MatCholeskyFactorNumeric(outA,inA,info);
1029: return 0;
1030: }
1032: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1033: {
1034: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1035: PetscInt i,nz,n;
1037: nz = baij->maxnz;
1038: n = mat->cmap->n;
1039: for (i=0; i<nz; i++) baij->j[i] = indices[i];
1041: baij->nz = nz;
1042: for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1044: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1045: return 0;
1046: }
1048: /*@
1049: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1050: in the matrix.
1052: Input Parameters:
1053: + mat - the SeqSBAIJ matrix
1054: - indices - the column indices
1056: Level: advanced
1058: Notes:
1059: This can be called if you have precomputed the nonzero structure of the
1060: matrix and want to provide it to the matrix object to improve the performance
1061: of the MatSetValues() operation.
1063: You MUST have set the correct numbers of nonzeros per row in the call to
1064: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1066: MUST be called before any calls to MatSetValues()
1068: .seealso: MatCreateSeqSBAIJ
1069: @*/
1070: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1071: {
1074: PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
1075: return 0;
1076: }
1078: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1079: {
1080: PetscBool isbaij;
1082: PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");
1084: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1085: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1086: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1087: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1092: PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);
1093: PetscObjectStateIncrease((PetscObject)B);
1094: } else {
1095: MatGetRowUpperTriangular(A);
1096: MatCopy_Basic(A,B,str);
1097: MatRestoreRowUpperTriangular(A);
1098: }
1099: return 0;
1100: }
1102: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1103: {
1104: MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);
1105: return 0;
1106: }
1108: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1109: {
1110: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1112: *array = a->a;
1113: return 0;
1114: }
1116: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1117: {
1118: *array = NULL;
1119: return 0;
1120: }
1122: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1123: {
1124: PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1125: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1126: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1128: /* Set the number of nonzeros in the new matrix */
1129: MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
1130: return 0;
1131: }
1133: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1134: {
1135: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1136: PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1137: PetscBLASInt one = 1;
1139: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
1140: PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE;
1141: if (e) {
1142: PetscArraycmp(x->i,y->i,x->mbs+1,&e);
1143: if (e) {
1144: PetscArraycmp(x->j,y->j,x->i[x->mbs],&e);
1145: if (e) str = SAME_NONZERO_PATTERN;
1146: }
1147: }
1149: }
1150: if (str == SAME_NONZERO_PATTERN) {
1151: PetscScalar alpha = a;
1152: PetscBLASInt bnz;
1153: PetscBLASIntCast(x->nz*bs2,&bnz);
1154: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1155: PetscObjectStateIncrease((PetscObject)Y);
1156: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1157: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);
1158: MatAXPY_Basic(Y,a,X,str);
1159: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);
1160: } else {
1161: Mat B;
1162: PetscInt *nnz;
1164: MatGetRowUpperTriangular(X);
1165: MatGetRowUpperTriangular(Y);
1166: PetscMalloc1(Y->rmap->N,&nnz);
1167: MatCreate(PetscObjectComm((PetscObject)Y),&B);
1168: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
1169: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
1170: MatSetBlockSizesFromMats(B,Y,Y);
1171: MatSetType(B,((PetscObject)Y)->type_name);
1172: MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);
1173: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1175: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1177: MatHeaderMerge(Y,&B);
1178: PetscFree(nnz);
1179: MatRestoreRowUpperTriangular(X);
1180: MatRestoreRowUpperTriangular(Y);
1181: }
1182: return 0;
1183: }
1185: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1186: {
1187: *flg = PETSC_TRUE;
1188: return 0;
1189: }
1191: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1192: {
1193: *flg = PETSC_TRUE;
1194: return 0;
1195: }
1197: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1198: {
1199: *flg = PETSC_FALSE;
1200: return 0;
1201: }
1203: PetscErrorCode MatConjugate_SeqSBAIJ(Mat A)
1204: {
1205: #if defined(PETSC_USE_COMPLEX)
1206: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1207: PetscInt i,nz = a->bs2*a->i[a->mbs];
1208: MatScalar *aa = a->a;
1210: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1211: #else
1212: #endif
1213: return 0;
1214: }
1216: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1217: {
1218: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1219: PetscInt i,nz = a->bs2*a->i[a->mbs];
1220: MatScalar *aa = a->a;
1222: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1223: return 0;
1224: }
1226: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1227: {
1228: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1229: PetscInt i,nz = a->bs2*a->i[a->mbs];
1230: MatScalar *aa = a->a;
1232: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1233: return 0;
1234: }
1236: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1237: {
1238: Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1239: PetscInt i,j,k,count;
1240: PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1241: PetscScalar zero = 0.0;
1242: MatScalar *aa;
1243: const PetscScalar *xx;
1244: PetscScalar *bb;
1245: PetscBool *zeroed,vecs = PETSC_FALSE;
1247: /* fix right hand side if needed */
1248: if (x && b) {
1249: VecGetArrayRead(x,&xx);
1250: VecGetArray(b,&bb);
1251: vecs = PETSC_TRUE;
1252: }
1254: /* zero the columns */
1255: PetscCalloc1(A->rmap->n,&zeroed);
1256: for (i=0; i<is_n; i++) {
1258: zeroed[is_idx[i]] = PETSC_TRUE;
1259: }
1260: if (vecs) {
1261: for (i=0; i<A->rmap->N; i++) {
1262: row = i/bs;
1263: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1264: for (k=0; k<bs; k++) {
1265: col = bs*baij->j[j] + k;
1266: if (col <= i) continue;
1267: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1268: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1269: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1270: }
1271: }
1272: }
1273: for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1274: }
1276: for (i=0; i<A->rmap->N; i++) {
1277: if (!zeroed[i]) {
1278: row = i/bs;
1279: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1280: for (k=0; k<bs; k++) {
1281: col = bs*baij->j[j] + k;
1282: if (zeroed[col]) {
1283: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1284: aa[0] = 0.0;
1285: }
1286: }
1287: }
1288: }
1289: }
1290: PetscFree(zeroed);
1291: if (vecs) {
1292: VecRestoreArrayRead(x,&xx);
1293: VecRestoreArray(b,&bb);
1294: }
1296: /* zero the rows */
1297: for (i=0; i<is_n; i++) {
1298: row = is_idx[i];
1299: count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1300: aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1301: for (k=0; k<count; k++) {
1302: aa[0] = zero;
1303: aa += bs;
1304: }
1305: if (diag != 0.0) {
1306: (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
1307: }
1308: }
1309: MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);
1310: return 0;
1311: }
1313: PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1314: {
1315: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1317: if (!Y->preallocated || !aij->nz) {
1318: MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
1319: }
1320: MatShift_Basic(Y,a);
1321: return 0;
1322: }
1324: /* -------------------------------------------------------------------*/
1325: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1326: MatGetRow_SeqSBAIJ,
1327: MatRestoreRow_SeqSBAIJ,
1328: MatMult_SeqSBAIJ_N,
1329: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1330: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1331: MatMultAdd_SeqSBAIJ_N,
1332: NULL,
1333: NULL,
1334: NULL,
1335: /* 10*/ NULL,
1336: NULL,
1337: MatCholeskyFactor_SeqSBAIJ,
1338: MatSOR_SeqSBAIJ,
1339: MatTranspose_SeqSBAIJ,
1340: /* 15*/ MatGetInfo_SeqSBAIJ,
1341: MatEqual_SeqSBAIJ,
1342: MatGetDiagonal_SeqSBAIJ,
1343: MatDiagonalScale_SeqSBAIJ,
1344: MatNorm_SeqSBAIJ,
1345: /* 20*/ NULL,
1346: MatAssemblyEnd_SeqSBAIJ,
1347: MatSetOption_SeqSBAIJ,
1348: MatZeroEntries_SeqSBAIJ,
1349: /* 24*/ NULL,
1350: NULL,
1351: NULL,
1352: NULL,
1353: NULL,
1354: /* 29*/ MatSetUp_SeqSBAIJ,
1355: NULL,
1356: NULL,
1357: NULL,
1358: NULL,
1359: /* 34*/ MatDuplicate_SeqSBAIJ,
1360: NULL,
1361: NULL,
1362: NULL,
1363: MatICCFactor_SeqSBAIJ,
1364: /* 39*/ MatAXPY_SeqSBAIJ,
1365: MatCreateSubMatrices_SeqSBAIJ,
1366: MatIncreaseOverlap_SeqSBAIJ,
1367: MatGetValues_SeqSBAIJ,
1368: MatCopy_SeqSBAIJ,
1369: /* 44*/ NULL,
1370: MatScale_SeqSBAIJ,
1371: MatShift_SeqSBAIJ,
1372: NULL,
1373: MatZeroRowsColumns_SeqSBAIJ,
1374: /* 49*/ NULL,
1375: MatGetRowIJ_SeqSBAIJ,
1376: MatRestoreRowIJ_SeqSBAIJ,
1377: NULL,
1378: NULL,
1379: /* 54*/ NULL,
1380: NULL,
1381: NULL,
1382: MatPermute_SeqSBAIJ,
1383: MatSetValuesBlocked_SeqSBAIJ,
1384: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1385: NULL,
1386: NULL,
1387: NULL,
1388: NULL,
1389: /* 64*/ NULL,
1390: NULL,
1391: NULL,
1392: NULL,
1393: NULL,
1394: /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1395: NULL,
1396: MatConvert_MPISBAIJ_Basic,
1397: NULL,
1398: NULL,
1399: /* 74*/ NULL,
1400: NULL,
1401: NULL,
1402: NULL,
1403: NULL,
1404: /* 79*/ NULL,
1405: NULL,
1406: NULL,
1407: MatGetInertia_SeqSBAIJ,
1408: MatLoad_SeqSBAIJ,
1409: /* 84*/ MatIsSymmetric_SeqSBAIJ,
1410: MatIsHermitian_SeqSBAIJ,
1411: MatIsStructurallySymmetric_SeqSBAIJ,
1412: NULL,
1413: NULL,
1414: /* 89*/ NULL,
1415: NULL,
1416: NULL,
1417: NULL,
1418: NULL,
1419: /* 94*/ NULL,
1420: NULL,
1421: NULL,
1422: NULL,
1423: NULL,
1424: /* 99*/ NULL,
1425: NULL,
1426: NULL,
1427: MatConjugate_SeqSBAIJ,
1428: NULL,
1429: /*104*/ NULL,
1430: MatRealPart_SeqSBAIJ,
1431: MatImaginaryPart_SeqSBAIJ,
1432: MatGetRowUpperTriangular_SeqSBAIJ,
1433: MatRestoreRowUpperTriangular_SeqSBAIJ,
1434: /*109*/ NULL,
1435: NULL,
1436: NULL,
1437: NULL,
1438: MatMissingDiagonal_SeqSBAIJ,
1439: /*114*/ NULL,
1440: NULL,
1441: NULL,
1442: NULL,
1443: NULL,
1444: /*119*/ NULL,
1445: NULL,
1446: NULL,
1447: NULL,
1448: NULL,
1449: /*124*/ NULL,
1450: NULL,
1451: NULL,
1452: NULL,
1453: NULL,
1454: /*129*/ NULL,
1455: NULL,
1456: NULL,
1457: NULL,
1458: NULL,
1459: /*134*/ NULL,
1460: NULL,
1461: NULL,
1462: NULL,
1463: NULL,
1464: /*139*/ MatSetBlockSizes_Default,
1465: NULL,
1466: NULL,
1467: NULL,
1468: NULL,
1469: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ,
1470: NULL,
1471: NULL,
1472: NULL
1473: };
1475: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1476: {
1477: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1478: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1482: /* allocate space for values if not already there */
1483: if (!aij->saved_values) {
1484: PetscMalloc1(nz+1,&aij->saved_values);
1485: }
1487: /* copy values over */
1488: PetscArraycpy(aij->saved_values,aij->a,nz);
1489: return 0;
1490: }
1492: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1493: {
1494: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1495: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1500: /* copy values over */
1501: PetscArraycpy(aij->a,aij->saved_values,nz);
1502: return 0;
1503: }
1505: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1506: {
1507: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1508: PetscInt i,mbs,nbs,bs2;
1509: PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1511: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1513: MatSetBlockSize(B,PetscAbs(bs));
1514: PetscLayoutSetUp(B->rmap);
1515: PetscLayoutSetUp(B->cmap);
1517: PetscLayoutGetBlockSize(B->rmap,&bs);
1519: B->preallocated = PETSC_TRUE;
1521: mbs = B->rmap->N/bs;
1522: nbs = B->cmap->n/bs;
1523: bs2 = bs*bs;
1527: if (nz == MAT_SKIP_ALLOCATION) {
1528: skipallocation = PETSC_TRUE;
1529: nz = 0;
1530: }
1532: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1534: if (nnz) {
1535: for (i=0; i<mbs; i++) {
1538: }
1539: }
1541: B->ops->mult = MatMult_SeqSBAIJ_N;
1542: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1543: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1544: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1546: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1547: if (!flg) {
1548: switch (bs) {
1549: case 1:
1550: B->ops->mult = MatMult_SeqSBAIJ_1;
1551: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1552: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1553: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1554: break;
1555: case 2:
1556: B->ops->mult = MatMult_SeqSBAIJ_2;
1557: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1558: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1559: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1560: break;
1561: case 3:
1562: B->ops->mult = MatMult_SeqSBAIJ_3;
1563: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1564: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1565: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1566: break;
1567: case 4:
1568: B->ops->mult = MatMult_SeqSBAIJ_4;
1569: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1570: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1571: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1572: break;
1573: case 5:
1574: B->ops->mult = MatMult_SeqSBAIJ_5;
1575: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1576: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1577: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1578: break;
1579: case 6:
1580: B->ops->mult = MatMult_SeqSBAIJ_6;
1581: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1582: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1583: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1584: break;
1585: case 7:
1586: B->ops->mult = MatMult_SeqSBAIJ_7;
1587: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1588: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1589: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1590: break;
1591: }
1592: }
1594: b->mbs = mbs;
1595: b->nbs = nbs;
1596: if (!skipallocation) {
1597: if (!b->imax) {
1598: PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
1600: b->free_imax_ilen = PETSC_TRUE;
1602: PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));
1603: }
1604: if (!nnz) {
1605: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1606: else if (nz <= 0) nz = 1;
1607: nz = PetscMin(nbs,nz);
1608: for (i=0; i<mbs; i++) b->imax[i] = nz;
1609: PetscIntMultError(nz,mbs,&nz);
1610: } else {
1611: PetscInt64 nz64 = 0;
1612: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1613: PetscIntCast(nz64,&nz);
1614: }
1615: /* b->ilen will count nonzeros in each block row so far. */
1616: for (i=0; i<mbs; i++) b->ilen[i] = 0;
1617: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1619: /* allocate the matrix space */
1620: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
1621: PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
1622: PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
1623: PetscArrayzero(b->a,nz*bs2);
1624: PetscArrayzero(b->j,nz);
1626: b->singlemalloc = PETSC_TRUE;
1628: /* pointer to beginning of each row */
1629: b->i[0] = 0;
1630: for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1632: b->free_a = PETSC_TRUE;
1633: b->free_ij = PETSC_TRUE;
1634: } else {
1635: b->free_a = PETSC_FALSE;
1636: b->free_ij = PETSC_FALSE;
1637: }
1639: b->bs2 = bs2;
1640: b->nz = 0;
1641: b->maxnz = nz;
1642: b->inew = NULL;
1643: b->jnew = NULL;
1644: b->anew = NULL;
1645: b->a2anew = NULL;
1646: b->permute = PETSC_FALSE;
1648: B->was_assembled = PETSC_FALSE;
1649: B->assembled = PETSC_FALSE;
1650: if (realalloc) MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1651: return 0;
1652: }
1654: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1655: {
1656: PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1657: PetscScalar *values=NULL;
1658: PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1661: PetscLayoutSetBlockSize(B->rmap,bs);
1662: PetscLayoutSetBlockSize(B->cmap,bs);
1663: PetscLayoutSetUp(B->rmap);
1664: PetscLayoutSetUp(B->cmap);
1665: PetscLayoutGetBlockSize(B->rmap,&bs);
1666: m = B->rmap->n/bs;
1669: PetscMalloc1(m+1,&nnz);
1670: for (i=0; i<m; i++) {
1671: nz = ii[i+1] - ii[i];
1673: anz = 0;
1674: for (j=0; j<nz; j++) {
1675: /* count only values on the diagonal or above */
1676: if (jj[ii[i] + j] >= i) {
1677: anz = nz - j;
1678: break;
1679: }
1680: }
1681: nz_max = PetscMax(nz_max,anz);
1682: nnz[i] = anz;
1683: }
1684: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1685: PetscFree(nnz);
1687: values = (PetscScalar*)V;
1688: if (!values) {
1689: PetscCalloc1(bs*bs*nz_max,&values);
1690: }
1691: for (i=0; i<m; i++) {
1692: PetscInt ncols = ii[i+1] - ii[i];
1693: const PetscInt *icols = jj + ii[i];
1694: if (!roworiented || bs == 1) {
1695: const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1696: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
1697: } else {
1698: for (j=0; j<ncols; j++) {
1699: const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1700: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
1701: }
1702: }
1703: }
1704: if (!V) PetscFree(values);
1705: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1706: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1707: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1708: return 0;
1709: }
1711: /*
1712: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1713: */
1714: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1715: {
1716: PetscBool flg = PETSC_FALSE;
1717: PetscInt bs = B->rmap->bs;
1719: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1720: if (flg) bs = 8;
1722: if (!natural) {
1723: switch (bs) {
1724: case 1:
1725: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1726: break;
1727: case 2:
1728: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1729: break;
1730: case 3:
1731: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1732: break;
1733: case 4:
1734: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1735: break;
1736: case 5:
1737: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1738: break;
1739: case 6:
1740: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1741: break;
1742: case 7:
1743: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1744: break;
1745: default:
1746: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1747: break;
1748: }
1749: } else {
1750: switch (bs) {
1751: case 1:
1752: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1753: break;
1754: case 2:
1755: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1756: break;
1757: case 3:
1758: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1759: break;
1760: case 4:
1761: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1762: break;
1763: case 5:
1764: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1765: break;
1766: case 6:
1767: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1768: break;
1769: case 7:
1770: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1771: break;
1772: default:
1773: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1774: break;
1775: }
1776: }
1777: return 0;
1778: }
1780: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*);
1781: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*);
1782: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type)
1783: {
1784: *type = MATSOLVERPETSC;
1785: return 0;
1786: }
1788: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1789: {
1790: PetscInt n = A->rmap->n;
1792: #if defined(PETSC_USE_COMPLEX)
1794: #endif
1796: MatCreate(PetscObjectComm((PetscObject)A),B);
1797: MatSetSizes(*B,n,n,n,n);
1798: if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1799: MatSetType(*B,MATSEQSBAIJ);
1800: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
1802: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1803: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1804: PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
1805: PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
1806: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1808: (*B)->factortype = ftype;
1809: (*B)->canuseordering = PETSC_TRUE;
1810: PetscFree((*B)->solvertype);
1811: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
1812: PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);
1813: return 0;
1814: }
1816: /*@C
1817: MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1819: Not Collective
1821: Input Parameter:
1822: . mat - a MATSEQSBAIJ matrix
1824: Output Parameter:
1825: . array - pointer to the data
1827: Level: intermediate
1829: .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1830: @*/
1831: PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1832: {
1833: PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
1834: return 0;
1835: }
1837: /*@C
1838: MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1840: Not Collective
1842: Input Parameters:
1843: + mat - a MATSEQSBAIJ matrix
1844: - array - pointer to the data
1846: Level: intermediate
1848: .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1849: @*/
1850: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1851: {
1852: PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
1853: return 0;
1854: }
1856: /*MC
1857: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1858: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1860: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1861: can call MatSetOption(Mat, MAT_HERMITIAN).
1863: Options Database Keys:
1864: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1866: Notes:
1867: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1868: stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1869: 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.
1871: The number of rows in the matrix must be less than or equal to the number of columns
1873: Level: beginner
1875: .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1876: M*/
1877: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1878: {
1879: Mat_SeqSBAIJ *b;
1881: PetscMPIInt size;
1882: PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1884: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1887: PetscNewLog(B,&b);
1888: B->data = (void*)b;
1889: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1891: B->ops->destroy = MatDestroy_SeqSBAIJ;
1892: B->ops->view = MatView_SeqSBAIJ;
1893: b->row = NULL;
1894: b->icol = NULL;
1895: b->reallocs = 0;
1896: b->saved_values = NULL;
1897: b->inode.limit = 5;
1898: b->inode.max_limit = 5;
1900: b->roworiented = PETSC_TRUE;
1901: b->nonew = 0;
1902: b->diag = NULL;
1903: b->solve_work = NULL;
1904: b->mult_work = NULL;
1905: B->spptr = NULL;
1906: B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1907: b->keepnonzeropattern = PETSC_FALSE;
1909: b->inew = NULL;
1910: b->jnew = NULL;
1911: b->anew = NULL;
1912: b->a2anew = NULL;
1913: b->permute = PETSC_FALSE;
1915: b->ignore_ltriangular = PETSC_TRUE;
1917: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);
1919: b->getrow_utriangular = PETSC_FALSE;
1921: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);
1923: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);
1924: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);
1925: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);
1926: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);
1927: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1928: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);
1929: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);
1930: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);
1931: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
1932: #if defined(PETSC_HAVE_ELEMENTAL)
1933: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);
1934: #endif
1935: #if defined(PETSC_HAVE_SCALAPACK)
1936: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);
1937: #endif
1939: B->symmetric = PETSC_TRUE;
1940: B->structurally_symmetric = PETSC_TRUE;
1941: B->symmetric_set = PETSC_TRUE;
1942: B->structurally_symmetric_set = PETSC_TRUE;
1943: B->symmetric_eternal = PETSC_TRUE;
1944: #if defined(PETSC_USE_COMPLEX)
1945: B->hermitian = PETSC_FALSE;
1946: B->hermitian_set = PETSC_FALSE;
1947: #else
1948: B->hermitian = PETSC_TRUE;
1949: B->hermitian_set = PETSC_TRUE;
1950: #endif
1952: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
1954: PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");
1955: PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);
1956: if (no_unroll) {
1957: PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");
1958: }
1959: PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);
1960: if (no_inode) PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");
1961: PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);
1962: PetscOptionsEnd();
1963: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
1964: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
1965: return 0;
1966: }
1968: /*@C
1969: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
1970: compressed row) format. For good matrix assembly performance the
1971: user should preallocate the matrix storage by setting the parameter nz
1972: (or the array nnz). By setting these parameters accurately, performance
1973: during matrix assembly can be increased by more than a factor of 50.
1975: Collective on Mat
1977: Input Parameters:
1978: + B - the symmetric matrix
1979: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
1980: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
1981: . nz - number of block nonzeros per block row (same for all rows)
1982: - nnz - array containing the number of block nonzeros in the upper triangular plus
1983: diagonal portion of each block (possibly different for each block row) or NULL
1985: Options Database Keys:
1986: + -mat_no_unroll - uses code that does not unroll the loops in the
1987: block calculations (much slower)
1988: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
1990: Level: intermediate
1992: Notes:
1993: Specify the preallocated storage with either nz or nnz (not both).
1994: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
1995: allocation. See Users-Manual: ch_mat for details.
1997: You can call MatGetInfo() to get information on how effective the preallocation was;
1998: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
1999: You can also run with the option -info and look for messages with the string
2000: malloc in them to see if additional memory allocation was needed.
2002: If the nnz parameter is given then the nz parameter is ignored
2004: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2005: @*/
2006: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2007: {
2011: PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
2012: return 0;
2013: }
2015: /*@C
2016: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2018: Input Parameters:
2019: + B - the matrix
2020: . bs - size of block, the blocks are ALWAYS square.
2021: . i - the indices into j for the start of each local row (starts with zero)
2022: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2023: - v - optional values in the matrix
2025: Level: advanced
2027: Notes:
2028: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2029: may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2030: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2031: MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2032: block column and the second index is over columns within a block.
2034: Any entries below the diagonal are ignored
2036: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2037: and usually the numerical values as well
2039: .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2040: @*/
2041: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2042: {
2046: PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
2047: return 0;
2048: }
2050: /*@C
2051: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2052: compressed row) format. For good matrix assembly performance the
2053: user should preallocate the matrix storage by setting the parameter nz
2054: (or the array nnz). By setting these parameters accurately, performance
2055: during matrix assembly can be increased by more than a factor of 50.
2057: Collective
2059: Input Parameters:
2060: + comm - MPI communicator, set to PETSC_COMM_SELF
2061: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2062: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2063: . m - number of rows, or number of columns
2064: . nz - number of block nonzeros per block row (same for all rows)
2065: - nnz - array containing the number of block nonzeros in the upper triangular plus
2066: diagonal portion of each block (possibly different for each block row) or NULL
2068: Output Parameter:
2069: . A - the symmetric matrix
2071: Options Database Keys:
2072: + -mat_no_unroll - uses code that does not unroll the loops in the
2073: block calculations (much slower)
2074: - -mat_block_size - size of the blocks to use
2076: Level: intermediate
2078: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2079: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2080: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2082: Notes:
2083: The number of rows and columns must be divisible by blocksize.
2084: This matrix type does not support complex Hermitian operation.
2086: Specify the preallocated storage with either nz or nnz (not both).
2087: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2088: allocation. See Users-Manual: ch_mat for details.
2090: If the nnz parameter is given then the nz parameter is ignored
2092: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2093: @*/
2094: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2095: {
2096: MatCreate(comm,A);
2097: MatSetSizes(*A,m,n,m,n);
2098: MatSetType(*A,MATSEQSBAIJ);
2099: MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
2100: return 0;
2101: }
2103: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2104: {
2105: Mat C;
2106: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2107: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2111: *B = NULL;
2112: MatCreate(PetscObjectComm((PetscObject)A),&C);
2113: MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
2114: MatSetBlockSizesFromMats(C,A,A);
2115: MatSetType(C,MATSEQSBAIJ);
2116: c = (Mat_SeqSBAIJ*)C->data;
2118: C->preallocated = PETSC_TRUE;
2119: C->factortype = A->factortype;
2120: c->row = NULL;
2121: c->icol = NULL;
2122: c->saved_values = NULL;
2123: c->keepnonzeropattern = a->keepnonzeropattern;
2124: C->assembled = PETSC_TRUE;
2126: PetscLayoutReference(A->rmap,&C->rmap);
2127: PetscLayoutReference(A->cmap,&C->cmap);
2128: c->bs2 = a->bs2;
2129: c->mbs = a->mbs;
2130: c->nbs = a->nbs;
2132: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2133: c->imax = a->imax;
2134: c->ilen = a->ilen;
2135: c->free_imax_ilen = PETSC_FALSE;
2136: } else {
2137: PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);
2138: PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));
2139: for (i=0; i<mbs; i++) {
2140: c->imax[i] = a->imax[i];
2141: c->ilen[i] = a->ilen[i];
2142: }
2143: c->free_imax_ilen = PETSC_TRUE;
2144: }
2146: /* allocate the matrix space */
2147: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2148: PetscMalloc1(bs2*nz,&c->a);
2149: PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));
2150: c->i = a->i;
2151: c->j = a->j;
2152: c->singlemalloc = PETSC_FALSE;
2153: c->free_a = PETSC_TRUE;
2154: c->free_ij = PETSC_FALSE;
2155: c->parent = A;
2156: PetscObjectReference((PetscObject)A);
2157: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2158: MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2159: } else {
2160: PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
2161: PetscArraycpy(c->i,a->i,mbs+1);
2162: PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
2163: c->singlemalloc = PETSC_TRUE;
2164: c->free_a = PETSC_TRUE;
2165: c->free_ij = PETSC_TRUE;
2166: }
2167: if (mbs > 0) {
2168: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2169: PetscArraycpy(c->j,a->j,nz);
2170: }
2171: if (cpvalues == MAT_COPY_VALUES) {
2172: PetscArraycpy(c->a,a->a,bs2*nz);
2173: } else {
2174: PetscArrayzero(c->a,bs2*nz);
2175: }
2176: if (a->jshort) {
2177: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2178: /* if the parent matrix is reassembled, this child matrix will never notice */
2179: PetscMalloc1(nz,&c->jshort);
2180: PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));
2181: PetscArraycpy(c->jshort,a->jshort,nz);
2183: c->free_jshort = PETSC_TRUE;
2184: }
2185: }
2187: c->roworiented = a->roworiented;
2188: c->nonew = a->nonew;
2190: if (a->diag) {
2191: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2192: c->diag = a->diag;
2193: c->free_diag = PETSC_FALSE;
2194: } else {
2195: PetscMalloc1(mbs,&c->diag);
2196: PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));
2197: for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2198: c->free_diag = PETSC_TRUE;
2199: }
2200: }
2201: c->nz = a->nz;
2202: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2203: c->solve_work = NULL;
2204: c->mult_work = NULL;
2206: *B = C;
2207: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
2208: return 0;
2209: }
2211: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2212: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2214: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)
2215: {
2216: PetscBool isbinary;
2218: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2220: MatLoad_SeqSBAIJ_Binary(mat,viewer);
2221: return 0;
2222: }
2224: /*@
2225: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2226: (upper triangular entries in CSR format) provided by the user.
2228: Collective
2230: Input Parameters:
2231: + comm - must be an MPI communicator of size 1
2232: . bs - size of block
2233: . m - number of rows
2234: . n - number of columns
2235: . 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
2236: . j - column indices
2237: - a - matrix values
2239: Output Parameter:
2240: . mat - the matrix
2242: Level: advanced
2244: Notes:
2245: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2246: once the matrix is destroyed
2248: You cannot set new nonzero locations into this matrix, that will generate an error.
2250: The i and j indices are 0 based
2252: 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
2253: it is the regular CSR format excluding the lower triangular elements.
2255: .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2257: @*/
2258: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2259: {
2260: PetscInt ii;
2261: Mat_SeqSBAIJ *sbaij;
2266: MatCreate(comm,mat);
2267: MatSetSizes(*mat,m,n,m,n);
2268: MatSetType(*mat,MATSEQSBAIJ);
2269: MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);
2270: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2271: PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);
2272: PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));
2274: sbaij->i = i;
2275: sbaij->j = j;
2276: sbaij->a = a;
2278: sbaij->singlemalloc = PETSC_FALSE;
2279: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2280: sbaij->free_a = PETSC_FALSE;
2281: sbaij->free_ij = PETSC_FALSE;
2282: sbaij->free_imax_ilen = PETSC_TRUE;
2284: for (ii=0; ii<m; ii++) {
2285: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2287: }
2288: if (PetscDefined(USE_DEBUG)) {
2289: for (ii=0; ii<sbaij->i[m]; ii++) {
2292: }
2293: }
2295: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2296: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2297: return 0;
2298: }
2300: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2301: {
2302: PetscMPIInt size;
2304: MPI_Comm_size(comm,&size);
2305: if (size == 1 && scall == MAT_REUSE_MATRIX) {
2306: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2307: } else {
2308: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);
2309: }
2310: return 0;
2311: }