Actual source code: sbaij.c
petsc-3.14.6 2021-03-30
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_NEW_DIAGONALS:
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 (idx) {PetscFree(*idx);}
313: if (v) {PetscFree(*v);}
314: return(0);
315: }
317: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
318: {
319: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
322: a->getrow_utriangular = PETSC_TRUE;
323: return(0);
324: }
326: PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
327: {
328: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
331: a->getrow_utriangular = PETSC_FALSE;
332: return(0);
333: }
335: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B)
336: {
340: if (reuse == MAT_INITIAL_MATRIX) {
341: MatDuplicate(A,MAT_COPY_VALUES,B);
342: } else if (reuse == MAT_REUSE_MATRIX) {
343: MatCopy(A,*B,SAME_NONZERO_PATTERN);
344: }
345: return(0);
346: }
348: PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
349: {
350: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
351: PetscErrorCode ierr;
352: PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
353: PetscViewerFormat format;
354: PetscInt *diag;
357: PetscViewerGetFormat(viewer,&format);
358: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
359: PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);
360: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
361: Mat aij;
362: const char *matname;
364: if (A->factortype && bs>1) {
365: PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");
366: return(0);
367: }
368: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);
369: PetscObjectGetName((PetscObject)A,&matname);
370: PetscObjectSetName((PetscObject)aij,matname);
371: MatView(aij,viewer);
372: MatDestroy(&aij);
373: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
374: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
375: for (i=0; i<a->mbs; i++) {
376: for (j=0; j<bs; j++) {
377: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
378: for (k=a->i[i]; k<a->i[i+1]; k++) {
379: for (l=0; l<bs; l++) {
380: #if defined(PETSC_USE_COMPLEX)
381: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
382: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
383: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
384: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
385: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
386: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
387: } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
388: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
389: }
390: #else
391: if (a->a[bs2*k + l*bs + j] != 0.0) {
392: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
393: }
394: #endif
395: }
396: }
397: PetscViewerASCIIPrintf(viewer,"\n");
398: }
399: }
400: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
401: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
402: return(0);
403: } else {
404: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
405: if (A->factortype) { /* for factored matrix */
406: if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet");
408: diag=a->diag;
409: for (i=0; i<a->mbs; i++) { /* for row block i */
410: PetscViewerASCIIPrintf(viewer,"row %D:",i);
411: /* diagonal entry */
412: #if defined(PETSC_USE_COMPLEX)
413: if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
414: 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]]));
415: } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
416: 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]]));
417: } else {
418: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));
419: }
420: #else
421: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));
422: #endif
423: /* off-diagonal entries */
424: for (k=a->i[i]; k<a->i[i+1]-1; k++) {
425: #if defined(PETSC_USE_COMPLEX)
426: if (PetscImaginaryPart(a->a[k]) > 0.0) {
427: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));
428: } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
429: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));
430: } else {
431: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));
432: }
433: #else
434: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);
435: #endif
436: }
437: PetscViewerASCIIPrintf(viewer,"\n");
438: }
440: } else { /* for non-factored matrix */
441: for (i=0; i<a->mbs; i++) { /* for row block i */
442: for (j=0; j<bs; j++) { /* for row bs*i + j */
443: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
444: for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */
445: for (l=0; l<bs; l++) { /* for column */
446: #if defined(PETSC_USE_COMPLEX)
447: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
448: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
449: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
450: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
451: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
452: (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
453: } else {
454: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
455: }
456: #else
457: PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
458: #endif
459: }
460: }
461: PetscViewerASCIIPrintf(viewer,"\n");
462: }
463: }
464: }
465: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
466: }
467: PetscViewerFlush(viewer);
468: return(0);
469: }
471: #include <petscdraw.h>
472: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
473: {
474: Mat A = (Mat) Aa;
475: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
477: PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
478: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
479: MatScalar *aa;
480: PetscViewer viewer;
483: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
484: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
486: /* loop over matrix elements drawing boxes */
488: PetscDrawCollectiveBegin(draw);
489: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
490: /* Blue for negative, Cyan for zero and Red for positive */
491: color = PETSC_DRAW_BLUE;
492: for (i=0,row=0; i<mbs; i++,row+=bs) {
493: for (j=a->i[i]; j<a->i[i+1]; j++) {
494: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
495: x_l = a->j[j]*bs; x_r = x_l + 1.0;
496: aa = a->a + j*bs2;
497: for (k=0; k<bs; k++) {
498: for (l=0; l<bs; l++) {
499: if (PetscRealPart(*aa++) >= 0.) continue;
500: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
501: }
502: }
503: }
504: }
505: color = PETSC_DRAW_CYAN;
506: for (i=0,row=0; i<mbs; i++,row+=bs) {
507: for (j=a->i[i]; j<a->i[i+1]; j++) {
508: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
509: x_l = a->j[j]*bs; x_r = x_l + 1.0;
510: aa = a->a + j*bs2;
511: for (k=0; k<bs; k++) {
512: for (l=0; l<bs; l++) {
513: if (PetscRealPart(*aa++) != 0.) continue;
514: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
515: }
516: }
517: }
518: }
519: color = PETSC_DRAW_RED;
520: for (i=0,row=0; i<mbs; i++,row+=bs) {
521: for (j=a->i[i]; j<a->i[i+1]; j++) {
522: y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
523: x_l = a->j[j]*bs; x_r = x_l + 1.0;
524: aa = a->a + j*bs2;
525: for (k=0; k<bs; k++) {
526: for (l=0; l<bs; l++) {
527: if (PetscRealPart(*aa++) <= 0.) continue;
528: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
529: }
530: }
531: }
532: }
533: PetscDrawCollectiveEnd(draw);
534: return(0);
535: }
537: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
538: {
540: PetscReal xl,yl,xr,yr,w,h;
541: PetscDraw draw;
542: PetscBool isnull;
545: PetscViewerDrawGetDraw(viewer,0,&draw);
546: PetscDrawIsNull(draw,&isnull);
547: if (isnull) return(0);
549: xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
550: xr += w; yr += h; xl = -w; yl = -h;
551: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
552: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
553: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
554: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
555: PetscDrawSave(draw);
556: return(0);
557: }
559: /* Used for both MPIBAIJ and MPISBAIJ matrices */
560: #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary
562: PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
563: {
565: PetscBool iascii,isbinary,isdraw;
568: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
569: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
570: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
571: if (iascii) {
572: MatView_SeqSBAIJ_ASCII(A,viewer);
573: } else if (isbinary) {
574: MatView_SeqSBAIJ_Binary(A,viewer);
575: } else if (isdraw) {
576: MatView_SeqSBAIJ_Draw(A,viewer);
577: } else {
578: Mat B;
579: const char *matname;
580: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
581: PetscObjectGetName((PetscObject)A,&matname);
582: PetscObjectSetName((PetscObject)B,matname);
583: MatView(B,viewer);
584: MatDestroy(&B);
585: }
586: return(0);
587: }
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: }
634: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
635: {
636: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
637: PetscErrorCode ierr;
638: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
639: PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
640: PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
641: PetscBool roworiented=a->roworiented;
642: const PetscScalar *value = v;
643: MatScalar *ap,*aa = a->a,*bap;
646: if (roworiented) stepval = (n-1)*bs;
647: else stepval = (m-1)*bs;
649: for (k=0; k<m; k++) { /* loop over added rows */
650: row = im[k];
651: if (row < 0) continue;
652: 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);
653: rp = aj + ai[row];
654: ap = aa + bs2*ai[row];
655: rmax = imax[row];
656: nrow = ailen[row];
657: low = 0;
658: high = nrow;
659: for (l=0; l<n; l++) { /* loop over added columns */
660: if (in[l] < 0) continue;
661: col = in[l];
662: 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);
663: if (col < row) {
664: if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
665: 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)");
666: }
667: if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
668: else value = v + l*(stepval+bs)*bs + k*bs;
670: if (col <= lastcol) low = 0;
671: else high = nrow;
673: lastcol = col;
674: while (high-low > 7) {
675: t = (low+high)/2;
676: if (rp[t] > col) high = t;
677: else low = t;
678: }
679: for (i=low; i<high; i++) {
680: if (rp[i] > col) break;
681: if (rp[i] == col) {
682: bap = ap + bs2*i;
683: if (roworiented) {
684: if (is == ADD_VALUES) {
685: for (ii=0; ii<bs; ii++,value+=stepval) {
686: for (jj=ii; jj<bs2; jj+=bs) {
687: bap[jj] += *value++;
688: }
689: }
690: } else {
691: for (ii=0; ii<bs; ii++,value+=stepval) {
692: for (jj=ii; jj<bs2; jj+=bs) {
693: bap[jj] = *value++;
694: }
695: }
696: }
697: } else {
698: if (is == ADD_VALUES) {
699: for (ii=0; ii<bs; ii++,value+=stepval) {
700: for (jj=0; jj<bs; jj++) {
701: *bap++ += *value++;
702: }
703: }
704: } else {
705: for (ii=0; ii<bs; ii++,value+=stepval) {
706: for (jj=0; jj<bs; jj++) {
707: *bap++ = *value++;
708: }
709: }
710: }
711: }
712: goto noinsert2;
713: }
714: }
715: if (nonew == 1) goto noinsert2;
716: 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);
717: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
718: N = nrow++ - 1; high++;
719: /* shift up all the later entries in this row */
720: PetscArraymove(rp+i+1,rp+i,N-i+1);
721: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
722: PetscArrayzero(ap+bs2*i,bs2);
723: rp[i] = col;
724: bap = ap + bs2*i;
725: if (roworiented) {
726: for (ii=0; ii<bs; ii++,value+=stepval) {
727: for (jj=ii; jj<bs2; jj+=bs) {
728: bap[jj] = *value++;
729: }
730: }
731: } else {
732: for (ii=0; ii<bs; ii++,value+=stepval) {
733: for (jj=0; jj<bs; jj++) {
734: *bap++ = *value++;
735: }
736: }
737: }
738: noinsert2:;
739: low = i;
740: }
741: ailen[row] = nrow;
742: }
743: return(0);
744: }
746: /*
747: This is not yet used
748: */
749: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
750: {
751: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
753: const PetscInt *ai = a->i, *aj = a->j,*cols;
754: PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
755: PetscBool flag;
758: PetscMalloc1(m,&ns);
759: while (i < m) {
760: nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
761: /* Limits the number of elements in a node to 'a->inode.limit' */
762: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
763: nzy = ai[j+1] - ai[j];
764: if (nzy != (nzx - j + i)) break;
765: PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);
766: if (!flag) break;
767: }
768: ns[node_count++] = blk_size;
770: i = j;
771: }
772: if (!a->inode.size && m && node_count > .9*m) {
773: PetscFree(ns);
774: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
775: } else {
776: a->inode.node_count = node_count;
778: PetscMalloc1(node_count,&a->inode.size);
779: PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));
780: PetscArraycpy(a->inode.size,ns,node_count);
781: PetscFree(ns);
782: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
784: /* count collections of adjacent columns in each inode */
785: row = 0;
786: cnt = 0;
787: for (i=0; i<node_count; i++) {
788: cols = aj + ai[row] + a->inode.size[i];
789: nz = ai[row+1] - ai[row] - a->inode.size[i];
790: for (j=1; j<nz; j++) {
791: if (cols[j] != cols[j-1]+1) cnt++;
792: }
793: cnt++;
794: row += a->inode.size[i];
795: }
796: PetscMalloc1(2*cnt,&counts);
797: cnt = 0;
798: row = 0;
799: for (i=0; i<node_count; i++) {
800: cols = aj + ai[row] + a->inode.size[i];
801: counts[2*cnt] = cols[0];
802: nz = ai[row+1] - ai[row] - a->inode.size[i];
803: cnt2 = 1;
804: for (j=1; j<nz; j++) {
805: if (cols[j] != cols[j-1]+1) {
806: counts[2*(cnt++)+1] = cnt2;
807: counts[2*cnt] = cols[j];
808: cnt2 = 1;
809: } else cnt2++;
810: }
811: counts[2*(cnt++)+1] = cnt2;
812: row += a->inode.size[i];
813: }
814: PetscIntView(2*cnt,counts,NULL);
815: }
816: return(0);
817: }
819: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
820: {
821: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
823: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
824: PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
825: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
826: MatScalar *aa = a->a,*ap;
829: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
831: if (m) rmax = ailen[0];
832: for (i=1; i<mbs; i++) {
833: /* move each row back by the amount of empty slots (fshift) before it*/
834: fshift += imax[i-1] - ailen[i-1];
835: rmax = PetscMax(rmax,ailen[i]);
836: if (fshift) {
837: ip = aj + ai[i];
838: ap = aa + bs2*ai[i];
839: N = ailen[i];
840: PetscArraymove(ip-fshift,ip,N);
841: PetscArraymove(ap-bs2*fshift,ap,bs2*N);
842: }
843: ai[i] = ai[i-1] + ailen[i-1];
844: }
845: if (mbs) {
846: fshift += imax[mbs-1] - ailen[mbs-1];
847: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
848: }
849: /* reset ilen and imax for each row */
850: for (i=0; i<mbs; i++) {
851: ailen[i] = imax[i] = ai[i+1] - ai[i];
852: }
853: a->nz = ai[mbs];
855: /* diagonals may have moved, reset it */
856: if (a->diag) {
857: PetscArraycpy(a->diag,ai,mbs);
858: }
859: 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);
861: 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);
862: PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
863: PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);
865: A->info.mallocs += a->reallocs;
866: a->reallocs = 0;
867: A->info.nz_unneeded = (PetscReal)fshift*bs2;
868: a->idiagvalid = PETSC_FALSE;
869: a->rmax = rmax;
871: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
872: if (a->jshort && a->free_jshort) {
873: /* when matrix data structure is changed, previous jshort must be replaced */
874: PetscFree(a->jshort);
875: }
876: PetscMalloc1(a->i[A->rmap->n],&a->jshort);
877: PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));
878: for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
879: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
880: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
881: a->free_jshort = PETSC_TRUE;
882: }
883: return(0);
884: }
886: /*
887: This function returns an array of flags which indicate the locations of contiguous
888: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
889: then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
890: Assume: sizes should be long enough to hold all the values.
891: */
892: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
893: {
894: PetscInt i,j,k,row;
895: PetscBool flg;
898: for (i=0,j=0; i<n; j++) {
899: row = idx[i];
900: if (row%bs!=0) { /* Not the begining of a block */
901: sizes[j] = 1;
902: i++;
903: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
904: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
905: i++;
906: } else { /* Begining of the block, so check if the complete block exists */
907: flg = PETSC_TRUE;
908: for (k=1; k<bs; k++) {
909: if (row+k != idx[i+k]) { /* break in the block */
910: flg = PETSC_FALSE;
911: break;
912: }
913: }
914: if (flg) { /* No break in the bs */
915: sizes[j] = bs;
916: i += bs;
917: } else {
918: sizes[j] = 1;
919: i++;
920: }
921: }
922: }
923: *bs_max = j;
924: return(0);
925: }
928: /* Only add/insert a(i,j) with i<=j (blocks).
929: Any a(i,j) with i>j input by user is ingored.
930: */
932: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
933: {
934: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
936: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
937: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
938: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
939: PetscInt ridx,cidx,bs2=a->bs2;
940: MatScalar *ap,value,*aa=a->a,*bap;
943: for (k=0; k<m; k++) { /* loop over added rows */
944: row = im[k]; /* row number */
945: brow = row/bs; /* block row number */
946: if (row < 0) continue;
947: 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);
948: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
949: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
950: rmax = imax[brow]; /* maximum space allocated for this row */
951: nrow = ailen[brow]; /* actual length of this row */
952: low = 0;
953: high = nrow;
954: for (l=0; l<n; l++) { /* loop over added columns */
955: if (in[l] < 0) continue;
956: 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);
957: col = in[l];
958: bcol = col/bs; /* block col number */
960: if (brow > bcol) {
961: if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
962: 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)");
963: }
965: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
966: if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
967: /* element value a(k,l) */
968: if (roworiented) value = v[l + k*n];
969: else value = v[k + l*m];
971: /* move pointer bap to a(k,l) quickly and add/insert value */
972: if (col <= lastcol) low = 0;
973: else high = nrow;
975: lastcol = col;
976: while (high-low > 7) {
977: t = (low+high)/2;
978: if (rp[t] > bcol) high = t;
979: else low = t;
980: }
981: for (i=low; i<high; i++) {
982: if (rp[i] > bcol) break;
983: if (rp[i] == bcol) {
984: bap = ap + bs2*i + bs*cidx + ridx;
985: if (is == ADD_VALUES) *bap += value;
986: else *bap = value;
987: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
988: if (brow == bcol && ridx < cidx) {
989: bap = ap + bs2*i + bs*ridx + cidx;
990: if (is == ADD_VALUES) *bap += value;
991: else *bap = value;
992: }
993: goto noinsert1;
994: }
995: }
997: if (nonew == 1) goto noinsert1;
998: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
999: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
1001: N = nrow++ - 1; high++;
1002: /* shift up all the later entries in this row */
1003: PetscArraymove(rp+i+1,rp+i,N-i+1);
1004: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
1005: PetscArrayzero(ap+bs2*i,bs2);
1006: rp[i] = bcol;
1007: ap[bs2*i + bs*cidx + ridx] = value;
1008: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1009: if (brow == bcol && ridx < cidx) {
1010: ap[bs2*i + bs*ridx + cidx] = value;
1011: }
1012: A->nonzerostate++;
1013: noinsert1:;
1014: low = i;
1015: }
1016: } /* end of loop over added columns */
1017: ailen[brow] = nrow;
1018: } /* end of loop over added rows */
1019: return(0);
1020: }
1022: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
1023: {
1024: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1025: Mat outA;
1027: PetscBool row_identity;
1030: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
1031: ISIdentity(row,&row_identity);
1032: if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1033: 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()! */
1035: outA = inA;
1036: inA->factortype = MAT_FACTOR_ICC;
1037: PetscFree(inA->solvertype);
1038: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
1040: MatMarkDiagonal_SeqSBAIJ(inA);
1041: MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);
1043: PetscObjectReference((PetscObject)row);
1044: ISDestroy(&a->row);
1045: a->row = row;
1046: PetscObjectReference((PetscObject)row);
1047: ISDestroy(&a->col);
1048: a->col = row;
1050: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1051: if (a->icol) {ISInvertPermutation(row,PETSC_DECIDE, &a->icol);}
1052: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
1054: if (!a->solve_work) {
1055: PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
1056: PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
1057: }
1059: MatCholeskyFactorNumeric(outA,inA,info);
1060: return(0);
1061: }
1063: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1064: {
1065: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1066: PetscInt i,nz,n;
1070: nz = baij->maxnz;
1071: n = mat->cmap->n;
1072: for (i=0; i<nz; i++) baij->j[i] = indices[i];
1074: baij->nz = nz;
1075: for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1077: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1078: return(0);
1079: }
1081: /*@
1082: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1083: in the matrix.
1085: Input Parameters:
1086: + mat - the SeqSBAIJ matrix
1087: - indices - the column indices
1089: Level: advanced
1091: Notes:
1092: This can be called if you have precomputed the nonzero structure of the
1093: matrix and want to provide it to the matrix object to improve the performance
1094: of the MatSetValues() operation.
1096: You MUST have set the correct numbers of nonzeros per row in the call to
1097: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1099: MUST be called before any calls to MatSetValues()
1101: .seealso: MatCreateSeqSBAIJ
1102: @*/
1103: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1104: {
1110: PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
1111: return(0);
1112: }
1114: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1115: {
1117: PetscBool isbaij;
1120: PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");
1121: if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name);
1122: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1123: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1124: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1125: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1127: 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");
1128: if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different");
1129: if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size");
1130: PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);
1131: PetscObjectStateIncrease((PetscObject)B);
1132: } else {
1133: MatGetRowUpperTriangular(A);
1134: MatCopy_Basic(A,B,str);
1135: MatRestoreRowUpperTriangular(A);
1136: }
1137: return(0);
1138: }
1140: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1141: {
1145: MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);
1146: return(0);
1147: }
1149: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1150: {
1151: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1154: *array = a->a;
1155: return(0);
1156: }
1158: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1159: {
1161: *array = NULL;
1162: return(0);
1163: }
1165: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1166: {
1167: PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1168: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1169: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1173: /* Set the number of nonzeros in the new matrix */
1174: MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
1175: return(0);
1176: }
1178: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1179: {
1180: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1182: PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1183: PetscBLASInt one = 1;
1186: if (str == SAME_NONZERO_PATTERN) {
1187: PetscScalar alpha = a;
1188: PetscBLASInt bnz;
1189: PetscBLASIntCast(x->nz*bs2,&bnz);
1190: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1191: PetscObjectStateIncrease((PetscObject)Y);
1192: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1193: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);
1194: MatAXPY_Basic(Y,a,X,str);
1195: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);
1196: } else {
1197: Mat B;
1198: PetscInt *nnz;
1199: if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
1200: MatGetRowUpperTriangular(X);
1201: MatGetRowUpperTriangular(Y);
1202: PetscMalloc1(Y->rmap->N,&nnz);
1203: MatCreate(PetscObjectComm((PetscObject)Y),&B);
1204: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
1205: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
1206: MatSetBlockSizesFromMats(B,Y,Y);
1207: MatSetType(B,((PetscObject)Y)->type_name);
1208: MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);
1209: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1211: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1213: MatHeaderReplace(Y,&B);
1214: PetscFree(nnz);
1215: MatRestoreRowUpperTriangular(X);
1216: MatRestoreRowUpperTriangular(Y);
1217: }
1218: return(0);
1219: }
1221: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1222: {
1224: *flg = PETSC_TRUE;
1225: return(0);
1226: }
1228: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1229: {
1231: *flg = PETSC_TRUE;
1232: return(0);
1233: }
1235: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1236: {
1238: *flg = PETSC_FALSE;
1239: return(0);
1240: }
1242: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1243: {
1244: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1245: PetscInt i,nz = a->bs2*a->i[a->mbs];
1246: MatScalar *aa = a->a;
1249: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1250: return(0);
1251: }
1253: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1254: {
1255: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1256: PetscInt i,nz = a->bs2*a->i[a->mbs];
1257: MatScalar *aa = a->a;
1260: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1261: return(0);
1262: }
1264: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1265: {
1266: Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1267: PetscErrorCode ierr;
1268: PetscInt i,j,k,count;
1269: PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1270: PetscScalar zero = 0.0;
1271: MatScalar *aa;
1272: const PetscScalar *xx;
1273: PetscScalar *bb;
1274: PetscBool *zeroed,vecs = PETSC_FALSE;
1277: /* fix right hand side if needed */
1278: if (x && b) {
1279: VecGetArrayRead(x,&xx);
1280: VecGetArray(b,&bb);
1281: vecs = PETSC_TRUE;
1282: }
1284: /* zero the columns */
1285: PetscCalloc1(A->rmap->n,&zeroed);
1286: for (i=0; i<is_n; i++) {
1287: 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]);
1288: zeroed[is_idx[i]] = PETSC_TRUE;
1289: }
1290: if (vecs) {
1291: for (i=0; i<A->rmap->N; i++) {
1292: row = i/bs;
1293: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1294: for (k=0; k<bs; k++) {
1295: col = bs*baij->j[j] + k;
1296: if (col <= i) continue;
1297: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1298: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1299: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1300: }
1301: }
1302: }
1303: for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1304: }
1306: for (i=0; i<A->rmap->N; i++) {
1307: if (!zeroed[i]) {
1308: row = i/bs;
1309: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1310: for (k=0; k<bs; k++) {
1311: col = bs*baij->j[j] + k;
1312: if (zeroed[col]) {
1313: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1314: aa[0] = 0.0;
1315: }
1316: }
1317: }
1318: }
1319: }
1320: PetscFree(zeroed);
1321: if (vecs) {
1322: VecRestoreArrayRead(x,&xx);
1323: VecRestoreArray(b,&bb);
1324: }
1326: /* zero the rows */
1327: for (i=0; i<is_n; i++) {
1328: row = is_idx[i];
1329: count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1330: aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1331: for (k=0; k<count; k++) {
1332: aa[0] = zero;
1333: aa += bs;
1334: }
1335: if (diag != 0.0) {
1336: (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
1337: }
1338: }
1339: MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);
1340: return(0);
1341: }
1343: PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1344: {
1346: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1349: if (!Y->preallocated || !aij->nz) {
1350: MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
1351: }
1352: MatShift_Basic(Y,a);
1353: return(0);
1354: }
1356: /* -------------------------------------------------------------------*/
1357: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1358: MatGetRow_SeqSBAIJ,
1359: MatRestoreRow_SeqSBAIJ,
1360: MatMult_SeqSBAIJ_N,
1361: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1362: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1363: MatMultAdd_SeqSBAIJ_N,
1364: NULL,
1365: NULL,
1366: NULL,
1367: /* 10*/ NULL,
1368: NULL,
1369: MatCholeskyFactor_SeqSBAIJ,
1370: MatSOR_SeqSBAIJ,
1371: MatTranspose_SeqSBAIJ,
1372: /* 15*/ MatGetInfo_SeqSBAIJ,
1373: MatEqual_SeqSBAIJ,
1374: MatGetDiagonal_SeqSBAIJ,
1375: MatDiagonalScale_SeqSBAIJ,
1376: MatNorm_SeqSBAIJ,
1377: /* 20*/ NULL,
1378: MatAssemblyEnd_SeqSBAIJ,
1379: MatSetOption_SeqSBAIJ,
1380: MatZeroEntries_SeqSBAIJ,
1381: /* 24*/ NULL,
1382: NULL,
1383: NULL,
1384: NULL,
1385: NULL,
1386: /* 29*/ MatSetUp_SeqSBAIJ,
1387: NULL,
1388: NULL,
1389: NULL,
1390: NULL,
1391: /* 34*/ MatDuplicate_SeqSBAIJ,
1392: NULL,
1393: NULL,
1394: NULL,
1395: MatICCFactor_SeqSBAIJ,
1396: /* 39*/ MatAXPY_SeqSBAIJ,
1397: MatCreateSubMatrices_SeqSBAIJ,
1398: MatIncreaseOverlap_SeqSBAIJ,
1399: MatGetValues_SeqSBAIJ,
1400: MatCopy_SeqSBAIJ,
1401: /* 44*/ NULL,
1402: MatScale_SeqSBAIJ,
1403: MatShift_SeqSBAIJ,
1404: NULL,
1405: MatZeroRowsColumns_SeqSBAIJ,
1406: /* 49*/ NULL,
1407: MatGetRowIJ_SeqSBAIJ,
1408: MatRestoreRowIJ_SeqSBAIJ,
1409: NULL,
1410: NULL,
1411: /* 54*/ NULL,
1412: NULL,
1413: NULL,
1414: NULL,
1415: MatSetValuesBlocked_SeqSBAIJ,
1416: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1417: NULL,
1418: NULL,
1419: NULL,
1420: NULL,
1421: /* 64*/ NULL,
1422: NULL,
1423: NULL,
1424: NULL,
1425: NULL,
1426: /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1427: NULL,
1428: MatConvert_MPISBAIJ_Basic,
1429: NULL,
1430: NULL,
1431: /* 74*/ NULL,
1432: NULL,
1433: NULL,
1434: NULL,
1435: NULL,
1436: /* 79*/ NULL,
1437: NULL,
1438: NULL,
1439: MatGetInertia_SeqSBAIJ,
1440: MatLoad_SeqSBAIJ,
1441: /* 84*/ MatIsSymmetric_SeqSBAIJ,
1442: MatIsHermitian_SeqSBAIJ,
1443: MatIsStructurallySymmetric_SeqSBAIJ,
1444: NULL,
1445: NULL,
1446: /* 89*/ NULL,
1447: NULL,
1448: NULL,
1449: NULL,
1450: NULL,
1451: /* 94*/ NULL,
1452: NULL,
1453: NULL,
1454: NULL,
1455: NULL,
1456: /* 99*/ NULL,
1457: NULL,
1458: NULL,
1459: NULL,
1460: NULL,
1461: /*104*/ NULL,
1462: MatRealPart_SeqSBAIJ,
1463: MatImaginaryPart_SeqSBAIJ,
1464: MatGetRowUpperTriangular_SeqSBAIJ,
1465: MatRestoreRowUpperTriangular_SeqSBAIJ,
1466: /*109*/ NULL,
1467: NULL,
1468: NULL,
1469: NULL,
1470: MatMissingDiagonal_SeqSBAIJ,
1471: /*114*/ NULL,
1472: NULL,
1473: NULL,
1474: NULL,
1475: NULL,
1476: /*119*/ NULL,
1477: NULL,
1478: NULL,
1479: NULL,
1480: NULL,
1481: /*124*/ NULL,
1482: NULL,
1483: NULL,
1484: NULL,
1485: NULL,
1486: /*129*/ NULL,
1487: NULL,
1488: NULL,
1489: NULL,
1490: NULL,
1491: /*134*/ NULL,
1492: NULL,
1493: NULL,
1494: NULL,
1495: NULL,
1496: /*139*/ MatSetBlockSizes_Default,
1497: NULL,
1498: NULL,
1499: NULL,
1500: NULL,
1501: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ
1502: };
1504: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1505: {
1506: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1507: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1511: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1513: /* allocate space for values if not already there */
1514: if (!aij->saved_values) {
1515: PetscMalloc1(nz+1,&aij->saved_values);
1516: }
1518: /* copy values over */
1519: PetscArraycpy(aij->saved_values,aij->a,nz);
1520: return(0);
1521: }
1523: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1524: {
1525: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1527: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1530: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1531: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1533: /* copy values over */
1534: PetscArraycpy(aij->a,aij->saved_values,nz);
1535: return(0);
1536: }
1538: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1539: {
1540: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1542: PetscInt i,mbs,nbs,bs2;
1543: PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1546: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1548: MatSetBlockSize(B,PetscAbs(bs));
1549: PetscLayoutSetUp(B->rmap);
1550: PetscLayoutSetUp(B->cmap);
1551: 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);
1552: PetscLayoutGetBlockSize(B->rmap,&bs);
1554: B->preallocated = PETSC_TRUE;
1556: mbs = B->rmap->N/bs;
1557: nbs = B->cmap->n/bs;
1558: bs2 = bs*bs;
1560: 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");
1562: if (nz == MAT_SKIP_ALLOCATION) {
1563: skipallocation = PETSC_TRUE;
1564: nz = 0;
1565: }
1567: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1568: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1569: if (nnz) {
1570: for (i=0; i<mbs; i++) {
1571: 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]);
1572: 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);
1573: }
1574: }
1576: B->ops->mult = MatMult_SeqSBAIJ_N;
1577: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1578: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1579: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1581: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1582: if (!flg) {
1583: switch (bs) {
1584: case 1:
1585: B->ops->mult = MatMult_SeqSBAIJ_1;
1586: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1587: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1588: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1589: break;
1590: case 2:
1591: B->ops->mult = MatMult_SeqSBAIJ_2;
1592: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1593: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1594: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1595: break;
1596: case 3:
1597: B->ops->mult = MatMult_SeqSBAIJ_3;
1598: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1599: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1600: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1601: break;
1602: case 4:
1603: B->ops->mult = MatMult_SeqSBAIJ_4;
1604: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1605: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1606: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1607: break;
1608: case 5:
1609: B->ops->mult = MatMult_SeqSBAIJ_5;
1610: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1611: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1612: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1613: break;
1614: case 6:
1615: B->ops->mult = MatMult_SeqSBAIJ_6;
1616: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1617: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1618: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1619: break;
1620: case 7:
1621: B->ops->mult = MatMult_SeqSBAIJ_7;
1622: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1623: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1624: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1625: break;
1626: }
1627: }
1629: b->mbs = mbs;
1630: b->nbs = nbs;
1631: if (!skipallocation) {
1632: if (!b->imax) {
1633: PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
1635: b->free_imax_ilen = PETSC_TRUE;
1637: PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));
1638: }
1639: if (!nnz) {
1640: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1641: else if (nz <= 0) nz = 1;
1642: nz = PetscMin(nbs,nz);
1643: for (i=0; i<mbs; i++) b->imax[i] = nz;
1644: PetscIntMultError(nz,mbs,&nz);
1645: } else {
1646: PetscInt64 nz64 = 0;
1647: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1648: PetscIntCast(nz64,&nz);
1649: }
1650: /* b->ilen will count nonzeros in each block row so far. */
1651: for (i=0; i<mbs; i++) b->ilen[i] = 0;
1652: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1654: /* allocate the matrix space */
1655: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
1656: PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
1657: PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
1658: PetscArrayzero(b->a,nz*bs2);
1659: PetscArrayzero(b->j,nz);
1661: b->singlemalloc = PETSC_TRUE;
1663: /* pointer to beginning of each row */
1664: b->i[0] = 0;
1665: for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1667: b->free_a = PETSC_TRUE;
1668: b->free_ij = PETSC_TRUE;
1669: } else {
1670: b->free_a = PETSC_FALSE;
1671: b->free_ij = PETSC_FALSE;
1672: }
1674: b->bs2 = bs2;
1675: b->nz = 0;
1676: b->maxnz = nz;
1677: b->inew = NULL;
1678: b->jnew = NULL;
1679: b->anew = NULL;
1680: b->a2anew = NULL;
1681: b->permute = PETSC_FALSE;
1683: B->was_assembled = PETSC_FALSE;
1684: B->assembled = PETSC_FALSE;
1685: if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
1686: return(0);
1687: }
1689: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1690: {
1691: PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1692: PetscScalar *values=NULL;
1693: PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1697: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
1698: PetscLayoutSetBlockSize(B->rmap,bs);
1699: PetscLayoutSetBlockSize(B->cmap,bs);
1700: PetscLayoutSetUp(B->rmap);
1701: PetscLayoutSetUp(B->cmap);
1702: PetscLayoutGetBlockSize(B->rmap,&bs);
1703: m = B->rmap->n/bs;
1705: if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]);
1706: PetscMalloc1(m+1,&nnz);
1707: for (i=0; i<m; i++) {
1708: nz = ii[i+1] - ii[i];
1709: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz);
1710: anz = 0;
1711: for (j=0; j<nz; j++) {
1712: /* count only values on the diagonal or above */
1713: if (jj[ii[i] + j] >= i) {
1714: anz = nz - j;
1715: break;
1716: }
1717: }
1718: nz_max = PetscMax(nz_max,anz);
1719: nnz[i] = anz;
1720: }
1721: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1722: PetscFree(nnz);
1724: values = (PetscScalar*)V;
1725: if (!values) {
1726: PetscCalloc1(bs*bs*nz_max,&values);
1727: }
1728: for (i=0; i<m; i++) {
1729: PetscInt ncols = ii[i+1] - ii[i];
1730: const PetscInt *icols = jj + ii[i];
1731: if (!roworiented || bs == 1) {
1732: const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1733: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
1734: } else {
1735: for (j=0; j<ncols; j++) {
1736: const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1737: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
1738: }
1739: }
1740: }
1741: if (!V) { PetscFree(values); }
1742: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1743: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1744: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1745: return(0);
1746: }
1748: /*
1749: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1750: */
1751: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1752: {
1754: PetscBool flg = PETSC_FALSE;
1755: PetscInt bs = B->rmap->bs;
1758: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1759: if (flg) bs = 8;
1761: if (!natural) {
1762: switch (bs) {
1763: case 1:
1764: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1765: break;
1766: case 2:
1767: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1768: break;
1769: case 3:
1770: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1771: break;
1772: case 4:
1773: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1774: break;
1775: case 5:
1776: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1777: break;
1778: case 6:
1779: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1780: break;
1781: case 7:
1782: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1783: break;
1784: default:
1785: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1786: break;
1787: }
1788: } else {
1789: switch (bs) {
1790: case 1:
1791: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1792: break;
1793: case 2:
1794: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1795: break;
1796: case 3:
1797: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1798: break;
1799: case 4:
1800: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1801: break;
1802: case 5:
1803: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1804: break;
1805: case 6:
1806: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1807: break;
1808: case 7:
1809: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1810: break;
1811: default:
1812: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1813: break;
1814: }
1815: }
1816: return(0);
1817: }
1819: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
1820: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*);
1822: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1823: {
1824: PetscInt n = A->rmap->n;
1828: #if defined(PETSC_USE_COMPLEX)
1829: 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");
1830: #endif
1832: MatCreate(PetscObjectComm((PetscObject)A),B);
1833: MatSetSizes(*B,n,n,n,n);
1834: if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1835: MatSetType(*B,MATSEQSBAIJ);
1836: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
1838: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1839: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1840: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1842: (*B)->factortype = ftype;
1843: (*B)->useordering = PETSC_TRUE;
1844: PetscFree((*B)->solvertype);
1845: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
1846: return(0);
1847: }
1849: /*@C
1850: MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1852: Not Collective
1854: Input Parameter:
1855: . mat - a MATSEQSBAIJ matrix
1857: Output Parameter:
1858: . array - pointer to the data
1860: Level: intermediate
1862: .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1863: @*/
1864: PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1865: {
1869: PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
1870: return(0);
1871: }
1873: /*@C
1874: MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1876: Not Collective
1878: Input Parameters:
1879: + mat - a MATSEQSBAIJ matrix
1880: - array - pointer to the data
1882: Level: intermediate
1884: .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1885: @*/
1886: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1887: {
1891: PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
1892: return(0);
1893: }
1895: /*MC
1896: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1897: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1899: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1900: can call MatSetOption(Mat, MAT_HERMITIAN).
1902: Options Database Keys:
1903: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1905: Notes:
1906: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1907: stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1908: 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.
1910: The number of rows in the matrix must be less than or equal to the number of columns
1912: Level: beginner
1914: .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1915: M*/
1916: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1917: {
1918: Mat_SeqSBAIJ *b;
1920: PetscMPIInt size;
1921: PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1924: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1925: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1927: PetscNewLog(B,&b);
1928: B->data = (void*)b;
1929: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1931: B->ops->destroy = MatDestroy_SeqSBAIJ;
1932: B->ops->view = MatView_SeqSBAIJ;
1933: b->row = NULL;
1934: b->icol = NULL;
1935: b->reallocs = 0;
1936: b->saved_values = NULL;
1937: b->inode.limit = 5;
1938: b->inode.max_limit = 5;
1940: b->roworiented = PETSC_TRUE;
1941: b->nonew = 0;
1942: b->diag = NULL;
1943: b->solve_work = NULL;
1944: b->mult_work = NULL;
1945: B->spptr = NULL;
1946: B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1947: b->keepnonzeropattern = PETSC_FALSE;
1949: b->inew = NULL;
1950: b->jnew = NULL;
1951: b->anew = NULL;
1952: b->a2anew = NULL;
1953: b->permute = PETSC_FALSE;
1955: b->ignore_ltriangular = PETSC_TRUE;
1957: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);
1959: b->getrow_utriangular = PETSC_FALSE;
1961: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);
1963: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);
1964: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);
1965: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);
1966: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);
1967: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1968: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);
1969: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);
1970: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);
1971: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
1972: #if defined(PETSC_HAVE_ELEMENTAL)
1973: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);
1974: #endif
1975: #if defined(PETSC_HAVE_SCALAPACK)
1976: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);
1977: #endif
1979: B->symmetric = PETSC_TRUE;
1980: B->structurally_symmetric = PETSC_TRUE;
1981: B->symmetric_set = PETSC_TRUE;
1982: B->structurally_symmetric_set = PETSC_TRUE;
1983: B->symmetric_eternal = PETSC_TRUE;
1984: #if defined(PETSC_USE_COMPLEX)
1985: B->hermitian = PETSC_FALSE;
1986: B->hermitian_set = PETSC_FALSE;
1987: #else
1988: B->hermitian = PETSC_TRUE;
1989: B->hermitian_set = PETSC_TRUE;
1990: #endif
1992: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
1994: PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");
1995: PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);
1996: if (no_unroll) {
1997: PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");
1998: }
1999: PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);
2000: if (no_inode) {
2001: PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");
2002: }
2003: PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);
2004: PetscOptionsEnd();
2005: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
2006: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
2007: return(0);
2008: }
2010: /*@C
2011: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
2012: compressed row) format. For good matrix assembly performance the
2013: user should preallocate the matrix storage by setting the parameter nz
2014: (or the array nnz). By setting these parameters accurately, performance
2015: during matrix assembly can be increased by more than a factor of 50.
2017: Collective on Mat
2019: Input Parameters:
2020: + B - the symmetric matrix
2021: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2022: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2023: . nz - number of block nonzeros per block row (same for all rows)
2024: - nnz - array containing the number of block nonzeros in the upper triangular plus
2025: diagonal portion of each block (possibly different for each block row) or NULL
2027: Options Database Keys:
2028: + -mat_no_unroll - uses code that does not unroll the loops in the
2029: block calculations (much slower)
2030: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2032: Level: intermediate
2034: Notes:
2035: Specify the preallocated storage with either nz or nnz (not both).
2036: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2037: allocation. See Users-Manual: ch_mat for details.
2039: You can call MatGetInfo() to get information on how effective the preallocation was;
2040: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2041: You can also run with the option -info and look for messages with the string
2042: malloc in them to see if additional memory allocation was needed.
2044: If the nnz parameter is given then the nz parameter is ignored
2047: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2048: @*/
2049: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2050: {
2057: PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
2058: return(0);
2059: }
2061: /*@C
2062: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2064: Input Parameters:
2065: + B - the matrix
2066: . bs - size of block, the blocks are ALWAYS square.
2067: . i - the indices into j for the start of each local row (starts with zero)
2068: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2069: - v - optional values in the matrix
2071: Level: advanced
2073: Notes:
2074: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2075: may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2076: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2077: MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2078: block column and the second index is over columns within a block.
2080: Any entries below the diagonal are ignored
2082: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2083: and usually the numerical values as well
2085: .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2086: @*/
2087: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2088: {
2095: PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
2096: return(0);
2097: }
2099: /*@C
2100: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2101: compressed row) format. For good matrix assembly performance the
2102: user should preallocate the matrix storage by setting the parameter nz
2103: (or the array nnz). By setting these parameters accurately, performance
2104: during matrix assembly can be increased by more than a factor of 50.
2106: Collective
2108: Input Parameters:
2109: + comm - MPI communicator, set to PETSC_COMM_SELF
2110: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2111: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2112: . m - number of rows, or number of columns
2113: . nz - number of block nonzeros per block row (same for all rows)
2114: - nnz - array containing the number of block nonzeros in the upper triangular plus
2115: diagonal portion of each block (possibly different for each block row) or NULL
2117: Output Parameter:
2118: . A - the symmetric matrix
2120: Options Database Keys:
2121: + -mat_no_unroll - uses code that does not unroll the loops in the
2122: block calculations (much slower)
2123: - -mat_block_size - size of the blocks to use
2125: Level: intermediate
2127: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2128: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2129: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2131: Notes:
2132: The number of rows and columns must be divisible by blocksize.
2133: This matrix type does not support complex Hermitian operation.
2135: Specify the preallocated storage with either nz or nnz (not both).
2136: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2137: allocation. See Users-Manual: ch_mat for details.
2139: If the nnz parameter is given then the nz parameter is ignored
2141: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2142: @*/
2143: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2144: {
2148: MatCreate(comm,A);
2149: MatSetSizes(*A,m,n,m,n);
2150: MatSetType(*A,MATSEQSBAIJ);
2151: MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
2152: return(0);
2153: }
2155: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2156: {
2157: Mat C;
2158: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2160: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2163: if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");
2165: *B = NULL;
2166: MatCreate(PetscObjectComm((PetscObject)A),&C);
2167: MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
2168: MatSetBlockSizesFromMats(C,A,A);
2169: MatSetType(C,MATSEQSBAIJ);
2170: c = (Mat_SeqSBAIJ*)C->data;
2172: C->preallocated = PETSC_TRUE;
2173: C->factortype = A->factortype;
2174: c->row = NULL;
2175: c->icol = NULL;
2176: c->saved_values = NULL;
2177: c->keepnonzeropattern = a->keepnonzeropattern;
2178: C->assembled = PETSC_TRUE;
2180: PetscLayoutReference(A->rmap,&C->rmap);
2181: PetscLayoutReference(A->cmap,&C->cmap);
2182: c->bs2 = a->bs2;
2183: c->mbs = a->mbs;
2184: c->nbs = a->nbs;
2186: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2187: c->imax = a->imax;
2188: c->ilen = a->ilen;
2189: c->free_imax_ilen = PETSC_FALSE;
2190: } else {
2191: PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);
2192: PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));
2193: for (i=0; i<mbs; i++) {
2194: c->imax[i] = a->imax[i];
2195: c->ilen[i] = a->ilen[i];
2196: }
2197: c->free_imax_ilen = PETSC_TRUE;
2198: }
2200: /* allocate the matrix space */
2201: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2202: PetscMalloc1(bs2*nz,&c->a);
2203: PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));
2204: c->i = a->i;
2205: c->j = a->j;
2206: c->singlemalloc = PETSC_FALSE;
2207: c->free_a = PETSC_TRUE;
2208: c->free_ij = PETSC_FALSE;
2209: c->parent = A;
2210: PetscObjectReference((PetscObject)A);
2211: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2212: MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2213: } else {
2214: PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
2215: PetscArraycpy(c->i,a->i,mbs+1);
2216: PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
2217: c->singlemalloc = PETSC_TRUE;
2218: c->free_a = PETSC_TRUE;
2219: c->free_ij = PETSC_TRUE;
2220: }
2221: if (mbs > 0) {
2222: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2223: PetscArraycpy(c->j,a->j,nz);
2224: }
2225: if (cpvalues == MAT_COPY_VALUES) {
2226: PetscArraycpy(c->a,a->a,bs2*nz);
2227: } else {
2228: PetscArrayzero(c->a,bs2*nz);
2229: }
2230: if (a->jshort) {
2231: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2232: /* if the parent matrix is reassembled, this child matrix will never notice */
2233: PetscMalloc1(nz,&c->jshort);
2234: PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));
2235: PetscArraycpy(c->jshort,a->jshort,nz);
2237: c->free_jshort = PETSC_TRUE;
2238: }
2239: }
2241: c->roworiented = a->roworiented;
2242: c->nonew = a->nonew;
2244: if (a->diag) {
2245: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2246: c->diag = a->diag;
2247: c->free_diag = PETSC_FALSE;
2248: } else {
2249: PetscMalloc1(mbs,&c->diag);
2250: PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));
2251: for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2252: c->free_diag = PETSC_TRUE;
2253: }
2254: }
2255: c->nz = a->nz;
2256: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2257: c->solve_work = NULL;
2258: c->mult_work = NULL;
2260: *B = C;
2261: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
2262: return(0);
2263: }
2265: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2266: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2268: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)
2269: {
2271: PetscBool isbinary;
2274: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2275: 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);
2276: MatLoad_SeqSBAIJ_Binary(mat,viewer);
2277: return(0);
2278: }
2280: /*@
2281: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2282: (upper triangular entries in CSR format) provided by the user.
2284: Collective
2286: Input Parameters:
2287: + comm - must be an MPI communicator of size 1
2288: . bs - size of block
2289: . m - number of rows
2290: . n - number of columns
2291: . 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
2292: . j - column indices
2293: - a - matrix values
2295: Output Parameter:
2296: . mat - the matrix
2298: Level: advanced
2300: Notes:
2301: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2302: once the matrix is destroyed
2304: You cannot set new nonzero locations into this matrix, that will generate an error.
2306: The i and j indices are 0 based
2308: 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
2309: it is the regular CSR format excluding the lower triangular elements.
2311: .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2313: @*/
2314: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2315: {
2317: PetscInt ii;
2318: Mat_SeqSBAIJ *sbaij;
2321: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2322: if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2324: MatCreate(comm,mat);
2325: MatSetSizes(*mat,m,n,m,n);
2326: MatSetType(*mat,MATSEQSBAIJ);
2327: MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);
2328: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2329: PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);
2330: PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));
2332: sbaij->i = i;
2333: sbaij->j = j;
2334: sbaij->a = a;
2336: sbaij->singlemalloc = PETSC_FALSE;
2337: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2338: sbaij->free_a = PETSC_FALSE;
2339: sbaij->free_ij = PETSC_FALSE;
2340: sbaij->free_imax_ilen = PETSC_TRUE;
2342: for (ii=0; ii<m; ii++) {
2343: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2344: 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]);
2345: }
2346: if (PetscDefined(USE_DEBUG)) {
2347: for (ii=0; ii<sbaij->i[m]; ii++) {
2348: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2349: 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]);
2350: }
2351: }
2353: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2354: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2355: return(0);
2356: }
2358: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2359: {
2361: PetscMPIInt size;
2364: MPI_Comm_size(comm,&size);
2365: if (size == 1 && scall == MAT_REUSE_MATRIX) {
2366: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2367: } else {
2368: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);
2369: }
2370: return(0);
2371: }