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;
29: PetscInt *diag,*ii = a->i,i;
32: MatMarkDiagonal_SeqSBAIJ(A);
33: *missing = PETSC_FALSE;
34: if (A->rmap->n > 0 && !ii) {
35: *missing = PETSC_TRUE;
36: if (dd) *dd = 0;
37: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
38: } else {
39: diag = a->diag;
40: for (i=0; i<a->mbs; i++) {
41: if (diag[i] >= ii[i+1]) {
42: *missing = PETSC_TRUE;
43: if (dd) *dd = i;
44: break;
45: }
46: }
47: }
48: return(0);
49: }
51: PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
52: {
53: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
55: PetscInt i,j;
58: if (!a->diag) {
59: PetscMalloc1(a->mbs,&a->diag);
60: PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt));
61: a->free_diag = PETSC_TRUE;
62: }
63: for (i=0; i<a->mbs; i++) {
64: a->diag[i] = a->i[i+1];
65: for (j=a->i[i]; j<a->i[i+1]; j++) {
66: if (a->j[j] == i) {
67: a->diag[i] = j;
68: break;
69: }
70: }
71: }
72: return(0);
73: }
75: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done)
76: {
77: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
79: PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt;
80: PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja;
83: *nn = n;
84: if (!ia) return(0);
85: if (symmetric) {
86: MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja);
87: nz = tia[n];
88: } else {
89: tia = a->i; tja = a->j;
90: }
92: if (!blockcompressed && bs > 1) {
93: (*nn) *= bs;
94: /* malloc & create the natural set of indices */
95: PetscMalloc1((n+1)*bs,ia);
96: if (n) {
97: (*ia)[0] = oshift;
98: for (j=1; j<bs; j++) {
99: (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1];
100: }
101: }
103: for (i=1; i<n; i++) {
104: (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1];
105: for (j=1; j<bs; j++) {
106: (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1];
107: }
108: }
109: if (n) {
110: (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1];
111: }
113: if (inja) {
114: PetscMalloc1(nz*bs*bs,ja);
115: cnt = 0;
116: for (i=0; i<n; i++) {
117: for (j=0; j<bs; j++) {
118: for (k=tia[i]; k<tia[i+1]; k++) {
119: for (l=0; l<bs; l++) {
120: (*ja)[cnt++] = bs*tja[k] + l;
121: }
122: }
123: }
124: }
125: }
127: if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
128: PetscFree(tia);
129: PetscFree(tja);
130: }
131: } else if (oshift == 1) {
132: if (symmetric) {
133: nz = tia[A->rmap->n/bs];
134: /* add 1 to i and j indices */
135: for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1;
136: *ia = tia;
137: if (ja) {
138: for (i=0; i<nz; i++) tja[i] = tja[i] + 1;
139: *ja = tja;
140: }
141: } else {
142: nz = a->i[A->rmap->n/bs];
143: /* malloc space and add 1 to i and j indices */
144: PetscMalloc1(A->rmap->n/bs+1,ia);
145: for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1;
146: if (ja) {
147: PetscMalloc1(nz,ja);
148: for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1;
149: }
150: }
151: } else {
152: *ia = tia;
153: if (ja) *ja = tja;
154: }
155: return(0);
156: }
158: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
159: {
163: if (!ia) return(0);
164: if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
165: PetscFree(*ia);
166: if (ja) {PetscFree(*ja);}
167: }
168: return(0);
169: }
171: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
172: {
173: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
177: #if defined(PETSC_USE_LOG)
178: PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz);
179: #endif
180: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
181: if (a->free_diag) {PetscFree(a->diag);}
182: ISDestroy(&a->row);
183: ISDestroy(&a->col);
184: ISDestroy(&a->icol);
185: PetscFree(a->idiag);
186: PetscFree(a->inode.size);
187: if (a->free_imax_ilen) {PetscFree2(a->imax,a->ilen);}
188: PetscFree(a->solve_work);
189: PetscFree(a->sor_work);
190: PetscFree(a->solves_work);
191: PetscFree(a->mult_work);
192: PetscFree(a->saved_values);
193: if (a->free_jshort) {PetscFree(a->jshort);}
194: PetscFree(a->inew);
195: MatDestroy(&a->parent);
196: PetscFree(A->data);
198: PetscObjectChangeTypeName((PetscObject)A,NULL);
199: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
200: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
201: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL);
202: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL);
203: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL);
204: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL);
205: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL);
206: #if defined(PETSC_HAVE_ELEMENTAL)
207: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL);
208: #endif
209: #if defined(PETSC_HAVE_SCALAPACK)
210: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_scalapack_C",NULL);
211: #endif
212: return(0);
213: }
215: PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg)
216: {
217: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
218: #if defined(PETSC_USE_COMPLEX)
219: PetscInt bs;
220: #endif
224: #if defined(PETSC_USE_COMPLEX)
225: MatGetBlockSize(A,&bs);
226: #endif
227: switch (op) {
228: case MAT_ROW_ORIENTED:
229: a->roworiented = flg;
230: break;
231: case MAT_KEEP_NONZERO_PATTERN:
232: a->keepnonzeropattern = flg;
233: break;
234: case MAT_NEW_NONZERO_LOCATIONS:
235: a->nonew = (flg ? 0 : 1);
236: break;
237: case MAT_NEW_NONZERO_LOCATION_ERR:
238: a->nonew = (flg ? -1 : 0);
239: break;
240: case MAT_NEW_NONZERO_ALLOCATION_ERR:
241: a->nonew = (flg ? -2 : 0);
242: break;
243: case MAT_UNUSED_NONZERO_LOCATION_ERR:
244: a->nounused = (flg ? -1 : 0);
245: break;
246: case MAT_FORCE_DIAGONAL_ENTRIES:
247: case MAT_IGNORE_OFF_PROC_ENTRIES:
248: case MAT_USE_HASH_TABLE:
249: case MAT_SORTED_FULL:
250: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
251: break;
252: case MAT_HERMITIAN:
253: #if defined(PETSC_USE_COMPLEX)
254: if (flg) { /* disable transpose ops */
255: if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1");
256: A->ops->multtranspose = NULL;
257: A->ops->multtransposeadd = NULL;
258: A->symmetric = PETSC_FALSE;
259: }
260: #endif
261: break;
262: case MAT_SYMMETRIC:
263: case MAT_SPD:
264: #if defined(PETSC_USE_COMPLEX)
265: if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */
266: A->ops->multtranspose = A->ops->mult;
267: A->ops->multtransposeadd = A->ops->multadd;
268: }
269: #endif
270: break;
271: /* These options are handled directly by MatSetOption() */
272: case MAT_STRUCTURALLY_SYMMETRIC:
273: case MAT_SYMMETRY_ETERNAL:
274: case MAT_STRUCTURE_ONLY:
275: /* These options are handled directly by MatSetOption() */
276: break;
277: case MAT_IGNORE_LOWER_TRIANGULAR:
278: a->ignore_ltriangular = flg;
279: break;
280: case MAT_ERROR_LOWER_TRIANGULAR:
281: a->ignore_ltriangular = flg;
282: break;
283: case MAT_GETROW_UPPERTRIANGULAR:
284: a->getrow_utriangular = flg;
285: break;
286: case MAT_SUBMAT_SINGLEIS:
287: break;
288: default:
289: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
290: }
291: return(0);
292: }
294: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
295: {
296: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
300: 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()");
302: /* Get the upper triangular part of the row */
303: MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);
304: return(0);
305: }
307: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
308: {
312: if (nz) *nz = 0;
313: if (idx) {PetscFree(*idx);}
314: if (v) {PetscFree(*v);}
315: return(0);
316: }
318: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
319: {
320: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
323: a->getrow_utriangular = PETSC_TRUE;
324: return(0);
325: }
327: PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
328: {
329: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
332: a->getrow_utriangular = PETSC_FALSE;
333: return(0);
334: }
336: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B)
337: {
341: if (reuse == MAT_INITIAL_MATRIX) {
342: MatDuplicate(A,MAT_COPY_VALUES,B);
343: } else if (reuse == MAT_REUSE_MATRIX) {
344: MatCopy(A,*B,SAME_NONZERO_PATTERN);
345: }
346: return(0);
347: }
349: PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
350: {
351: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
352: PetscErrorCode ierr;
353: PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
354: PetscViewerFormat format;
355: PetscInt *diag;
358: PetscViewerGetFormat(viewer,&format);
359: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
360: PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);
361: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
362: Mat aij;
363: const char *matname;
365: if (A->factortype && bs>1) {
366: PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");
367: return(0);
368: }
369: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);
370: PetscObjectGetName((PetscObject)A,&matname);
371: PetscObjectSetName((PetscObject)aij,matname);
372: MatView(aij,viewer);
373: MatDestroy(&aij);
374: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
375: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
376: for (i=0; i<a->mbs; i++) {
377: for (j=0; j<bs; j++) {
378: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
379: for (k=a->i[i]; k<a->i[i+1]; k++) {
380: for (l=0; l<bs; l++) {
381: #if defined(PETSC_USE_COMPLEX)
382: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
383: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
384: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
385: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
386: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
387: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
388: } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
389: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
390: }
391: #else
392: if (a->a[bs2*k + l*bs + j] != 0.0) {
393: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
394: }
395: #endif
396: }
397: }
398: PetscViewerASCIIPrintf(viewer,"\n");
399: }
400: }
401: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
402: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
403: return(0);
404: } else {
405: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
406: if (A->factortype) { /* for factored matrix */
407: if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet");
409: diag=a->diag;
410: for (i=0; i<a->mbs; i++) { /* for row block i */
411: PetscViewerASCIIPrintf(viewer,"row %D:",i);
412: /* diagonal entry */
413: #if defined(PETSC_USE_COMPLEX)
414: if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
415: 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]]));
416: } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
417: 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]]));
418: } else {
419: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));
420: }
421: #else
422: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));
423: #endif
424: /* off-diagonal entries */
425: for (k=a->i[i]; k<a->i[i+1]-1; k++) {
426: #if defined(PETSC_USE_COMPLEX)
427: if (PetscImaginaryPart(a->a[k]) > 0.0) {
428: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));
429: } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
430: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));
431: } else {
432: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));
433: }
434: #else
435: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);
436: #endif
437: }
438: PetscViewerASCIIPrintf(viewer,"\n");
439: }
441: } else { /* for non-factored matrix */
442: for (i=0; i<a->mbs; i++) { /* for row block i */
443: for (j=0; j<bs; j++) { /* for row bs*i + j */
444: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
445: for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */
446: for (l=0; l<bs; l++) { /* for column */
447: #if defined(PETSC_USE_COMPLEX)
448: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
449: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
450: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
451: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
452: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
453: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
454: } else {
455: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
456: }
457: #else
458: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
459: #endif
460: }
461: }
462: PetscViewerASCIIPrintf(viewer,"\n");
463: }
464: }
465: }
466: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
467: }
468: PetscViewerFlush(viewer);
469: return(0);
470: }
472: #include <petscdraw.h>
473: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
474: {
475: Mat A = (Mat) Aa;
476: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
478: PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
479: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
480: MatScalar *aa;
481: PetscViewer viewer;
484: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
485: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
487: /* loop over matrix elements drawing boxes */
489: PetscDrawCollectiveBegin(draw);
490: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
491: /* Blue for negative, Cyan for zero and Red for positive */
492: color = PETSC_DRAW_BLUE;
493: for (i=0,row=0; i<mbs; i++,row+=bs) {
494: for (j=a->i[i]; j<a->i[i+1]; j++) {
495: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
496: x_l = a->j[j]*bs; x_r = x_l + 1.0;
497: aa = a->a + j*bs2;
498: for (k=0; k<bs; k++) {
499: for (l=0; l<bs; l++) {
500: if (PetscRealPart(*aa++) >= 0.) continue;
501: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
502: }
503: }
504: }
505: }
506: color = PETSC_DRAW_CYAN;
507: for (i=0,row=0; i<mbs; i++,row+=bs) {
508: for (j=a->i[i]; j<a->i[i+1]; j++) {
509: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
510: x_l = a->j[j]*bs; x_r = x_l + 1.0;
511: aa = a->a + j*bs2;
512: for (k=0; k<bs; k++) {
513: for (l=0; l<bs; l++) {
514: if (PetscRealPart(*aa++) != 0.) continue;
515: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
516: }
517: }
518: }
519: }
520: color = PETSC_DRAW_RED;
521: for (i=0,row=0; i<mbs; i++,row+=bs) {
522: for (j=a->i[i]; j<a->i[i+1]; j++) {
523: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
524: x_l = a->j[j]*bs; x_r = x_l + 1.0;
525: aa = a->a + j*bs2;
526: for (k=0; k<bs; k++) {
527: for (l=0; l<bs; l++) {
528: if (PetscRealPart(*aa++) <= 0.) continue;
529: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
530: }
531: }
532: }
533: }
534: PetscDrawCollectiveEnd(draw);
535: return(0);
536: }
538: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
539: {
541: PetscReal xl,yl,xr,yr,w,h;
542: PetscDraw draw;
543: PetscBool isnull;
546: PetscViewerDrawGetDraw(viewer,0,&draw);
547: PetscDrawIsNull(draw,&isnull);
548: if (isnull) return(0);
550: xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
551: xr += w; yr += h; xl = -w; yl = -h;
552: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
553: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
554: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
555: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
556: PetscDrawSave(draw);
557: return(0);
558: }
560: /* Used for both MPIBAIJ and MPISBAIJ matrices */
561: #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary
563: PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
564: {
566: PetscBool iascii,isbinary,isdraw;
569: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
570: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
571: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
572: if (iascii) {
573: MatView_SeqSBAIJ_ASCII(A,viewer);
574: } else if (isbinary) {
575: MatView_SeqSBAIJ_Binary(A,viewer);
576: } else if (isdraw) {
577: MatView_SeqSBAIJ_Draw(A,viewer);
578: } else {
579: Mat B;
580: const char *matname;
581: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
582: PetscObjectGetName((PetscObject)A,&matname);
583: PetscObjectSetName((PetscObject)B,matname);
584: MatView(B,viewer);
585: MatDestroy(&B);
586: }
587: return(0);
588: }
590: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
591: {
592: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
593: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
594: PetscInt *ai = a->i,*ailen = a->ilen;
595: PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
596: MatScalar *ap,*aa = a->a;
599: for (k=0; k<m; k++) { /* loop over rows */
600: row = im[k]; brow = row/bs;
601: if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
602: 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);
603: rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
604: nrow = ailen[brow];
605: for (l=0; l<n; l++) { /* loop over columns */
606: if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
607: 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);
608: col = in[l];
609: bcol = col/bs;
610: cidx = col%bs;
611: ridx = row%bs;
612: high = nrow;
613: low = 0; /* assume unsorted */
614: while (high-low > 5) {
615: t = (low+high)/2;
616: if (rp[t] > bcol) high = t;
617: else low = t;
618: }
619: for (i=low; i<high; i++) {
620: if (rp[i] > bcol) break;
621: if (rp[i] == bcol) {
622: *v++ = ap[bs2*i+bs*cidx+ridx];
623: goto finished;
624: }
625: }
626: *v++ = 0.0;
627: finished:;
628: }
629: }
630: return(0);
631: }
633: PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B)
634: {
635: Mat C;
639: MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C);
640: MatPermute(C,rowp,colp,B);
641: MatDestroy(&C);
642: if (rowp == colp) {
643: MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B);
644: }
645: return(0);
646: }
648: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
649: {
650: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
651: PetscErrorCode ierr;
652: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
653: PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
654: PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
655: PetscBool roworiented=a->roworiented;
656: const PetscScalar *value = v;
657: MatScalar *ap,*aa = a->a,*bap;
660: if (roworiented) stepval = (n-1)*bs;
661: else stepval = (m-1)*bs;
663: for (k=0; k<m; k++) { /* loop over added rows */
664: row = im[k];
665: if (row < 0) continue;
666: if (PetscUnlikelyDebug(row >= a->mbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index row too large %D max %D",row,a->mbs-1);
667: rp = aj + ai[row];
668: ap = aa + bs2*ai[row];
669: rmax = imax[row];
670: nrow = ailen[row];
671: low = 0;
672: high = nrow;
673: for (l=0; l<n; l++) { /* loop over added columns */
674: if (in[l] < 0) continue;
675: col = in[l];
676: if (PetscUnlikelyDebug(col >= a->nbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index column too large %D max %D",col,a->nbs-1);
677: if (col < row) {
678: if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
679: 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)");
680: }
681: if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
682: else value = v + l*(stepval+bs)*bs + k*bs;
684: if (col <= lastcol) low = 0;
685: else high = nrow;
687: lastcol = col;
688: while (high-low > 7) {
689: t = (low+high)/2;
690: if (rp[t] > col) high = t;
691: else low = t;
692: }
693: for (i=low; i<high; i++) {
694: if (rp[i] > col) break;
695: if (rp[i] == col) {
696: bap = ap + bs2*i;
697: if (roworiented) {
698: if (is == ADD_VALUES) {
699: for (ii=0; ii<bs; ii++,value+=stepval) {
700: for (jj=ii; jj<bs2; jj+=bs) {
701: bap[jj] += *value++;
702: }
703: }
704: } else {
705: for (ii=0; ii<bs; ii++,value+=stepval) {
706: for (jj=ii; jj<bs2; jj+=bs) {
707: bap[jj] = *value++;
708: }
709: }
710: }
711: } else {
712: if (is == ADD_VALUES) {
713: for (ii=0; ii<bs; ii++,value+=stepval) {
714: for (jj=0; jj<bs; jj++) {
715: *bap++ += *value++;
716: }
717: }
718: } else {
719: for (ii=0; ii<bs; ii++,value+=stepval) {
720: for (jj=0; jj<bs; jj++) {
721: *bap++ = *value++;
722: }
723: }
724: }
725: }
726: goto noinsert2;
727: }
728: }
729: if (nonew == 1) goto noinsert2;
730: 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);
731: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
732: N = nrow++ - 1; high++;
733: /* shift up all the later entries in this row */
734: PetscArraymove(rp+i+1,rp+i,N-i+1);
735: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
736: PetscArrayzero(ap+bs2*i,bs2);
737: rp[i] = col;
738: bap = ap + bs2*i;
739: if (roworiented) {
740: for (ii=0; ii<bs; ii++,value+=stepval) {
741: for (jj=ii; jj<bs2; jj+=bs) {
742: bap[jj] = *value++;
743: }
744: }
745: } else {
746: for (ii=0; ii<bs; ii++,value+=stepval) {
747: for (jj=0; jj<bs; jj++) {
748: *bap++ = *value++;
749: }
750: }
751: }
752: noinsert2:;
753: low = i;
754: }
755: ailen[row] = nrow;
756: }
757: return(0);
758: }
760: /*
761: This is not yet used
762: */
763: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
764: {
765: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
767: const PetscInt *ai = a->i, *aj = a->j,*cols;
768: PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
769: PetscBool flag;
772: PetscMalloc1(m,&ns);
773: while (i < m) {
774: nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
775: /* Limits the number of elements in a node to 'a->inode.limit' */
776: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
777: nzy = ai[j+1] - ai[j];
778: if (nzy != (nzx - j + i)) break;
779: PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);
780: if (!flag) break;
781: }
782: ns[node_count++] = blk_size;
784: i = j;
785: }
786: if (!a->inode.size && m && node_count > .9*m) {
787: PetscFree(ns);
788: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
789: } else {
790: a->inode.node_count = node_count;
792: PetscMalloc1(node_count,&a->inode.size);
793: PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));
794: PetscArraycpy(a->inode.size,ns,node_count);
795: PetscFree(ns);
796: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
798: /* count collections of adjacent columns in each inode */
799: row = 0;
800: cnt = 0;
801: for (i=0; i<node_count; i++) {
802: cols = aj + ai[row] + a->inode.size[i];
803: nz = ai[row+1] - ai[row] - a->inode.size[i];
804: for (j=1; j<nz; j++) {
805: if (cols[j] != cols[j-1]+1) cnt++;
806: }
807: cnt++;
808: row += a->inode.size[i];
809: }
810: PetscMalloc1(2*cnt,&counts);
811: cnt = 0;
812: row = 0;
813: for (i=0; i<node_count; i++) {
814: cols = aj + ai[row] + a->inode.size[i];
815: counts[2*cnt] = cols[0];
816: nz = ai[row+1] - ai[row] - a->inode.size[i];
817: cnt2 = 1;
818: for (j=1; j<nz; j++) {
819: if (cols[j] != cols[j-1]+1) {
820: counts[2*(cnt++)+1] = cnt2;
821: counts[2*cnt] = cols[j];
822: cnt2 = 1;
823: } else cnt2++;
824: }
825: counts[2*(cnt++)+1] = cnt2;
826: row += a->inode.size[i];
827: }
828: PetscIntView(2*cnt,counts,NULL);
829: }
830: return(0);
831: }
833: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
834: {
835: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
837: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
838: PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
839: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
840: MatScalar *aa = a->a,*ap;
843: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
845: if (m) rmax = ailen[0];
846: for (i=1; i<mbs; i++) {
847: /* move each row back by the amount of empty slots (fshift) before it*/
848: fshift += imax[i-1] - ailen[i-1];
849: rmax = PetscMax(rmax,ailen[i]);
850: if (fshift) {
851: ip = aj + ai[i];
852: ap = aa + bs2*ai[i];
853: N = ailen[i];
854: PetscArraymove(ip-fshift,ip,N);
855: PetscArraymove(ap-bs2*fshift,ap,bs2*N);
856: }
857: ai[i] = ai[i-1] + ailen[i-1];
858: }
859: if (mbs) {
860: fshift += imax[mbs-1] - ailen[mbs-1];
861: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
862: }
863: /* reset ilen and imax for each row */
864: for (i=0; i<mbs; i++) {
865: ailen[i] = imax[i] = ai[i+1] - ai[i];
866: }
867: a->nz = ai[mbs];
869: /* diagonals may have moved, reset it */
870: if (a->diag) {
871: PetscArraycpy(a->diag,ai,mbs);
872: }
873: 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);
875: 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);
876: PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
877: PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);
879: A->info.mallocs += a->reallocs;
880: a->reallocs = 0;
881: A->info.nz_unneeded = (PetscReal)fshift*bs2;
882: a->idiagvalid = PETSC_FALSE;
883: a->rmax = rmax;
885: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
886: if (a->jshort && a->free_jshort) {
887: /* when matrix data structure is changed, previous jshort must be replaced */
888: PetscFree(a->jshort);
889: }
890: PetscMalloc1(a->i[A->rmap->n],&a->jshort);
891: PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));
892: for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
893: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
894: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
895: a->free_jshort = PETSC_TRUE;
896: }
897: return(0);
898: }
900: /*
901: This function returns an array of flags which indicate the locations of contiguous
902: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
903: then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)]
904: Assume: sizes should be long enough to hold all the values.
905: */
906: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
907: {
908: PetscInt i,j,k,row;
909: PetscBool flg;
912: for (i=0,j=0; i<n; j++) {
913: row = idx[i];
914: if (row%bs!=0) { /* Not the beginning of a block */
915: sizes[j] = 1;
916: i++;
917: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
918: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
919: i++;
920: } else { /* Begining of the block, so check if the complete block exists */
921: flg = PETSC_TRUE;
922: for (k=1; k<bs; k++) {
923: if (row+k != idx[i+k]) { /* break in the block */
924: flg = PETSC_FALSE;
925: break;
926: }
927: }
928: if (flg) { /* No break in the bs */
929: sizes[j] = bs;
930: i += bs;
931: } else {
932: sizes[j] = 1;
933: i++;
934: }
935: }
936: }
937: *bs_max = j;
938: return(0);
939: }
941: /* Only add/insert a(i,j) with i<=j (blocks).
942: Any a(i,j) with i>j input by user is ingored.
943: */
945: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
946: {
947: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
949: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
950: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
951: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
952: PetscInt ridx,cidx,bs2=a->bs2;
953: MatScalar *ap,value,*aa=a->a,*bap;
956: for (k=0; k<m; k++) { /* loop over added rows */
957: row = im[k]; /* row number */
958: brow = row/bs; /* block row number */
959: if (row < 0) continue;
960: if (PetscUnlikelyDebug(row >= A->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1);
961: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
962: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
963: rmax = imax[brow]; /* maximum space allocated for this row */
964: nrow = ailen[brow]; /* actual length of this row */
965: low = 0;
966: high = nrow;
967: for (l=0; l<n; l++) { /* loop over added columns */
968: if (in[l] < 0) continue;
969: if (PetscUnlikelyDebug(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);
970: col = in[l];
971: bcol = col/bs; /* block col number */
973: if (brow > bcol) {
974: if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
975: 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)");
976: }
978: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
979: if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
980: /* element value a(k,l) */
981: if (roworiented) value = v[l + k*n];
982: else value = v[k + l*m];
984: /* move pointer bap to a(k,l) quickly and add/insert value */
985: if (col <= lastcol) low = 0;
986: else high = nrow;
988: lastcol = col;
989: while (high-low > 7) {
990: t = (low+high)/2;
991: if (rp[t] > bcol) high = t;
992: else low = t;
993: }
994: for (i=low; i<high; i++) {
995: if (rp[i] > bcol) break;
996: if (rp[i] == bcol) {
997: bap = ap + bs2*i + bs*cidx + ridx;
998: if (is == ADD_VALUES) *bap += value;
999: else *bap = value;
1000: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1001: if (brow == bcol && ridx < cidx) {
1002: bap = ap + bs2*i + bs*ridx + cidx;
1003: if (is == ADD_VALUES) *bap += value;
1004: else *bap = value;
1005: }
1006: goto noinsert1;
1007: }
1008: }
1010: if (nonew == 1) goto noinsert1;
1011: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
1012: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
1014: N = nrow++ - 1; high++;
1015: /* shift up all the later entries in this row */
1016: PetscArraymove(rp+i+1,rp+i,N-i+1);
1017: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
1018: PetscArrayzero(ap+bs2*i,bs2);
1019: rp[i] = bcol;
1020: ap[bs2*i + bs*cidx + ridx] = value;
1021: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1022: if (brow == bcol && ridx < cidx) {
1023: ap[bs2*i + bs*ridx + cidx] = value;
1024: }
1025: A->nonzerostate++;
1026: noinsert1:;
1027: low = i;
1028: }
1029: } /* end of loop over added columns */
1030: ailen[brow] = nrow;
1031: } /* end of loop over added rows */
1032: return(0);
1033: }
1035: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
1036: {
1037: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1038: Mat outA;
1040: PetscBool row_identity;
1043: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
1044: ISIdentity(row,&row_identity);
1045: if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1046: 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()! */
1048: outA = inA;
1049: inA->factortype = MAT_FACTOR_ICC;
1050: PetscFree(inA->solvertype);
1051: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
1053: MatMarkDiagonal_SeqSBAIJ(inA);
1054: MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);
1056: PetscObjectReference((PetscObject)row);
1057: ISDestroy(&a->row);
1058: a->row = row;
1059: PetscObjectReference((PetscObject)row);
1060: ISDestroy(&a->col);
1061: a->col = row;
1063: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1064: if (a->icol) {ISInvertPermutation(row,PETSC_DECIDE, &a->icol);}
1065: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
1067: if (!a->solve_work) {
1068: PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
1069: PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
1070: }
1072: MatCholeskyFactorNumeric(outA,inA,info);
1073: return(0);
1074: }
1076: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1077: {
1078: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1079: PetscInt i,nz,n;
1083: nz = baij->maxnz;
1084: n = mat->cmap->n;
1085: for (i=0; i<nz; i++) baij->j[i] = indices[i];
1087: baij->nz = nz;
1088: for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1090: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1091: return(0);
1092: }
1094: /*@
1095: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1096: in the matrix.
1098: Input Parameters:
1099: + mat - the SeqSBAIJ matrix
1100: - indices - the column indices
1102: Level: advanced
1104: Notes:
1105: This can be called if you have precomputed the nonzero structure of the
1106: matrix and want to provide it to the matrix object to improve the performance
1107: of the MatSetValues() operation.
1109: You MUST have set the correct numbers of nonzeros per row in the call to
1110: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1112: MUST be called before any calls to MatSetValues()
1114: .seealso: MatCreateSeqSBAIJ
1115: @*/
1116: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1117: {
1123: PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
1124: return(0);
1125: }
1127: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1128: {
1130: PetscBool isbaij;
1133: PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");
1134: if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name);
1135: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1136: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1137: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1138: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1140: 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");
1141: if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different");
1142: if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size");
1143: PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);
1144: PetscObjectStateIncrease((PetscObject)B);
1145: } else {
1146: MatGetRowUpperTriangular(A);
1147: MatCopy_Basic(A,B,str);
1148: MatRestoreRowUpperTriangular(A);
1149: }
1150: return(0);
1151: }
1153: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1154: {
1158: MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);
1159: return(0);
1160: }
1162: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1163: {
1164: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1167: *array = a->a;
1168: return(0);
1169: }
1171: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1172: {
1174: *array = NULL;
1175: return(0);
1176: }
1178: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1179: {
1180: PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1181: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1182: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1186: /* Set the number of nonzeros in the new matrix */
1187: MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
1188: return(0);
1189: }
1191: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1192: {
1193: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1195: PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1196: PetscBLASInt one = 1;
1199: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
1200: PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE;
1201: if (e) {
1202: PetscArraycmp(x->i,y->i,x->mbs+1,&e);
1203: if (e) {
1204: PetscArraycmp(x->j,y->j,x->i[x->mbs],&e);
1205: if (e) str = SAME_NONZERO_PATTERN;
1206: }
1207: }
1208: if (!e && str == SAME_NONZERO_PATTERN) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN");
1209: }
1210: if (str == SAME_NONZERO_PATTERN) {
1211: PetscScalar alpha = a;
1212: PetscBLASInt bnz;
1213: PetscBLASIntCast(x->nz*bs2,&bnz);
1214: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1215: PetscObjectStateIncrease((PetscObject)Y);
1216: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1217: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);
1218: MatAXPY_Basic(Y,a,X,str);
1219: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);
1220: } else {
1221: Mat B;
1222: PetscInt *nnz;
1223: if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
1224: MatGetRowUpperTriangular(X);
1225: MatGetRowUpperTriangular(Y);
1226: PetscMalloc1(Y->rmap->N,&nnz);
1227: MatCreate(PetscObjectComm((PetscObject)Y),&B);
1228: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
1229: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
1230: MatSetBlockSizesFromMats(B,Y,Y);
1231: MatSetType(B,((PetscObject)Y)->type_name);
1232: MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);
1233: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1235: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1237: MatHeaderReplace(Y,&B);
1238: PetscFree(nnz);
1239: MatRestoreRowUpperTriangular(X);
1240: MatRestoreRowUpperTriangular(Y);
1241: }
1242: return(0);
1243: }
1245: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1246: {
1248: *flg = PETSC_TRUE;
1249: return(0);
1250: }
1252: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1253: {
1255: *flg = PETSC_TRUE;
1256: return(0);
1257: }
1259: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1260: {
1262: *flg = PETSC_FALSE;
1263: return(0);
1264: }
1266: PetscErrorCode MatConjugate_SeqSBAIJ(Mat A)
1267: {
1268: #if defined(PETSC_USE_COMPLEX)
1269: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1270: PetscInt i,nz = a->bs2*a->i[a->mbs];
1271: MatScalar *aa = a->a;
1274: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1275: #else
1277: #endif
1278: return(0);
1279: }
1281: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1282: {
1283: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1284: PetscInt i,nz = a->bs2*a->i[a->mbs];
1285: MatScalar *aa = a->a;
1288: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1289: return(0);
1290: }
1292: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1293: {
1294: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1295: PetscInt i,nz = a->bs2*a->i[a->mbs];
1296: MatScalar *aa = a->a;
1299: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1300: return(0);
1301: }
1303: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1304: {
1305: Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1306: PetscErrorCode ierr;
1307: PetscInt i,j,k,count;
1308: PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1309: PetscScalar zero = 0.0;
1310: MatScalar *aa;
1311: const PetscScalar *xx;
1312: PetscScalar *bb;
1313: PetscBool *zeroed,vecs = PETSC_FALSE;
1316: /* fix right hand side if needed */
1317: if (x && b) {
1318: VecGetArrayRead(x,&xx);
1319: VecGetArray(b,&bb);
1320: vecs = PETSC_TRUE;
1321: }
1323: /* zero the columns */
1324: PetscCalloc1(A->rmap->n,&zeroed);
1325: for (i=0; i<is_n; i++) {
1326: 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]);
1327: zeroed[is_idx[i]] = PETSC_TRUE;
1328: }
1329: if (vecs) {
1330: for (i=0; i<A->rmap->N; i++) {
1331: row = i/bs;
1332: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1333: for (k=0; k<bs; k++) {
1334: col = bs*baij->j[j] + k;
1335: if (col <= i) continue;
1336: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1337: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1338: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1339: }
1340: }
1341: }
1342: for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1343: }
1345: for (i=0; i<A->rmap->N; i++) {
1346: if (!zeroed[i]) {
1347: row = i/bs;
1348: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1349: for (k=0; k<bs; k++) {
1350: col = bs*baij->j[j] + k;
1351: if (zeroed[col]) {
1352: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1353: aa[0] = 0.0;
1354: }
1355: }
1356: }
1357: }
1358: }
1359: PetscFree(zeroed);
1360: if (vecs) {
1361: VecRestoreArrayRead(x,&xx);
1362: VecRestoreArray(b,&bb);
1363: }
1365: /* zero the rows */
1366: for (i=0; i<is_n; i++) {
1367: row = is_idx[i];
1368: count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1369: aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1370: for (k=0; k<count; k++) {
1371: aa[0] = zero;
1372: aa += bs;
1373: }
1374: if (diag != 0.0) {
1375: (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
1376: }
1377: }
1378: MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);
1379: return(0);
1380: }
1382: PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1383: {
1385: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1388: if (!Y->preallocated || !aij->nz) {
1389: MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
1390: }
1391: MatShift_Basic(Y,a);
1392: return(0);
1393: }
1395: /* -------------------------------------------------------------------*/
1396: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1397: MatGetRow_SeqSBAIJ,
1398: MatRestoreRow_SeqSBAIJ,
1399: MatMult_SeqSBAIJ_N,
1400: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1401: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1402: MatMultAdd_SeqSBAIJ_N,
1403: NULL,
1404: NULL,
1405: NULL,
1406: /* 10*/ NULL,
1407: NULL,
1408: MatCholeskyFactor_SeqSBAIJ,
1409: MatSOR_SeqSBAIJ,
1410: MatTranspose_SeqSBAIJ,
1411: /* 15*/ MatGetInfo_SeqSBAIJ,
1412: MatEqual_SeqSBAIJ,
1413: MatGetDiagonal_SeqSBAIJ,
1414: MatDiagonalScale_SeqSBAIJ,
1415: MatNorm_SeqSBAIJ,
1416: /* 20*/ NULL,
1417: MatAssemblyEnd_SeqSBAIJ,
1418: MatSetOption_SeqSBAIJ,
1419: MatZeroEntries_SeqSBAIJ,
1420: /* 24*/ NULL,
1421: NULL,
1422: NULL,
1423: NULL,
1424: NULL,
1425: /* 29*/ MatSetUp_SeqSBAIJ,
1426: NULL,
1427: NULL,
1428: NULL,
1429: NULL,
1430: /* 34*/ MatDuplicate_SeqSBAIJ,
1431: NULL,
1432: NULL,
1433: NULL,
1434: MatICCFactor_SeqSBAIJ,
1435: /* 39*/ MatAXPY_SeqSBAIJ,
1436: MatCreateSubMatrices_SeqSBAIJ,
1437: MatIncreaseOverlap_SeqSBAIJ,
1438: MatGetValues_SeqSBAIJ,
1439: MatCopy_SeqSBAIJ,
1440: /* 44*/ NULL,
1441: MatScale_SeqSBAIJ,
1442: MatShift_SeqSBAIJ,
1443: NULL,
1444: MatZeroRowsColumns_SeqSBAIJ,
1445: /* 49*/ NULL,
1446: MatGetRowIJ_SeqSBAIJ,
1447: MatRestoreRowIJ_SeqSBAIJ,
1448: NULL,
1449: NULL,
1450: /* 54*/ NULL,
1451: NULL,
1452: NULL,
1453: MatPermute_SeqSBAIJ,
1454: MatSetValuesBlocked_SeqSBAIJ,
1455: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1456: NULL,
1457: NULL,
1458: NULL,
1459: NULL,
1460: /* 64*/ NULL,
1461: NULL,
1462: NULL,
1463: NULL,
1464: NULL,
1465: /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1466: NULL,
1467: MatConvert_MPISBAIJ_Basic,
1468: NULL,
1469: NULL,
1470: /* 74*/ NULL,
1471: NULL,
1472: NULL,
1473: NULL,
1474: NULL,
1475: /* 79*/ NULL,
1476: NULL,
1477: NULL,
1478: MatGetInertia_SeqSBAIJ,
1479: MatLoad_SeqSBAIJ,
1480: /* 84*/ MatIsSymmetric_SeqSBAIJ,
1481: MatIsHermitian_SeqSBAIJ,
1482: MatIsStructurallySymmetric_SeqSBAIJ,
1483: NULL,
1484: NULL,
1485: /* 89*/ NULL,
1486: NULL,
1487: NULL,
1488: NULL,
1489: NULL,
1490: /* 94*/ NULL,
1491: NULL,
1492: NULL,
1493: NULL,
1494: NULL,
1495: /* 99*/ NULL,
1496: NULL,
1497: NULL,
1498: MatConjugate_SeqSBAIJ,
1499: NULL,
1500: /*104*/ NULL,
1501: MatRealPart_SeqSBAIJ,
1502: MatImaginaryPart_SeqSBAIJ,
1503: MatGetRowUpperTriangular_SeqSBAIJ,
1504: MatRestoreRowUpperTriangular_SeqSBAIJ,
1505: /*109*/ NULL,
1506: NULL,
1507: NULL,
1508: NULL,
1509: MatMissingDiagonal_SeqSBAIJ,
1510: /*114*/ NULL,
1511: NULL,
1512: NULL,
1513: NULL,
1514: NULL,
1515: /*119*/ NULL,
1516: NULL,
1517: NULL,
1518: NULL,
1519: NULL,
1520: /*124*/ NULL,
1521: NULL,
1522: NULL,
1523: NULL,
1524: NULL,
1525: /*129*/ NULL,
1526: NULL,
1527: NULL,
1528: NULL,
1529: NULL,
1530: /*134*/ NULL,
1531: NULL,
1532: NULL,
1533: NULL,
1534: NULL,
1535: /*139*/ MatSetBlockSizes_Default,
1536: NULL,
1537: NULL,
1538: NULL,
1539: NULL,
1540: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ
1541: };
1543: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1544: {
1545: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1546: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1550: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1552: /* allocate space for values if not already there */
1553: if (!aij->saved_values) {
1554: PetscMalloc1(nz+1,&aij->saved_values);
1555: }
1557: /* copy values over */
1558: PetscArraycpy(aij->saved_values,aij->a,nz);
1559: return(0);
1560: }
1562: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1563: {
1564: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1566: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1569: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1570: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1572: /* copy values over */
1573: PetscArraycpy(aij->a,aij->saved_values,nz);
1574: return(0);
1575: }
1577: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1578: {
1579: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1581: PetscInt i,mbs,nbs,bs2;
1582: PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1585: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1587: MatSetBlockSize(B,PetscAbs(bs));
1588: PetscLayoutSetUp(B->rmap);
1589: PetscLayoutSetUp(B->cmap);
1590: 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);
1591: PetscLayoutGetBlockSize(B->rmap,&bs);
1593: B->preallocated = PETSC_TRUE;
1595: mbs = B->rmap->N/bs;
1596: nbs = B->cmap->n/bs;
1597: bs2 = bs*bs;
1599: 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");
1601: if (nz == MAT_SKIP_ALLOCATION) {
1602: skipallocation = PETSC_TRUE;
1603: nz = 0;
1604: }
1606: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1607: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1608: if (nnz) {
1609: for (i=0; i<mbs; i++) {
1610: 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]);
1611: 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);
1612: }
1613: }
1615: B->ops->mult = MatMult_SeqSBAIJ_N;
1616: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1617: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1618: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1620: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1621: if (!flg) {
1622: switch (bs) {
1623: case 1:
1624: B->ops->mult = MatMult_SeqSBAIJ_1;
1625: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1626: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1627: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1628: break;
1629: case 2:
1630: B->ops->mult = MatMult_SeqSBAIJ_2;
1631: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1632: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1633: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1634: break;
1635: case 3:
1636: B->ops->mult = MatMult_SeqSBAIJ_3;
1637: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1638: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1639: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1640: break;
1641: case 4:
1642: B->ops->mult = MatMult_SeqSBAIJ_4;
1643: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1644: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1645: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1646: break;
1647: case 5:
1648: B->ops->mult = MatMult_SeqSBAIJ_5;
1649: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1650: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1651: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1652: break;
1653: case 6:
1654: B->ops->mult = MatMult_SeqSBAIJ_6;
1655: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1656: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1657: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1658: break;
1659: case 7:
1660: B->ops->mult = MatMult_SeqSBAIJ_7;
1661: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1662: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1663: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1664: break;
1665: }
1666: }
1668: b->mbs = mbs;
1669: b->nbs = nbs;
1670: if (!skipallocation) {
1671: if (!b->imax) {
1672: PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
1674: b->free_imax_ilen = PETSC_TRUE;
1676: PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));
1677: }
1678: if (!nnz) {
1679: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1680: else if (nz <= 0) nz = 1;
1681: nz = PetscMin(nbs,nz);
1682: for (i=0; i<mbs; i++) b->imax[i] = nz;
1683: PetscIntMultError(nz,mbs,&nz);
1684: } else {
1685: PetscInt64 nz64 = 0;
1686: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1687: PetscIntCast(nz64,&nz);
1688: }
1689: /* b->ilen will count nonzeros in each block row so far. */
1690: for (i=0; i<mbs; i++) b->ilen[i] = 0;
1691: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1693: /* allocate the matrix space */
1694: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
1695: PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
1696: PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
1697: PetscArrayzero(b->a,nz*bs2);
1698: PetscArrayzero(b->j,nz);
1700: b->singlemalloc = PETSC_TRUE;
1702: /* pointer to beginning of each row */
1703: b->i[0] = 0;
1704: for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1706: b->free_a = PETSC_TRUE;
1707: b->free_ij = PETSC_TRUE;
1708: } else {
1709: b->free_a = PETSC_FALSE;
1710: b->free_ij = PETSC_FALSE;
1711: }
1713: b->bs2 = bs2;
1714: b->nz = 0;
1715: b->maxnz = nz;
1716: b->inew = NULL;
1717: b->jnew = NULL;
1718: b->anew = NULL;
1719: b->a2anew = NULL;
1720: b->permute = PETSC_FALSE;
1722: B->was_assembled = PETSC_FALSE;
1723: B->assembled = PETSC_FALSE;
1724: if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
1725: return(0);
1726: }
1728: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1729: {
1730: PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1731: PetscScalar *values=NULL;
1732: PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1736: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
1737: PetscLayoutSetBlockSize(B->rmap,bs);
1738: PetscLayoutSetBlockSize(B->cmap,bs);
1739: PetscLayoutSetUp(B->rmap);
1740: PetscLayoutSetUp(B->cmap);
1741: PetscLayoutGetBlockSize(B->rmap,&bs);
1742: m = B->rmap->n/bs;
1744: if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]);
1745: PetscMalloc1(m+1,&nnz);
1746: for (i=0; i<m; i++) {
1747: nz = ii[i+1] - ii[i];
1748: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz);
1749: anz = 0;
1750: for (j=0; j<nz; j++) {
1751: /* count only values on the diagonal or above */
1752: if (jj[ii[i] + j] >= i) {
1753: anz = nz - j;
1754: break;
1755: }
1756: }
1757: nz_max = PetscMax(nz_max,anz);
1758: nnz[i] = anz;
1759: }
1760: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1761: PetscFree(nnz);
1763: values = (PetscScalar*)V;
1764: if (!values) {
1765: PetscCalloc1(bs*bs*nz_max,&values);
1766: }
1767: for (i=0; i<m; i++) {
1768: PetscInt ncols = ii[i+1] - ii[i];
1769: const PetscInt *icols = jj + ii[i];
1770: if (!roworiented || bs == 1) {
1771: const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1772: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
1773: } else {
1774: for (j=0; j<ncols; j++) {
1775: const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1776: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
1777: }
1778: }
1779: }
1780: if (!V) { PetscFree(values); }
1781: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1782: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1783: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1784: return(0);
1785: }
1787: /*
1788: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1789: */
1790: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1791: {
1793: PetscBool flg = PETSC_FALSE;
1794: PetscInt bs = B->rmap->bs;
1797: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1798: if (flg) bs = 8;
1800: if (!natural) {
1801: switch (bs) {
1802: case 1:
1803: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1804: break;
1805: case 2:
1806: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1807: break;
1808: case 3:
1809: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1810: break;
1811: case 4:
1812: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1813: break;
1814: case 5:
1815: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1816: break;
1817: case 6:
1818: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1819: break;
1820: case 7:
1821: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1822: break;
1823: default:
1824: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1825: break;
1826: }
1827: } else {
1828: switch (bs) {
1829: case 1:
1830: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1831: break;
1832: case 2:
1833: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1834: break;
1835: case 3:
1836: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1837: break;
1838: case 4:
1839: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1840: break;
1841: case 5:
1842: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1843: break;
1844: case 6:
1845: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1846: break;
1847: case 7:
1848: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1849: break;
1850: default:
1851: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1852: break;
1853: }
1854: }
1855: return(0);
1856: }
1858: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*);
1859: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*);
1860: static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type)
1861: {
1863: *type = MATSOLVERPETSC;
1864: return(0);
1865: }
1867: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1868: {
1869: PetscInt n = A->rmap->n;
1873: #if defined(PETSC_USE_COMPLEX)
1874: if (A->hermitian && !A->symmetric && (ftype == MAT_FACTOR_CHOLESKY||ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported");
1875: #endif
1877: MatCreate(PetscObjectComm((PetscObject)A),B);
1878: MatSetSizes(*B,n,n,n,n);
1879: if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1880: MatSetType(*B,MATSEQSBAIJ);
1881: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
1883: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1884: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1885: PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
1886: PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
1887: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1889: (*B)->factortype = ftype;
1890: (*B)->canuseordering = PETSC_TRUE;
1891: PetscFree((*B)->solvertype);
1892: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
1893: PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);
1894: return(0);
1895: }
1897: /*@C
1898: MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1900: Not Collective
1902: Input Parameter:
1903: . mat - a MATSEQSBAIJ matrix
1905: Output Parameter:
1906: . array - pointer to the data
1908: Level: intermediate
1910: .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1911: @*/
1912: PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1913: {
1917: PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
1918: return(0);
1919: }
1921: /*@C
1922: MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1924: Not Collective
1926: Input Parameters:
1927: + mat - a MATSEQSBAIJ matrix
1928: - array - pointer to the data
1930: Level: intermediate
1932: .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1933: @*/
1934: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1935: {
1939: PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
1940: return(0);
1941: }
1943: /*MC
1944: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1945: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1947: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1948: can call MatSetOption(Mat, MAT_HERMITIAN).
1950: Options Database Keys:
1951: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1953: Notes:
1954: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1955: stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1956: 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.
1958: The number of rows in the matrix must be less than or equal to the number of columns
1960: Level: beginner
1962: .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1963: M*/
1964: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1965: {
1966: Mat_SeqSBAIJ *b;
1968: PetscMPIInt size;
1969: PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1972: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1973: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1975: PetscNewLog(B,&b);
1976: B->data = (void*)b;
1977: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1979: B->ops->destroy = MatDestroy_SeqSBAIJ;
1980: B->ops->view = MatView_SeqSBAIJ;
1981: b->row = NULL;
1982: b->icol = NULL;
1983: b->reallocs = 0;
1984: b->saved_values = NULL;
1985: b->inode.limit = 5;
1986: b->inode.max_limit = 5;
1988: b->roworiented = PETSC_TRUE;
1989: b->nonew = 0;
1990: b->diag = NULL;
1991: b->solve_work = NULL;
1992: b->mult_work = NULL;
1993: B->spptr = NULL;
1994: B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1995: b->keepnonzeropattern = PETSC_FALSE;
1997: b->inew = NULL;
1998: b->jnew = NULL;
1999: b->anew = NULL;
2000: b->a2anew = NULL;
2001: b->permute = PETSC_FALSE;
2003: b->ignore_ltriangular = PETSC_TRUE;
2005: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);
2007: b->getrow_utriangular = PETSC_FALSE;
2009: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);
2011: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);
2012: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);
2013: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);
2014: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);
2015: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
2016: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);
2017: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);
2018: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);
2019: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
2020: #if defined(PETSC_HAVE_ELEMENTAL)
2021: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);
2022: #endif
2023: #if defined(PETSC_HAVE_SCALAPACK)
2024: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);
2025: #endif
2027: B->symmetric = PETSC_TRUE;
2028: B->structurally_symmetric = PETSC_TRUE;
2029: B->symmetric_set = PETSC_TRUE;
2030: B->structurally_symmetric_set = PETSC_TRUE;
2031: B->symmetric_eternal = PETSC_TRUE;
2032: #if defined(PETSC_USE_COMPLEX)
2033: B->hermitian = PETSC_FALSE;
2034: B->hermitian_set = PETSC_FALSE;
2035: #else
2036: B->hermitian = PETSC_TRUE;
2037: B->hermitian_set = PETSC_TRUE;
2038: #endif
2040: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
2042: PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");
2043: PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);
2044: if (no_unroll) {
2045: PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");
2046: }
2047: PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);
2048: if (no_inode) {
2049: PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");
2050: }
2051: PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);
2052: PetscOptionsEnd();
2053: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
2054: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
2055: return(0);
2056: }
2058: /*@C
2059: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
2060: compressed row) format. For good matrix assembly performance the
2061: user should preallocate the matrix storage by setting the parameter nz
2062: (or the array nnz). By setting these parameters accurately, performance
2063: during matrix assembly can be increased by more than a factor of 50.
2065: Collective on Mat
2067: Input Parameters:
2068: + B - the symmetric matrix
2069: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2070: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2071: . nz - number of block nonzeros per block row (same for all rows)
2072: - nnz - array containing the number of block nonzeros in the upper triangular plus
2073: diagonal portion of each block (possibly different for each block row) or NULL
2075: Options Database Keys:
2076: + -mat_no_unroll - uses code that does not unroll the loops in the
2077: block calculations (much slower)
2078: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2080: Level: intermediate
2082: Notes:
2083: Specify the preallocated storage with either nz or nnz (not both).
2084: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2085: allocation. See Users-Manual: ch_mat for details.
2087: You can call MatGetInfo() to get information on how effective the preallocation was;
2088: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2089: You can also run with the option -info and look for messages with the string
2090: malloc in them to see if additional memory allocation was needed.
2092: If the nnz parameter is given then the nz parameter is ignored
2094: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2095: @*/
2096: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2097: {
2104: PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
2105: return(0);
2106: }
2108: /*@C
2109: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2111: Input Parameters:
2112: + B - the matrix
2113: . bs - size of block, the blocks are ALWAYS square.
2114: . i - the indices into j for the start of each local row (starts with zero)
2115: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2116: - v - optional values in the matrix
2118: Level: advanced
2120: Notes:
2121: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2122: may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2123: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2124: MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2125: block column and the second index is over columns within a block.
2127: Any entries below the diagonal are ignored
2129: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2130: and usually the numerical values as well
2132: .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2133: @*/
2134: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2135: {
2142: PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
2143: return(0);
2144: }
2146: /*@C
2147: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2148: compressed row) format. For good matrix assembly performance the
2149: user should preallocate the matrix storage by setting the parameter nz
2150: (or the array nnz). By setting these parameters accurately, performance
2151: during matrix assembly can be increased by more than a factor of 50.
2153: Collective
2155: Input Parameters:
2156: + comm - MPI communicator, set to PETSC_COMM_SELF
2157: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2158: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2159: . m - number of rows, or number of columns
2160: . nz - number of block nonzeros per block row (same for all rows)
2161: - nnz - array containing the number of block nonzeros in the upper triangular plus
2162: diagonal portion of each block (possibly different for each block row) or NULL
2164: Output Parameter:
2165: . A - the symmetric matrix
2167: Options Database Keys:
2168: + -mat_no_unroll - uses code that does not unroll the loops in the
2169: block calculations (much slower)
2170: - -mat_block_size - size of the blocks to use
2172: Level: intermediate
2174: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2175: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2176: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2178: Notes:
2179: The number of rows and columns must be divisible by blocksize.
2180: This matrix type does not support complex Hermitian operation.
2182: Specify the preallocated storage with either nz or nnz (not both).
2183: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2184: allocation. See Users-Manual: ch_mat for details.
2186: If the nnz parameter is given then the nz parameter is ignored
2188: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2189: @*/
2190: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2191: {
2195: MatCreate(comm,A);
2196: MatSetSizes(*A,m,n,m,n);
2197: MatSetType(*A,MATSEQSBAIJ);
2198: MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
2199: return(0);
2200: }
2202: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2203: {
2204: Mat C;
2205: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2207: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2210: if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");
2212: *B = NULL;
2213: MatCreate(PetscObjectComm((PetscObject)A),&C);
2214: MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
2215: MatSetBlockSizesFromMats(C,A,A);
2216: MatSetType(C,MATSEQSBAIJ);
2217: c = (Mat_SeqSBAIJ*)C->data;
2219: C->preallocated = PETSC_TRUE;
2220: C->factortype = A->factortype;
2221: c->row = NULL;
2222: c->icol = NULL;
2223: c->saved_values = NULL;
2224: c->keepnonzeropattern = a->keepnonzeropattern;
2225: C->assembled = PETSC_TRUE;
2227: PetscLayoutReference(A->rmap,&C->rmap);
2228: PetscLayoutReference(A->cmap,&C->cmap);
2229: c->bs2 = a->bs2;
2230: c->mbs = a->mbs;
2231: c->nbs = a->nbs;
2233: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2234: c->imax = a->imax;
2235: c->ilen = a->ilen;
2236: c->free_imax_ilen = PETSC_FALSE;
2237: } else {
2238: PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);
2239: PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));
2240: for (i=0; i<mbs; i++) {
2241: c->imax[i] = a->imax[i];
2242: c->ilen[i] = a->ilen[i];
2243: }
2244: c->free_imax_ilen = PETSC_TRUE;
2245: }
2247: /* allocate the matrix space */
2248: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2249: PetscMalloc1(bs2*nz,&c->a);
2250: PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));
2251: c->i = a->i;
2252: c->j = a->j;
2253: c->singlemalloc = PETSC_FALSE;
2254: c->free_a = PETSC_TRUE;
2255: c->free_ij = PETSC_FALSE;
2256: c->parent = A;
2257: PetscObjectReference((PetscObject)A);
2258: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2259: MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2260: } else {
2261: PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
2262: PetscArraycpy(c->i,a->i,mbs+1);
2263: PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
2264: c->singlemalloc = PETSC_TRUE;
2265: c->free_a = PETSC_TRUE;
2266: c->free_ij = PETSC_TRUE;
2267: }
2268: if (mbs > 0) {
2269: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2270: PetscArraycpy(c->j,a->j,nz);
2271: }
2272: if (cpvalues == MAT_COPY_VALUES) {
2273: PetscArraycpy(c->a,a->a,bs2*nz);
2274: } else {
2275: PetscArrayzero(c->a,bs2*nz);
2276: }
2277: if (a->jshort) {
2278: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2279: /* if the parent matrix is reassembled, this child matrix will never notice */
2280: PetscMalloc1(nz,&c->jshort);
2281: PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));
2282: PetscArraycpy(c->jshort,a->jshort,nz);
2284: c->free_jshort = PETSC_TRUE;
2285: }
2286: }
2288: c->roworiented = a->roworiented;
2289: c->nonew = a->nonew;
2291: if (a->diag) {
2292: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2293: c->diag = a->diag;
2294: c->free_diag = PETSC_FALSE;
2295: } else {
2296: PetscMalloc1(mbs,&c->diag);
2297: PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));
2298: for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2299: c->free_diag = PETSC_TRUE;
2300: }
2301: }
2302: c->nz = a->nz;
2303: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2304: c->solve_work = NULL;
2305: c->mult_work = NULL;
2307: *B = C;
2308: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
2309: return(0);
2310: }
2312: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2313: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2315: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)
2316: {
2318: PetscBool isbinary;
2321: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2322: if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)mat)->type_name);
2323: MatLoad_SeqSBAIJ_Binary(mat,viewer);
2324: return(0);
2325: }
2327: /*@
2328: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2329: (upper triangular entries in CSR format) provided by the user.
2331: Collective
2333: Input Parameters:
2334: + comm - must be an MPI communicator of size 1
2335: . bs - size of block
2336: . m - number of rows
2337: . n - number of columns
2338: . 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
2339: . j - column indices
2340: - a - matrix values
2342: Output Parameter:
2343: . mat - the matrix
2345: Level: advanced
2347: Notes:
2348: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2349: once the matrix is destroyed
2351: You cannot set new nonzero locations into this matrix, that will generate an error.
2353: The i and j indices are 0 based
2355: 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
2356: it is the regular CSR format excluding the lower triangular elements.
2358: .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2360: @*/
2361: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2362: {
2364: PetscInt ii;
2365: Mat_SeqSBAIJ *sbaij;
2368: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2369: if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2371: MatCreate(comm,mat);
2372: MatSetSizes(*mat,m,n,m,n);
2373: MatSetType(*mat,MATSEQSBAIJ);
2374: MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);
2375: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2376: PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);
2377: PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));
2379: sbaij->i = i;
2380: sbaij->j = j;
2381: sbaij->a = a;
2383: sbaij->singlemalloc = PETSC_FALSE;
2384: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2385: sbaij->free_a = PETSC_FALSE;
2386: sbaij->free_ij = PETSC_FALSE;
2387: sbaij->free_imax_ilen = PETSC_TRUE;
2389: for (ii=0; ii<m; ii++) {
2390: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2391: if (PetscUnlikelyDebug(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]);
2392: }
2393: if (PetscDefined(USE_DEBUG)) {
2394: for (ii=0; ii<sbaij->i[m]; ii++) {
2395: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2396: 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]);
2397: }
2398: }
2400: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2401: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2402: return(0);
2403: }
2405: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2406: {
2408: PetscMPIInt size;
2411: MPI_Comm_size(comm,&size);
2412: if (size == 1 && scall == MAT_REUSE_MATRIX) {
2413: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2414: } else {
2415: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);
2416: }
2417: return(0);
2418: }