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: }
591: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
592: {
593: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
594: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
595: PetscInt *ai = a->i,*ailen = a->ilen;
596: PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
597: MatScalar *ap,*aa = a->a;
600: for (k=0; k<m; k++) { /* loop over rows */
601: row = im[k]; brow = row/bs;
602: if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
603: 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);
604: rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
605: nrow = ailen[brow];
606: for (l=0; l<n; l++) { /* loop over columns */
607: if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
608: 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);
609: col = in[l];
610: bcol = col/bs;
611: cidx = col%bs;
612: ridx = row%bs;
613: high = nrow;
614: low = 0; /* assume unsorted */
615: while (high-low > 5) {
616: t = (low+high)/2;
617: if (rp[t] > bcol) high = t;
618: else low = t;
619: }
620: for (i=low; i<high; i++) {
621: if (rp[i] > bcol) break;
622: if (rp[i] == bcol) {
623: *v++ = ap[bs2*i+bs*cidx+ridx];
624: goto finished;
625: }
626: }
627: *v++ = 0.0;
628: finished:;
629: }
630: }
631: return(0);
632: }
634: PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B)
635: {
636: Mat C;
640: MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C);
641: MatPermute(C,rowp,colp,B);
642: MatDestroy(&C);
643: if (rowp == colp) {
644: MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B);
645: }
646: return(0);
647: }
649: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
650: {
651: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
652: PetscErrorCode ierr;
653: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
654: PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
655: PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
656: PetscBool roworiented=a->roworiented;
657: const PetscScalar *value = v;
658: MatScalar *ap,*aa = a->a,*bap;
661: if (roworiented) stepval = (n-1)*bs;
662: else stepval = (m-1)*bs;
664: for (k=0; k<m; k++) { /* loop over added rows */
665: row = im[k];
666: if (row < 0) continue;
667: 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);
668: rp = aj + ai[row];
669: ap = aa + bs2*ai[row];
670: rmax = imax[row];
671: nrow = ailen[row];
672: low = 0;
673: high = nrow;
674: for (l=0; l<n; l++) { /* loop over added columns */
675: if (in[l] < 0) continue;
676: col = in[l];
677: 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);
678: if (col < row) {
679: if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
680: 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)");
681: }
682: if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
683: else value = v + l*(stepval+bs)*bs + k*bs;
685: if (col <= lastcol) low = 0;
686: else high = nrow;
688: lastcol = col;
689: while (high-low > 7) {
690: t = (low+high)/2;
691: if (rp[t] > col) high = t;
692: else low = t;
693: }
694: for (i=low; i<high; i++) {
695: if (rp[i] > col) break;
696: if (rp[i] == col) {
697: bap = ap + bs2*i;
698: if (roworiented) {
699: if (is == ADD_VALUES) {
700: for (ii=0; ii<bs; ii++,value+=stepval) {
701: for (jj=ii; jj<bs2; jj+=bs) {
702: bap[jj] += *value++;
703: }
704: }
705: } else {
706: for (ii=0; ii<bs; ii++,value+=stepval) {
707: for (jj=ii; jj<bs2; jj+=bs) {
708: bap[jj] = *value++;
709: }
710: }
711: }
712: } else {
713: if (is == ADD_VALUES) {
714: for (ii=0; ii<bs; ii++,value+=stepval) {
715: for (jj=0; jj<bs; jj++) {
716: *bap++ += *value++;
717: }
718: }
719: } else {
720: for (ii=0; ii<bs; ii++,value+=stepval) {
721: for (jj=0; jj<bs; jj++) {
722: *bap++ = *value++;
723: }
724: }
725: }
726: }
727: goto noinsert2;
728: }
729: }
730: if (nonew == 1) goto noinsert2;
731: 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);
732: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
733: N = nrow++ - 1; high++;
734: /* shift up all the later entries in this row */
735: PetscArraymove(rp+i+1,rp+i,N-i+1);
736: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
737: PetscArrayzero(ap+bs2*i,bs2);
738: rp[i] = col;
739: bap = ap + bs2*i;
740: if (roworiented) {
741: for (ii=0; ii<bs; ii++,value+=stepval) {
742: for (jj=ii; jj<bs2; jj+=bs) {
743: bap[jj] = *value++;
744: }
745: }
746: } else {
747: for (ii=0; ii<bs; ii++,value+=stepval) {
748: for (jj=0; jj<bs; jj++) {
749: *bap++ = *value++;
750: }
751: }
752: }
753: noinsert2:;
754: low = i;
755: }
756: ailen[row] = nrow;
757: }
758: return(0);
759: }
761: /*
762: This is not yet used
763: */
764: PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
765: {
766: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
768: const PetscInt *ai = a->i, *aj = a->j,*cols;
769: PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
770: PetscBool flag;
773: PetscMalloc1(m,&ns);
774: while (i < m) {
775: nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
776: /* Limits the number of elements in a node to 'a->inode.limit' */
777: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
778: nzy = ai[j+1] - ai[j];
779: if (nzy != (nzx - j + i)) break;
780: PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);
781: if (!flag) break;
782: }
783: ns[node_count++] = blk_size;
785: i = j;
786: }
787: if (!a->inode.size && m && node_count > .9*m) {
788: PetscFree(ns);
789: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
790: } else {
791: a->inode.node_count = node_count;
793: PetscMalloc1(node_count,&a->inode.size);
794: PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));
795: PetscArraycpy(a->inode.size,ns,node_count);
796: PetscFree(ns);
797: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
799: /* count collections of adjacent columns in each inode */
800: row = 0;
801: cnt = 0;
802: for (i=0; i<node_count; i++) {
803: cols = aj + ai[row] + a->inode.size[i];
804: nz = ai[row+1] - ai[row] - a->inode.size[i];
805: for (j=1; j<nz; j++) {
806: if (cols[j] != cols[j-1]+1) cnt++;
807: }
808: cnt++;
809: row += a->inode.size[i];
810: }
811: PetscMalloc1(2*cnt,&counts);
812: cnt = 0;
813: row = 0;
814: for (i=0; i<node_count; i++) {
815: cols = aj + ai[row] + a->inode.size[i];
816: counts[2*cnt] = cols[0];
817: nz = ai[row+1] - ai[row] - a->inode.size[i];
818: cnt2 = 1;
819: for (j=1; j<nz; j++) {
820: if (cols[j] != cols[j-1]+1) {
821: counts[2*(cnt++)+1] = cnt2;
822: counts[2*cnt] = cols[j];
823: cnt2 = 1;
824: } else cnt2++;
825: }
826: counts[2*(cnt++)+1] = cnt2;
827: row += a->inode.size[i];
828: }
829: PetscIntView(2*cnt,counts,NULL);
830: }
831: return(0);
832: }
834: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
835: {
836: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
838: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
839: PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
840: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
841: MatScalar *aa = a->a,*ap;
844: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
846: if (m) rmax = ailen[0];
847: for (i=1; i<mbs; i++) {
848: /* move each row back by the amount of empty slots (fshift) before it*/
849: fshift += imax[i-1] - ailen[i-1];
850: rmax = PetscMax(rmax,ailen[i]);
851: if (fshift) {
852: ip = aj + ai[i];
853: ap = aa + bs2*ai[i];
854: N = ailen[i];
855: PetscArraymove(ip-fshift,ip,N);
856: PetscArraymove(ap-bs2*fshift,ap,bs2*N);
857: }
858: ai[i] = ai[i-1] + ailen[i-1];
859: }
860: if (mbs) {
861: fshift += imax[mbs-1] - ailen[mbs-1];
862: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
863: }
864: /* reset ilen and imax for each row */
865: for (i=0; i<mbs; i++) {
866: ailen[i] = imax[i] = ai[i+1] - ai[i];
867: }
868: a->nz = ai[mbs];
870: /* diagonals may have moved, reset it */
871: if (a->diag) {
872: PetscArraycpy(a->diag,ai,mbs);
873: }
874: 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);
876: 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);
877: PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
878: PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);
880: A->info.mallocs += a->reallocs;
881: a->reallocs = 0;
882: A->info.nz_unneeded = (PetscReal)fshift*bs2;
883: a->idiagvalid = PETSC_FALSE;
884: a->rmax = rmax;
886: if (A->cmap->n < 65536 && A->cmap->bs == 1) {
887: if (a->jshort && a->free_jshort) {
888: /* when matrix data structure is changed, previous jshort must be replaced */
889: PetscFree(a->jshort);
890: }
891: PetscMalloc1(a->i[A->rmap->n],&a->jshort);
892: PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));
893: for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
894: A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
895: A->ops->sor = MatSOR_SeqSBAIJ_ushort;
896: a->free_jshort = PETSC_TRUE;
897: }
898: return(0);
899: }
901: /*
902: This function returns an array of flags which indicate the locations of contiguous
903: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
904: then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
905: Assume: sizes should be long enough to hold all the values.
906: */
907: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
908: {
909: PetscInt i,j,k,row;
910: PetscBool flg;
913: for (i=0,j=0; i<n; j++) {
914: row = idx[i];
915: if (row%bs!=0) { /* Not the begining of a block */
916: sizes[j] = 1;
917: i++;
918: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
919: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
920: i++;
921: } else { /* Begining of the block, so check if the complete block exists */
922: flg = PETSC_TRUE;
923: for (k=1; k<bs; k++) {
924: if (row+k != idx[i+k]) { /* break in the block */
925: flg = PETSC_FALSE;
926: break;
927: }
928: }
929: if (flg) { /* No break in the bs */
930: sizes[j] = bs;
931: i += bs;
932: } else {
933: sizes[j] = 1;
934: i++;
935: }
936: }
937: }
938: *bs_max = j;
939: return(0);
940: }
943: /* Only add/insert a(i,j) with i<=j (blocks).
944: Any a(i,j) with i>j input by user is ingored.
945: */
947: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
948: {
949: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
951: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
952: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
953: PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
954: PetscInt ridx,cidx,bs2=a->bs2;
955: MatScalar *ap,value,*aa=a->a,*bap;
958: for (k=0; k<m; k++) { /* loop over added rows */
959: row = im[k]; /* row number */
960: brow = row/bs; /* block row number */
961: if (row < 0) continue;
962: 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);
963: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
964: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
965: rmax = imax[brow]; /* maximum space allocated for this row */
966: nrow = ailen[brow]; /* actual length of this row */
967: low = 0;
968: high = nrow;
969: for (l=0; l<n; l++) { /* loop over added columns */
970: if (in[l] < 0) continue;
971: 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);
972: col = in[l];
973: bcol = col/bs; /* block col number */
975: if (brow > bcol) {
976: if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
977: 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)");
978: }
980: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
981: if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
982: /* element value a(k,l) */
983: if (roworiented) value = v[l + k*n];
984: else value = v[k + l*m];
986: /* move pointer bap to a(k,l) quickly and add/insert value */
987: if (col <= lastcol) low = 0;
988: else high = nrow;
990: lastcol = col;
991: while (high-low > 7) {
992: t = (low+high)/2;
993: if (rp[t] > bcol) high = t;
994: else low = t;
995: }
996: for (i=low; i<high; i++) {
997: if (rp[i] > bcol) break;
998: if (rp[i] == bcol) {
999: bap = ap + bs2*i + bs*cidx + ridx;
1000: if (is == ADD_VALUES) *bap += value;
1001: else *bap = value;
1002: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1003: if (brow == bcol && ridx < cidx) {
1004: bap = ap + bs2*i + bs*ridx + cidx;
1005: if (is == ADD_VALUES) *bap += value;
1006: else *bap = value;
1007: }
1008: goto noinsert1;
1009: }
1010: }
1012: if (nonew == 1) goto noinsert1;
1013: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
1014: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
1016: N = nrow++ - 1; high++;
1017: /* shift up all the later entries in this row */
1018: PetscArraymove(rp+i+1,rp+i,N-i+1);
1019: PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
1020: PetscArrayzero(ap+bs2*i,bs2);
1021: rp[i] = bcol;
1022: ap[bs2*i + bs*cidx + ridx] = value;
1023: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1024: if (brow == bcol && ridx < cidx) {
1025: ap[bs2*i + bs*ridx + cidx] = value;
1026: }
1027: A->nonzerostate++;
1028: noinsert1:;
1029: low = i;
1030: }
1031: } /* end of loop over added columns */
1032: ailen[brow] = nrow;
1033: } /* end of loop over added rows */
1034: return(0);
1035: }
1037: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
1038: {
1039: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1040: Mat outA;
1042: PetscBool row_identity;
1045: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
1046: ISIdentity(row,&row_identity);
1047: if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1048: 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()! */
1050: outA = inA;
1051: inA->factortype = MAT_FACTOR_ICC;
1052: PetscFree(inA->solvertype);
1053: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
1055: MatMarkDiagonal_SeqSBAIJ(inA);
1056: MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);
1058: PetscObjectReference((PetscObject)row);
1059: ISDestroy(&a->row);
1060: a->row = row;
1061: PetscObjectReference((PetscObject)row);
1062: ISDestroy(&a->col);
1063: a->col = row;
1065: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1066: if (a->icol) {ISInvertPermutation(row,PETSC_DECIDE, &a->icol);}
1067: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
1069: if (!a->solve_work) {
1070: PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
1071: PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
1072: }
1074: MatCholeskyFactorNumeric(outA,inA,info);
1075: return(0);
1076: }
1078: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1079: {
1080: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1081: PetscInt i,nz,n;
1085: nz = baij->maxnz;
1086: n = mat->cmap->n;
1087: for (i=0; i<nz; i++) baij->j[i] = indices[i];
1089: baij->nz = nz;
1090: for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1092: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1093: return(0);
1094: }
1096: /*@
1097: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1098: in the matrix.
1100: Input Parameters:
1101: + mat - the SeqSBAIJ matrix
1102: - indices - the column indices
1104: Level: advanced
1106: Notes:
1107: This can be called if you have precomputed the nonzero structure of the
1108: matrix and want to provide it to the matrix object to improve the performance
1109: of the MatSetValues() operation.
1111: You MUST have set the correct numbers of nonzeros per row in the call to
1112: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1114: MUST be called before any calls to MatSetValues()
1116: .seealso: MatCreateSeqSBAIJ
1117: @*/
1118: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1119: {
1125: PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
1126: return(0);
1127: }
1129: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1130: {
1132: PetscBool isbaij;
1135: PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");
1136: if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name);
1137: /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1138: if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1139: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1140: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1142: 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");
1143: if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different");
1144: if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size");
1145: PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);
1146: PetscObjectStateIncrease((PetscObject)B);
1147: } else {
1148: MatGetRowUpperTriangular(A);
1149: MatCopy_Basic(A,B,str);
1150: MatRestoreRowUpperTriangular(A);
1151: }
1152: return(0);
1153: }
1155: PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1156: {
1160: MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);
1161: return(0);
1162: }
1164: static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1165: {
1166: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1169: *array = a->a;
1170: return(0);
1171: }
1173: static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1174: {
1176: *array = NULL;
1177: return(0);
1178: }
1180: PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1181: {
1182: PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1183: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1184: Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1188: /* Set the number of nonzeros in the new matrix */
1189: MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
1190: return(0);
1191: }
1193: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1194: {
1195: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1197: PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1198: PetscBLASInt one = 1;
1201: if (str == SAME_NONZERO_PATTERN) {
1202: PetscScalar alpha = a;
1203: PetscBLASInt bnz;
1204: PetscBLASIntCast(x->nz*bs2,&bnz);
1205: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1206: PetscObjectStateIncrease((PetscObject)Y);
1207: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1208: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);
1209: MatAXPY_Basic(Y,a,X,str);
1210: MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);
1211: } else {
1212: Mat B;
1213: PetscInt *nnz;
1214: if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
1215: MatGetRowUpperTriangular(X);
1216: MatGetRowUpperTriangular(Y);
1217: PetscMalloc1(Y->rmap->N,&nnz);
1218: MatCreate(PetscObjectComm((PetscObject)Y),&B);
1219: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
1220: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
1221: MatSetBlockSizesFromMats(B,Y,Y);
1222: MatSetType(B,((PetscObject)Y)->type_name);
1223: MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);
1224: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1226: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1228: MatHeaderReplace(Y,&B);
1229: PetscFree(nnz);
1230: MatRestoreRowUpperTriangular(X);
1231: MatRestoreRowUpperTriangular(Y);
1232: }
1233: return(0);
1234: }
1236: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1237: {
1239: *flg = PETSC_TRUE;
1240: return(0);
1241: }
1243: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1244: {
1246: *flg = PETSC_TRUE;
1247: return(0);
1248: }
1250: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1251: {
1253: *flg = PETSC_FALSE;
1254: return(0);
1255: }
1257: PetscErrorCode MatConjugate_SeqSBAIJ(Mat A)
1258: {
1259: #if defined(PETSC_USE_COMPLEX)
1260: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1261: PetscInt i,nz = a->bs2*a->i[a->mbs];
1262: MatScalar *aa = a->a;
1265: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1266: #else
1268: #endif
1269: return(0);
1270: }
1272: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1273: {
1274: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1275: PetscInt i,nz = a->bs2*a->i[a->mbs];
1276: MatScalar *aa = a->a;
1279: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1280: return(0);
1281: }
1283: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1284: {
1285: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1286: PetscInt i,nz = a->bs2*a->i[a->mbs];
1287: MatScalar *aa = a->a;
1290: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1291: return(0);
1292: }
1294: PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1295: {
1296: Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1297: PetscErrorCode ierr;
1298: PetscInt i,j,k,count;
1299: PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1300: PetscScalar zero = 0.0;
1301: MatScalar *aa;
1302: const PetscScalar *xx;
1303: PetscScalar *bb;
1304: PetscBool *zeroed,vecs = PETSC_FALSE;
1307: /* fix right hand side if needed */
1308: if (x && b) {
1309: VecGetArrayRead(x,&xx);
1310: VecGetArray(b,&bb);
1311: vecs = PETSC_TRUE;
1312: }
1314: /* zero the columns */
1315: PetscCalloc1(A->rmap->n,&zeroed);
1316: for (i=0; i<is_n; i++) {
1317: 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]);
1318: zeroed[is_idx[i]] = PETSC_TRUE;
1319: }
1320: if (vecs) {
1321: for (i=0; i<A->rmap->N; i++) {
1322: row = i/bs;
1323: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1324: for (k=0; k<bs; k++) {
1325: col = bs*baij->j[j] + k;
1326: if (col <= i) continue;
1327: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1328: if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1329: if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1330: }
1331: }
1332: }
1333: for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1334: }
1336: for (i=0; i<A->rmap->N; i++) {
1337: if (!zeroed[i]) {
1338: row = i/bs;
1339: for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1340: for (k=0; k<bs; k++) {
1341: col = bs*baij->j[j] + k;
1342: if (zeroed[col]) {
1343: aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1344: aa[0] = 0.0;
1345: }
1346: }
1347: }
1348: }
1349: }
1350: PetscFree(zeroed);
1351: if (vecs) {
1352: VecRestoreArrayRead(x,&xx);
1353: VecRestoreArray(b,&bb);
1354: }
1356: /* zero the rows */
1357: for (i=0; i<is_n; i++) {
1358: row = is_idx[i];
1359: count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1360: aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1361: for (k=0; k<count; k++) {
1362: aa[0] = zero;
1363: aa += bs;
1364: }
1365: if (diag != 0.0) {
1366: (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
1367: }
1368: }
1369: MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);
1370: return(0);
1371: }
1373: PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1374: {
1376: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1379: if (!Y->preallocated || !aij->nz) {
1380: MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
1381: }
1382: MatShift_Basic(Y,a);
1383: return(0);
1384: }
1386: /* -------------------------------------------------------------------*/
1387: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1388: MatGetRow_SeqSBAIJ,
1389: MatRestoreRow_SeqSBAIJ,
1390: MatMult_SeqSBAIJ_N,
1391: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1392: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1393: MatMultAdd_SeqSBAIJ_N,
1394: NULL,
1395: NULL,
1396: NULL,
1397: /* 10*/ NULL,
1398: NULL,
1399: MatCholeskyFactor_SeqSBAIJ,
1400: MatSOR_SeqSBAIJ,
1401: MatTranspose_SeqSBAIJ,
1402: /* 15*/ MatGetInfo_SeqSBAIJ,
1403: MatEqual_SeqSBAIJ,
1404: MatGetDiagonal_SeqSBAIJ,
1405: MatDiagonalScale_SeqSBAIJ,
1406: MatNorm_SeqSBAIJ,
1407: /* 20*/ NULL,
1408: MatAssemblyEnd_SeqSBAIJ,
1409: MatSetOption_SeqSBAIJ,
1410: MatZeroEntries_SeqSBAIJ,
1411: /* 24*/ NULL,
1412: NULL,
1413: NULL,
1414: NULL,
1415: NULL,
1416: /* 29*/ MatSetUp_SeqSBAIJ,
1417: NULL,
1418: NULL,
1419: NULL,
1420: NULL,
1421: /* 34*/ MatDuplicate_SeqSBAIJ,
1422: NULL,
1423: NULL,
1424: NULL,
1425: MatICCFactor_SeqSBAIJ,
1426: /* 39*/ MatAXPY_SeqSBAIJ,
1427: MatCreateSubMatrices_SeqSBAIJ,
1428: MatIncreaseOverlap_SeqSBAIJ,
1429: MatGetValues_SeqSBAIJ,
1430: MatCopy_SeqSBAIJ,
1431: /* 44*/ NULL,
1432: MatScale_SeqSBAIJ,
1433: MatShift_SeqSBAIJ,
1434: NULL,
1435: MatZeroRowsColumns_SeqSBAIJ,
1436: /* 49*/ NULL,
1437: MatGetRowIJ_SeqSBAIJ,
1438: MatRestoreRowIJ_SeqSBAIJ,
1439: NULL,
1440: NULL,
1441: /* 54*/ NULL,
1442: NULL,
1443: NULL,
1444: MatPermute_SeqSBAIJ,
1445: MatSetValuesBlocked_SeqSBAIJ,
1446: /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1447: NULL,
1448: NULL,
1449: NULL,
1450: NULL,
1451: /* 64*/ NULL,
1452: NULL,
1453: NULL,
1454: NULL,
1455: NULL,
1456: /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1457: NULL,
1458: MatConvert_MPISBAIJ_Basic,
1459: NULL,
1460: NULL,
1461: /* 74*/ NULL,
1462: NULL,
1463: NULL,
1464: NULL,
1465: NULL,
1466: /* 79*/ NULL,
1467: NULL,
1468: NULL,
1469: MatGetInertia_SeqSBAIJ,
1470: MatLoad_SeqSBAIJ,
1471: /* 84*/ MatIsSymmetric_SeqSBAIJ,
1472: MatIsHermitian_SeqSBAIJ,
1473: MatIsStructurallySymmetric_SeqSBAIJ,
1474: NULL,
1475: NULL,
1476: /* 89*/ NULL,
1477: NULL,
1478: NULL,
1479: NULL,
1480: NULL,
1481: /* 94*/ NULL,
1482: NULL,
1483: NULL,
1484: NULL,
1485: NULL,
1486: /* 99*/ NULL,
1487: NULL,
1488: NULL,
1489: MatConjugate_SeqSBAIJ,
1490: NULL,
1491: /*104*/ NULL,
1492: MatRealPart_SeqSBAIJ,
1493: MatImaginaryPart_SeqSBAIJ,
1494: MatGetRowUpperTriangular_SeqSBAIJ,
1495: MatRestoreRowUpperTriangular_SeqSBAIJ,
1496: /*109*/ NULL,
1497: NULL,
1498: NULL,
1499: NULL,
1500: MatMissingDiagonal_SeqSBAIJ,
1501: /*114*/ NULL,
1502: NULL,
1503: NULL,
1504: NULL,
1505: NULL,
1506: /*119*/ NULL,
1507: NULL,
1508: NULL,
1509: NULL,
1510: NULL,
1511: /*124*/ NULL,
1512: NULL,
1513: NULL,
1514: NULL,
1515: NULL,
1516: /*129*/ NULL,
1517: NULL,
1518: NULL,
1519: NULL,
1520: NULL,
1521: /*134*/ NULL,
1522: NULL,
1523: NULL,
1524: NULL,
1525: NULL,
1526: /*139*/ MatSetBlockSizes_Default,
1527: NULL,
1528: NULL,
1529: NULL,
1530: NULL,
1531: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ
1532: };
1534: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1535: {
1536: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1537: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1541: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1543: /* allocate space for values if not already there */
1544: if (!aij->saved_values) {
1545: PetscMalloc1(nz+1,&aij->saved_values);
1546: }
1548: /* copy values over */
1549: PetscArraycpy(aij->saved_values,aij->a,nz);
1550: return(0);
1551: }
1553: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1554: {
1555: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1557: PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1560: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1561: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1563: /* copy values over */
1564: PetscArraycpy(aij->a,aij->saved_values,nz);
1565: return(0);
1566: }
1568: static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1569: {
1570: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1572: PetscInt i,mbs,nbs,bs2;
1573: PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1576: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1578: MatSetBlockSize(B,PetscAbs(bs));
1579: PetscLayoutSetUp(B->rmap);
1580: PetscLayoutSetUp(B->cmap);
1581: 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);
1582: PetscLayoutGetBlockSize(B->rmap,&bs);
1584: B->preallocated = PETSC_TRUE;
1586: mbs = B->rmap->N/bs;
1587: nbs = B->cmap->n/bs;
1588: bs2 = bs*bs;
1590: 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");
1592: if (nz == MAT_SKIP_ALLOCATION) {
1593: skipallocation = PETSC_TRUE;
1594: nz = 0;
1595: }
1597: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1598: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1599: if (nnz) {
1600: for (i=0; i<mbs; i++) {
1601: 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]);
1602: 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);
1603: }
1604: }
1606: B->ops->mult = MatMult_SeqSBAIJ_N;
1607: B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1608: B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1609: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1611: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1612: if (!flg) {
1613: switch (bs) {
1614: case 1:
1615: B->ops->mult = MatMult_SeqSBAIJ_1;
1616: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1617: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1618: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1619: break;
1620: case 2:
1621: B->ops->mult = MatMult_SeqSBAIJ_2;
1622: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1623: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1624: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1625: break;
1626: case 3:
1627: B->ops->mult = MatMult_SeqSBAIJ_3;
1628: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1629: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1630: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1631: break;
1632: case 4:
1633: B->ops->mult = MatMult_SeqSBAIJ_4;
1634: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1635: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1636: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1637: break;
1638: case 5:
1639: B->ops->mult = MatMult_SeqSBAIJ_5;
1640: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1641: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1642: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1643: break;
1644: case 6:
1645: B->ops->mult = MatMult_SeqSBAIJ_6;
1646: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1647: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1648: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1649: break;
1650: case 7:
1651: B->ops->mult = MatMult_SeqSBAIJ_7;
1652: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1653: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1654: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1655: break;
1656: }
1657: }
1659: b->mbs = mbs;
1660: b->nbs = nbs;
1661: if (!skipallocation) {
1662: if (!b->imax) {
1663: PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
1665: b->free_imax_ilen = PETSC_TRUE;
1667: PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));
1668: }
1669: if (!nnz) {
1670: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1671: else if (nz <= 0) nz = 1;
1672: nz = PetscMin(nbs,nz);
1673: for (i=0; i<mbs; i++) b->imax[i] = nz;
1674: PetscIntMultError(nz,mbs,&nz);
1675: } else {
1676: PetscInt64 nz64 = 0;
1677: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1678: PetscIntCast(nz64,&nz);
1679: }
1680: /* b->ilen will count nonzeros in each block row so far. */
1681: for (i=0; i<mbs; i++) b->ilen[i] = 0;
1682: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1684: /* allocate the matrix space */
1685: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
1686: PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
1687: PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
1688: PetscArrayzero(b->a,nz*bs2);
1689: PetscArrayzero(b->j,nz);
1691: b->singlemalloc = PETSC_TRUE;
1693: /* pointer to beginning of each row */
1694: b->i[0] = 0;
1695: for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1697: b->free_a = PETSC_TRUE;
1698: b->free_ij = PETSC_TRUE;
1699: } else {
1700: b->free_a = PETSC_FALSE;
1701: b->free_ij = PETSC_FALSE;
1702: }
1704: b->bs2 = bs2;
1705: b->nz = 0;
1706: b->maxnz = nz;
1707: b->inew = NULL;
1708: b->jnew = NULL;
1709: b->anew = NULL;
1710: b->a2anew = NULL;
1711: b->permute = PETSC_FALSE;
1713: B->was_assembled = PETSC_FALSE;
1714: B->assembled = PETSC_FALSE;
1715: if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
1716: return(0);
1717: }
1719: PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1720: {
1721: PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1722: PetscScalar *values=NULL;
1723: PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1727: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
1728: PetscLayoutSetBlockSize(B->rmap,bs);
1729: PetscLayoutSetBlockSize(B->cmap,bs);
1730: PetscLayoutSetUp(B->rmap);
1731: PetscLayoutSetUp(B->cmap);
1732: PetscLayoutGetBlockSize(B->rmap,&bs);
1733: m = B->rmap->n/bs;
1735: if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]);
1736: PetscMalloc1(m+1,&nnz);
1737: for (i=0; i<m; i++) {
1738: nz = ii[i+1] - ii[i];
1739: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz);
1740: anz = 0;
1741: for (j=0; j<nz; j++) {
1742: /* count only values on the diagonal or above */
1743: if (jj[ii[i] + j] >= i) {
1744: anz = nz - j;
1745: break;
1746: }
1747: }
1748: nz_max = PetscMax(nz_max,anz);
1749: nnz[i] = anz;
1750: }
1751: MatSeqSBAIJSetPreallocation(B,bs,0,nnz);
1752: PetscFree(nnz);
1754: values = (PetscScalar*)V;
1755: if (!values) {
1756: PetscCalloc1(bs*bs*nz_max,&values);
1757: }
1758: for (i=0; i<m; i++) {
1759: PetscInt ncols = ii[i+1] - ii[i];
1760: const PetscInt *icols = jj + ii[i];
1761: if (!roworiented || bs == 1) {
1762: const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1763: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
1764: } else {
1765: for (j=0; j<ncols; j++) {
1766: const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1767: MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
1768: }
1769: }
1770: }
1771: if (!V) { PetscFree(values); }
1772: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1773: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1774: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
1775: return(0);
1776: }
1778: /*
1779: This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1780: */
1781: PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1782: {
1784: PetscBool flg = PETSC_FALSE;
1785: PetscInt bs = B->rmap->bs;
1788: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);
1789: if (flg) bs = 8;
1791: if (!natural) {
1792: switch (bs) {
1793: case 1:
1794: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1795: break;
1796: case 2:
1797: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1798: break;
1799: case 3:
1800: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1801: break;
1802: case 4:
1803: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1804: break;
1805: case 5:
1806: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1807: break;
1808: case 6:
1809: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1810: break;
1811: case 7:
1812: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1813: break;
1814: default:
1815: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1816: break;
1817: }
1818: } else {
1819: switch (bs) {
1820: case 1:
1821: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1822: break;
1823: case 2:
1824: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1825: break;
1826: case 3:
1827: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1828: break;
1829: case 4:
1830: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1831: break;
1832: case 5:
1833: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1834: break;
1835: case 6:
1836: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1837: break;
1838: case 7:
1839: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1840: break;
1841: default:
1842: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1843: break;
1844: }
1845: }
1846: return(0);
1847: }
1849: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*);
1850: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*);
1852: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1853: {
1854: PetscInt n = A->rmap->n;
1858: #if defined(PETSC_USE_COMPLEX)
1859: 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");
1860: #endif
1862: MatCreate(PetscObjectComm((PetscObject)A),B);
1863: MatSetSizes(*B,n,n,n,n);
1864: if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1865: MatSetType(*B,MATSEQSBAIJ);
1866: MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);
1868: (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1869: (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1870: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1872: (*B)->factortype = ftype;
1873: (*B)->useordering = PETSC_TRUE;
1874: PetscFree((*B)->solvertype);
1875: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
1876: return(0);
1877: }
1879: /*@C
1880: MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1882: Not Collective
1884: Input Parameter:
1885: . mat - a MATSEQSBAIJ matrix
1887: Output Parameter:
1888: . array - pointer to the data
1890: Level: intermediate
1892: .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1893: @*/
1894: PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1895: {
1899: PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
1900: return(0);
1901: }
1903: /*@C
1904: MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1906: Not Collective
1908: Input Parameters:
1909: + mat - a MATSEQSBAIJ matrix
1910: - array - pointer to the data
1912: Level: intermediate
1914: .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1915: @*/
1916: PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1917: {
1921: PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
1922: return(0);
1923: }
1925: /*MC
1926: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1927: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1929: For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1930: can call MatSetOption(Mat, MAT_HERMITIAN).
1932: Options Database Keys:
1933: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1935: Notes:
1936: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1937: stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1938: 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.
1940: The number of rows in the matrix must be less than or equal to the number of columns
1942: Level: beginner
1944: .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1945: M*/
1946: PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1947: {
1948: Mat_SeqSBAIJ *b;
1950: PetscMPIInt size;
1951: PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1954: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
1955: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1957: PetscNewLog(B,&b);
1958: B->data = (void*)b;
1959: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1961: B->ops->destroy = MatDestroy_SeqSBAIJ;
1962: B->ops->view = MatView_SeqSBAIJ;
1963: b->row = NULL;
1964: b->icol = NULL;
1965: b->reallocs = 0;
1966: b->saved_values = NULL;
1967: b->inode.limit = 5;
1968: b->inode.max_limit = 5;
1970: b->roworiented = PETSC_TRUE;
1971: b->nonew = 0;
1972: b->diag = NULL;
1973: b->solve_work = NULL;
1974: b->mult_work = NULL;
1975: B->spptr = NULL;
1976: B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1977: b->keepnonzeropattern = PETSC_FALSE;
1979: b->inew = NULL;
1980: b->jnew = NULL;
1981: b->anew = NULL;
1982: b->a2anew = NULL;
1983: b->permute = PETSC_FALSE;
1985: b->ignore_ltriangular = PETSC_TRUE;
1987: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);
1989: b->getrow_utriangular = PETSC_FALSE;
1991: PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);
1993: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);
1994: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);
1995: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);
1996: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);
1997: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1998: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);
1999: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);
2000: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);
2001: PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);
2002: #if defined(PETSC_HAVE_ELEMENTAL)
2003: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);
2004: #endif
2005: #if defined(PETSC_HAVE_SCALAPACK)
2006: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);
2007: #endif
2009: B->symmetric = PETSC_TRUE;
2010: B->structurally_symmetric = PETSC_TRUE;
2011: B->symmetric_set = PETSC_TRUE;
2012: B->structurally_symmetric_set = PETSC_TRUE;
2013: B->symmetric_eternal = PETSC_TRUE;
2014: #if defined(PETSC_USE_COMPLEX)
2015: B->hermitian = PETSC_FALSE;
2016: B->hermitian_set = PETSC_FALSE;
2017: #else
2018: B->hermitian = PETSC_TRUE;
2019: B->hermitian_set = PETSC_TRUE;
2020: #endif
2022: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
2024: PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");
2025: PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);
2026: if (no_unroll) {
2027: PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");
2028: }
2029: PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);
2030: if (no_inode) {
2031: PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");
2032: }
2033: PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);
2034: PetscOptionsEnd();
2035: b->inode.use = (PetscBool)(!(no_unroll || no_inode));
2036: if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
2037: return(0);
2038: }
2040: /*@C
2041: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
2042: compressed row) format. For good matrix assembly performance the
2043: user should preallocate the matrix storage by setting the parameter nz
2044: (or the array nnz). By setting these parameters accurately, performance
2045: during matrix assembly can be increased by more than a factor of 50.
2047: Collective on Mat
2049: Input Parameters:
2050: + B - the symmetric matrix
2051: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2052: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2053: . nz - number of block nonzeros per block row (same for all rows)
2054: - nnz - array containing the number of block nonzeros in the upper triangular plus
2055: diagonal portion of each block (possibly different for each block row) or NULL
2057: Options Database Keys:
2058: + -mat_no_unroll - uses code that does not unroll the loops in the
2059: block calculations (much slower)
2060: - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2062: Level: intermediate
2064: Notes:
2065: Specify the preallocated storage with either nz or nnz (not both).
2066: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2067: allocation. See Users-Manual: ch_mat for details.
2069: You can call MatGetInfo() to get information on how effective the preallocation was;
2070: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2071: You can also run with the option -info and look for messages with the string
2072: malloc in them to see if additional memory allocation was needed.
2074: If the nnz parameter is given then the nz parameter is ignored
2077: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2078: @*/
2079: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2080: {
2087: PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
2088: return(0);
2089: }
2091: /*@C
2092: MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2094: Input Parameters:
2095: + B - the matrix
2096: . bs - size of block, the blocks are ALWAYS square.
2097: . i - the indices into j for the start of each local row (starts with zero)
2098: . j - the column indices for each local row (starts with zero) these must be sorted for each row
2099: - v - optional values in the matrix
2101: Level: advanced
2103: Notes:
2104: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2105: may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2106: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2107: MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2108: block column and the second index is over columns within a block.
2110: Any entries below the diagonal are ignored
2112: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2113: and usually the numerical values as well
2115: .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2116: @*/
2117: PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2118: {
2125: PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
2126: return(0);
2127: }
2129: /*@C
2130: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2131: compressed row) format. For good matrix assembly performance the
2132: user should preallocate the matrix storage by setting the parameter nz
2133: (or the array nnz). By setting these parameters accurately, performance
2134: during matrix assembly can be increased by more than a factor of 50.
2136: Collective
2138: Input Parameters:
2139: + comm - MPI communicator, set to PETSC_COMM_SELF
2140: . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2141: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2142: . m - number of rows, or number of columns
2143: . nz - number of block nonzeros per block row (same for all rows)
2144: - nnz - array containing the number of block nonzeros in the upper triangular plus
2145: diagonal portion of each block (possibly different for each block row) or NULL
2147: Output Parameter:
2148: . A - the symmetric matrix
2150: Options Database Keys:
2151: + -mat_no_unroll - uses code that does not unroll the loops in the
2152: block calculations (much slower)
2153: - -mat_block_size - size of the blocks to use
2155: Level: intermediate
2157: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2158: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2159: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2161: Notes:
2162: The number of rows and columns must be divisible by blocksize.
2163: This matrix type does not support complex Hermitian operation.
2165: Specify the preallocated storage with either nz or nnz (not both).
2166: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2167: allocation. See Users-Manual: ch_mat for details.
2169: If the nnz parameter is given then the nz parameter is ignored
2171: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2172: @*/
2173: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2174: {
2178: MatCreate(comm,A);
2179: MatSetSizes(*A,m,n,m,n);
2180: MatSetType(*A,MATSEQSBAIJ);
2181: MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
2182: return(0);
2183: }
2185: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2186: {
2187: Mat C;
2188: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2190: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2193: if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");
2195: *B = NULL;
2196: MatCreate(PetscObjectComm((PetscObject)A),&C);
2197: MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
2198: MatSetBlockSizesFromMats(C,A,A);
2199: MatSetType(C,MATSEQSBAIJ);
2200: c = (Mat_SeqSBAIJ*)C->data;
2202: C->preallocated = PETSC_TRUE;
2203: C->factortype = A->factortype;
2204: c->row = NULL;
2205: c->icol = NULL;
2206: c->saved_values = NULL;
2207: c->keepnonzeropattern = a->keepnonzeropattern;
2208: C->assembled = PETSC_TRUE;
2210: PetscLayoutReference(A->rmap,&C->rmap);
2211: PetscLayoutReference(A->cmap,&C->cmap);
2212: c->bs2 = a->bs2;
2213: c->mbs = a->mbs;
2214: c->nbs = a->nbs;
2216: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2217: c->imax = a->imax;
2218: c->ilen = a->ilen;
2219: c->free_imax_ilen = PETSC_FALSE;
2220: } else {
2221: PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);
2222: PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));
2223: for (i=0; i<mbs; i++) {
2224: c->imax[i] = a->imax[i];
2225: c->ilen[i] = a->ilen[i];
2226: }
2227: c->free_imax_ilen = PETSC_TRUE;
2228: }
2230: /* allocate the matrix space */
2231: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2232: PetscMalloc1(bs2*nz,&c->a);
2233: PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));
2234: c->i = a->i;
2235: c->j = a->j;
2236: c->singlemalloc = PETSC_FALSE;
2237: c->free_a = PETSC_TRUE;
2238: c->free_ij = PETSC_FALSE;
2239: c->parent = A;
2240: PetscObjectReference((PetscObject)A);
2241: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2242: MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2243: } else {
2244: PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
2245: PetscArraycpy(c->i,a->i,mbs+1);
2246: PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
2247: c->singlemalloc = PETSC_TRUE;
2248: c->free_a = PETSC_TRUE;
2249: c->free_ij = PETSC_TRUE;
2250: }
2251: if (mbs > 0) {
2252: if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2253: PetscArraycpy(c->j,a->j,nz);
2254: }
2255: if (cpvalues == MAT_COPY_VALUES) {
2256: PetscArraycpy(c->a,a->a,bs2*nz);
2257: } else {
2258: PetscArrayzero(c->a,bs2*nz);
2259: }
2260: if (a->jshort) {
2261: /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2262: /* if the parent matrix is reassembled, this child matrix will never notice */
2263: PetscMalloc1(nz,&c->jshort);
2264: PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));
2265: PetscArraycpy(c->jshort,a->jshort,nz);
2267: c->free_jshort = PETSC_TRUE;
2268: }
2269: }
2271: c->roworiented = a->roworiented;
2272: c->nonew = a->nonew;
2274: if (a->diag) {
2275: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2276: c->diag = a->diag;
2277: c->free_diag = PETSC_FALSE;
2278: } else {
2279: PetscMalloc1(mbs,&c->diag);
2280: PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));
2281: for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2282: c->free_diag = PETSC_TRUE;
2283: }
2284: }
2285: c->nz = a->nz;
2286: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2287: c->solve_work = NULL;
2288: c->mult_work = NULL;
2290: *B = C;
2291: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
2292: return(0);
2293: }
2295: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2296: #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2298: PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)
2299: {
2301: PetscBool isbinary;
2304: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2305: 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);
2306: MatLoad_SeqSBAIJ_Binary(mat,viewer);
2307: return(0);
2308: }
2310: /*@
2311: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2312: (upper triangular entries in CSR format) provided by the user.
2314: Collective
2316: Input Parameters:
2317: + comm - must be an MPI communicator of size 1
2318: . bs - size of block
2319: . m - number of rows
2320: . n - number of columns
2321: . 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
2322: . j - column indices
2323: - a - matrix values
2325: Output Parameter:
2326: . mat - the matrix
2328: Level: advanced
2330: Notes:
2331: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2332: once the matrix is destroyed
2334: You cannot set new nonzero locations into this matrix, that will generate an error.
2336: The i and j indices are 0 based
2338: 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
2339: it is the regular CSR format excluding the lower triangular elements.
2341: .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2343: @*/
2344: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2345: {
2347: PetscInt ii;
2348: Mat_SeqSBAIJ *sbaij;
2351: if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2352: if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2354: MatCreate(comm,mat);
2355: MatSetSizes(*mat,m,n,m,n);
2356: MatSetType(*mat,MATSEQSBAIJ);
2357: MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);
2358: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2359: PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);
2360: PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));
2362: sbaij->i = i;
2363: sbaij->j = j;
2364: sbaij->a = a;
2366: sbaij->singlemalloc = PETSC_FALSE;
2367: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2368: sbaij->free_a = PETSC_FALSE;
2369: sbaij->free_ij = PETSC_FALSE;
2370: sbaij->free_imax_ilen = PETSC_TRUE;
2372: for (ii=0; ii<m; ii++) {
2373: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2374: 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]);
2375: }
2376: if (PetscDefined(USE_DEBUG)) {
2377: for (ii=0; ii<sbaij->i[m]; ii++) {
2378: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2379: 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]);
2380: }
2381: }
2383: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2384: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2385: return(0);
2386: }
2388: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2389: {
2391: PetscMPIInt size;
2394: MPI_Comm_size(comm,&size);
2395: if (size == 1 && scall == MAT_REUSE_MATRIX) {
2396: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2397: } else {
2398: MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);
2399: }
2400: return(0);
2401: }