Actual source code: sbaij.c
1: #define PETSCMAT_DLL
3: /*
4: Defines the basic matrix operations for the SBAIJ (compressed row)
5: matrix storage format.
6: */
7: #include src/mat/impls/baij/seq/baij.h
8: #include src/inline/spops.h
9: #include src/mat/impls/sbaij/seq/sbaij.h
11: #define CHUNKSIZE 10
13: /*
14: Checks for missing diagonals
15: */
18: PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A)
19: {
20: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
22: PetscInt *diag,*jj = a->j,i;
25: MatMarkDiagonal_SeqSBAIJ(A);
26: diag = a->diag;
27: for (i=0; i<a->mbs; i++) {
28: if (jj[diag[i]] != i) SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Matrix is missing diagonal number %D",i);
29: }
30: return(0);
31: }
35: PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
36: {
37: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
39: PetscInt i;
42: if (!a->diag) {
43: PetscMalloc(a->mbs*sizeof(PetscInt),&a->diag);
44: }
45: for (i=0; i<a->mbs; i++) a->diag[i] = a->i[i];
46: return(0);
47: }
51: static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done)
52: {
53: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
54: PetscInt n = a->mbs,i;
57: *nn = n;
58: if (!ia) return(0);
60: if (oshift == 1) {
61: /* temporarily add 1 to i and j indices */
62: PetscInt nz = a->i[n];
63: for (i=0; i<nz; i++) a->j[i]++;
64: for (i=0; i<n+1; i++) a->i[i]++;
65: *ia = a->i; *ja = a->j;
66: } else {
67: *ia = a->i; *ja = a->j;
68: }
69: return(0);
70: }
74: static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done)
75: {
76: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
77: PetscInt i,n = a->mbs;
80: if (!ia) return(0);
82: if (oshift == 1) {
83: PetscInt nz = a->i[n]-1;
84: for (i=0; i<nz; i++) a->j[i]--;
85: for (i=0; i<n+1; i++) a->i[i]--;
86: }
87: return(0);
88: }
92: PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
93: {
94: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
98: #if defined(PETSC_USE_LOG)
99: PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap.N,a->nz);
100: #endif
101: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
102: if (a->row) {ISDestroy(a->row);}
103: if (a->col){ISDestroy(a->col);}
104: if (a->icol) {ISDestroy(a->icol);}
105: PetscFree(a->diag);
106: PetscFree2(a->imax,a->ilen);
107: PetscFree(a->solve_work);
108: PetscFree(a->relax_work);
109: PetscFree(a->solves_work);
110: PetscFree(a->mult_work);
111: PetscFree(a->saved_values);
112: PetscFree(a->xtoy);
114: PetscFree(a->inew);
115: PetscFree(a);
117: PetscObjectChangeTypeName((PetscObject)A,0);
118: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);
119: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);
120: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);
121: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);
122: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);
123: PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);
124: return(0);
125: }
129: PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op)
130: {
131: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
135: switch (op) {
136: case MAT_ROW_ORIENTED:
137: a->roworiented = PETSC_TRUE;
138: break;
139: case MAT_COLUMN_ORIENTED:
140: a->roworiented = PETSC_FALSE;
141: break;
142: case MAT_COLUMNS_SORTED:
143: a->sorted = PETSC_TRUE;
144: break;
145: case MAT_COLUMNS_UNSORTED:
146: a->sorted = PETSC_FALSE;
147: break;
148: case MAT_KEEP_ZEROED_ROWS:
149: a->keepzeroedrows = PETSC_TRUE;
150: break;
151: case MAT_NO_NEW_NONZERO_LOCATIONS:
152: a->nonew = 1;
153: break;
154: case MAT_NEW_NONZERO_LOCATION_ERR:
155: a->nonew = -1;
156: break;
157: case MAT_NEW_NONZERO_ALLOCATION_ERR:
158: a->nonew = -2;
159: break;
160: case MAT_YES_NEW_NONZERO_LOCATIONS:
161: a->nonew = 0;
162: break;
163: case MAT_ROWS_SORTED:
164: case MAT_ROWS_UNSORTED:
165: case MAT_YES_NEW_DIAGONALS:
166: case MAT_IGNORE_OFF_PROC_ENTRIES:
167: case MAT_USE_HASH_TABLE:
168: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
169: break;
170: case MAT_NO_NEW_DIAGONALS:
171: SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS");
172: case MAT_NOT_SYMMETRIC:
173: case MAT_NOT_STRUCTURALLY_SYMMETRIC:
174: case MAT_HERMITIAN:
175: SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric");
176: case MAT_SYMMETRIC:
177: case MAT_STRUCTURALLY_SYMMETRIC:
178: case MAT_NOT_HERMITIAN:
179: case MAT_SYMMETRY_ETERNAL:
180: case MAT_NOT_SYMMETRY_ETERNAL:
181: PetscInfo1(A,"Option %s not relevent\n",MatOptions[op]);
182: break;
183: case MAT_IGNORE_LOWER_TRIANGULAR:
184: a->ignore_ltriangular = PETSC_TRUE;
185: break;
186: case MAT_ERROR_LOWER_TRIANGULAR:
187: a->ignore_ltriangular = PETSC_FALSE;
188: break;
189: case MAT_GETROW_UPPERTRIANGULAR:
190: a->getrow_utriangular = PETSC_TRUE;
191: break;
192: default:
193: SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op);
194: }
195: return(0);
196: }
200: PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v)
201: {
202: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
204: PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2;
205: MatScalar *aa,*aa_i;
206: PetscScalar *v_i;
209: if (A && !a->getrow_utriangular) SETERRQ(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) or MatGetRowUpperTriangular()");
210: /* Get the upper triangular part of the row */
211: bs = A->rmap.bs;
212: ai = a->i;
213: aj = a->j;
214: aa = a->a;
215: bs2 = a->bs2;
216:
217: if (row < 0 || row >= A->rmap.N) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row);
218:
219: bn = row/bs; /* Block number */
220: bp = row % bs; /* Block position */
221: M = ai[bn+1] - ai[bn];
222: *ncols = bs*M;
223:
224: if (v) {
225: *v = 0;
226: if (*ncols) {
227: PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);
228: for (i=0; i<M; i++) { /* for each block in the block row */
229: v_i = *v + i*bs;
230: aa_i = aa + bs2*(ai[bn] + i);
231: for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];}
232: }
233: }
234: }
235:
236: if (cols) {
237: *cols = 0;
238: if (*ncols) {
239: PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);
240: for (i=0; i<M; i++) { /* for each block in the block row */
241: cols_i = *cols + i*bs;
242: itmp = bs*aj[ai[bn] + i];
243: for (j=0; j<bs; j++) {cols_i[j] = itmp++;}
244: }
245: }
246: }
247:
248: /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */
249: /* this segment is currently removed, so only entries in the upper triangle are obtained */
250: #ifdef column_search
251: v_i = *v + M*bs;
252: cols_i = *cols + M*bs;
253: for (i=0; i<bn; i++){ /* for each block row */
254: M = ai[i+1] - ai[i];
255: for (j=0; j<M; j++){
256: itmp = aj[ai[i] + j]; /* block column value */
257: if (itmp == bn){
258: aa_i = aa + bs2*(ai[i] + j) + bs*bp;
259: for (k=0; k<bs; k++) {
260: *cols_i++ = i*bs+k;
261: *v_i++ = aa_i[k];
262: }
263: *ncols += bs;
264: break;
265: }
266: }
267: }
268: #endif
269: return(0);
270: }
274: PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
275: {
277:
279: if (idx) {PetscFree(*idx);}
280: if (v) {PetscFree(*v);}
281: return(0);
282: }
286: PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
287: {
288: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
291: a->getrow_utriangular = PETSC_TRUE;
292: return(0);
293: }
296: PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
297: {
298: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
301: a->getrow_utriangular = PETSC_FALSE;
302: return(0);
303: }
307: PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,Mat *B)
308: {
311: MatDuplicate(A,MAT_COPY_VALUES,B);
312: return(0);
313: }
317: static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
318: {
319: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
320: PetscErrorCode ierr;
321: PetscInt i,j,bs = A->rmap.bs,k,l,bs2=a->bs2;
322: const char *name;
323: PetscViewerFormat format;
324:
326: PetscObjectGetName((PetscObject)A,&name);
327: PetscViewerGetFormat(viewer,&format);
328: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
329: PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);
330: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
331: SETERRQ(PETSC_ERR_SUP,"Matlab format not supported");
332: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
333: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
334: for (i=0; i<a->mbs; i++) {
335: for (j=0; j<bs; j++) {
336: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
337: for (k=a->i[i]; k<a->i[i+1]; k++) {
338: for (l=0; l<bs; l++) {
339: #if defined(PETSC_USE_COMPLEX)
340: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
341: PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l,
342: PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
343: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
344: PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l,
345: PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
346: } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
347: PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));
348: }
349: #else
350: if (a->a[bs2*k + l*bs + j] != 0.0) {
351: PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);
352: }
353: #endif
354: }
355: }
356: PetscViewerASCIIPrintf(viewer,"\n");
357: }
358: }
359: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
360: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
361: return(0);
362: } else {
363: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
364: for (i=0; i<a->mbs; i++) {
365: for (j=0; j<bs; j++) {
366: PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
367: for (k=a->i[i]; k<a->i[i+1]; k++) {
368: for (l=0; l<bs; l++) {
369: #if defined(PETSC_USE_COMPLEX)
370: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
371: PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l,
372: PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
373: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
374: PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l,
375: PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
376: } else {
377: PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));
378: }
379: #else
380: PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);
381: #endif
382: }
383: }
384: PetscViewerASCIIPrintf(viewer,"\n");
385: }
386: }
387: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
388: }
389: PetscViewerFlush(viewer);
390: return(0);
391: }
395: static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
396: {
397: Mat A = (Mat) Aa;
398: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
400: PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap.bs,bs2=a->bs2;
401: PetscMPIInt rank;
402: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
403: MatScalar *aa;
404: MPI_Comm comm;
405: PetscViewer viewer;
406:
408: /*
409: This is nasty. If this is called from an originally parallel matrix
410: then all processes call this,but only the first has the matrix so the
411: rest should return immediately.
412: */
413: PetscObjectGetComm((PetscObject)draw,&comm);
414: MPI_Comm_rank(comm,&rank);
415: if (rank) return(0);
416:
417: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
418:
419: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
420: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
421:
422: /* loop over matrix elements drawing boxes */
423: color = PETSC_DRAW_BLUE;
424: for (i=0,row=0; i<mbs; i++,row+=bs) {
425: for (j=a->i[i]; j<a->i[i+1]; j++) {
426: y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0;
427: x_l = a->j[j]*bs; x_r = x_l + 1.0;
428: aa = a->a + j*bs2;
429: for (k=0; k<bs; k++) {
430: for (l=0; l<bs; l++) {
431: if (PetscRealPart(*aa++) >= 0.) continue;
432: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
433: }
434: }
435: }
436: }
437: color = PETSC_DRAW_CYAN;
438: for (i=0,row=0; i<mbs; i++,row+=bs) {
439: for (j=a->i[i]; j<a->i[i+1]; j++) {
440: y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0;
441: x_l = a->j[j]*bs; x_r = x_l + 1.0;
442: aa = a->a + j*bs2;
443: for (k=0; k<bs; k++) {
444: for (l=0; l<bs; l++) {
445: if (PetscRealPart(*aa++) != 0.) continue;
446: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
447: }
448: }
449: }
450: }
451:
452: color = PETSC_DRAW_RED;
453: for (i=0,row=0; i<mbs; i++,row+=bs) {
454: for (j=a->i[i]; j<a->i[i+1]; j++) {
455: y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0;
456: x_l = a->j[j]*bs; x_r = x_l + 1.0;
457: aa = a->a + j*bs2;
458: for (k=0; k<bs; k++) {
459: for (l=0; l<bs; l++) {
460: if (PetscRealPart(*aa++) <= 0.) continue;
461: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
462: }
463: }
464: }
465: }
466: return(0);
467: }
471: static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
472: {
474: PetscReal xl,yl,xr,yr,w,h;
475: PetscDraw draw;
476: PetscTruth isnull;
477:
479:
480: PetscViewerDrawGetDraw(viewer,0,&draw);
481: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
482:
483: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
484: xr = A->rmap.N; yr = A->rmap.N; h = yr/10.0; w = xr/10.0;
485: xr += w; yr += h; xl = -w; yl = -h;
486: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
487: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
488: PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
489: return(0);
490: }
494: PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
495: {
497: PetscTruth iascii,isdraw;
498:
500: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
501: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
502: if (iascii){
503: MatView_SeqSBAIJ_ASCII(A,viewer);
504: } else if (isdraw) {
505: MatView_SeqSBAIJ_Draw(A,viewer);
506: } else {
507: Mat B;
508: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
509: MatView(B,viewer);
510: MatDestroy(B);
511: }
512: return(0);
513: }
518: PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
519: {
520: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
521: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
522: PetscInt *ai = a->i,*ailen = a->ilen;
523: PetscInt brow,bcol,ridx,cidx,bs=A->rmap.bs,bs2=a->bs2;
524: MatScalar *ap,*aa = a->a,zero = 0.0;
525:
527: for (k=0; k<m; k++) { /* loop over rows */
528: row = im[k]; brow = row/bs;
529: if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row);
530: if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1);
531: rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ;
532: nrow = ailen[brow];
533: for (l=0; l<n; l++) { /* loop over columns */
534: if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]);
535: if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1);
536: col = in[l] ;
537: bcol = col/bs;
538: cidx = col%bs;
539: ridx = row%bs;
540: high = nrow;
541: low = 0; /* assume unsorted */
542: while (high-low > 5) {
543: t = (low+high)/2;
544: if (rp[t] > bcol) high = t;
545: else low = t;
546: }
547: for (i=low; i<high; i++) {
548: if (rp[i] > bcol) break;
549: if (rp[i] == bcol) {
550: *v++ = ap[bs2*i+bs*cidx+ridx];
551: goto finished;
552: }
553: }
554: *v++ = zero;
555: finished:;
556: }
557: }
558: return(0);
559: }
564: PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
565: {
566: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
567: PetscErrorCode ierr;
568: PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
569: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen;
570: PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap.bs,stepval;
571: PetscTruth roworiented=a->roworiented;
572: const MatScalar *value = v;
573: MatScalar *ap,*aa = a->a,*bap;
574:
576: if (roworiented) {
577: stepval = (n-1)*bs;
578: } else {
579: stepval = (m-1)*bs;
580: }
581: for (k=0; k<m; k++) { /* loop over added rows */
582: row = im[k];
583: if (row < 0) continue;
584: #if defined(PETSC_USE_DEBUG)
585: if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1);
586: #endif
587: rp = aj + ai[row];
588: ap = aa + bs2*ai[row];
589: rmax = imax[row];
590: nrow = ailen[row];
591: low = 0;
592: high = nrow;
593: for (l=0; l<n; l++) { /* loop over added columns */
594: if (in[l] < 0) continue;
595: col = in[l];
596: #if defined(PETSC_USE_DEBUG)
597: if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1);
598: #endif
599: if (col < row) continue; /* ignore lower triangular block */
600: if (roworiented) {
601: value = v + k*(stepval+bs)*bs + l*bs;
602: } else {
603: value = v + l*(stepval+bs)*bs + k*bs;
604: }
605: if (col <= lastcol) low = 0; else high = nrow;
606: lastcol = col;
607: while (high-low > 7) {
608: t = (low+high)/2;
609: if (rp[t] > col) high = t;
610: else low = t;
611: }
612: for (i=low; i<high; i++) {
613: if (rp[i] > col) break;
614: if (rp[i] == col) {
615: bap = ap + bs2*i;
616: if (roworiented) {
617: if (is == ADD_VALUES) {
618: for (ii=0; ii<bs; ii++,value+=stepval) {
619: for (jj=ii; jj<bs2; jj+=bs) {
620: bap[jj] += *value++;
621: }
622: }
623: } else {
624: for (ii=0; ii<bs; ii++,value+=stepval) {
625: for (jj=ii; jj<bs2; jj+=bs) {
626: bap[jj] = *value++;
627: }
628: }
629: }
630: } else {
631: if (is == ADD_VALUES) {
632: for (ii=0; ii<bs; ii++,value+=stepval) {
633: for (jj=0; jj<bs; jj++) {
634: *bap++ += *value++;
635: }
636: }
637: } else {
638: for (ii=0; ii<bs; ii++,value+=stepval) {
639: for (jj=0; jj<bs; jj++) {
640: *bap++ = *value++;
641: }
642: }
643: }
644: }
645: goto noinsert2;
646: }
647: }
648: if (nonew == 1) goto noinsert2;
649: if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
650: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
651: N = nrow++ - 1; high++;
652: /* shift up all the later entries in this row */
653: for (ii=N; ii>=i; ii--) {
654: rp[ii+1] = rp[ii];
655: PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));
656: }
657: if (N >= i) {
658: PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));
659: }
660: rp[i] = col;
661: bap = ap + bs2*i;
662: if (roworiented) {
663: for (ii=0; ii<bs; ii++,value+=stepval) {
664: for (jj=ii; jj<bs2; jj+=bs) {
665: bap[jj] = *value++;
666: }
667: }
668: } else {
669: for (ii=0; ii<bs; ii++,value+=stepval) {
670: for (jj=0; jj<bs; jj++) {
671: *bap++ = *value++;
672: }
673: }
674: }
675: noinsert2:;
676: low = i;
677: }
678: ailen[row] = nrow;
679: }
680: return(0);
681: }
685: PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
686: {
687: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
689: PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
690: PetscInt m = A->rmap.N,*ip,N,*ailen = a->ilen;
691: PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
692: MatScalar *aa = a->a,*ap;
693:
695: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
696:
697: if (m) rmax = ailen[0];
698: for (i=1; i<mbs; i++) {
699: /* move each row back by the amount of empty slots (fshift) before it*/
700: fshift += imax[i-1] - ailen[i-1];
701: rmax = PetscMax(rmax,ailen[i]);
702: if (fshift) {
703: ip = aj + ai[i]; ap = aa + bs2*ai[i];
704: N = ailen[i];
705: for (j=0; j<N; j++) {
706: ip[j-fshift] = ip[j];
707: PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));
708: }
709: }
710: ai[i] = ai[i-1] + ailen[i-1];
711: }
712: if (mbs) {
713: fshift += imax[mbs-1] - ailen[mbs-1];
714: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
715: }
716: /* reset ilen and imax for each row */
717: for (i=0; i<mbs; i++) {
718: ailen[i] = imax[i] = ai[i+1] - ai[i];
719: }
720: a->nz = ai[mbs];
721:
722: /* diagonals may have moved, reset it */
723: if (a->diag) {
724: PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));
725: }
726: 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);
727: PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
728: PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);
729: a->reallocs = 0;
730: A->info.nz_unneeded = (PetscReal)fshift*bs2;
731: return(0);
732: }
734: /*
735: This function returns an array of flags which indicate the locations of contiguous
736: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
737: then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
738: Assume: sizes should be long enough to hold all the values.
739: */
742: PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
743: {
744: PetscInt i,j,k,row;
745: PetscTruth flg;
746:
748: for (i=0,j=0; i<n; j++) {
749: row = idx[i];
750: if (row%bs!=0) { /* Not the begining of a block */
751: sizes[j] = 1;
752: i++;
753: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
754: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
755: i++;
756: } else { /* Begining of the block, so check if the complete block exists */
757: flg = PETSC_TRUE;
758: for (k=1; k<bs; k++) {
759: if (row+k != idx[i+k]) { /* break in the block */
760: flg = PETSC_FALSE;
761: break;
762: }
763: }
764: if (flg) { /* No break in the bs */
765: sizes[j] = bs;
766: i+= bs;
767: } else {
768: sizes[j] = 1;
769: i++;
770: }
771: }
772: }
773: *bs_max = j;
774: return(0);
775: }
778: /* Only add/insert a(i,j) with i<=j (blocks).
779: Any a(i,j) with i>j input by user is ingored.
780: */
784: PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
785: {
786: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
788: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
789: PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
790: PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap.bs,brow,bcol;
791: PetscInt ridx,cidx,bs2=a->bs2;
792: MatScalar *ap,value,*aa=a->a,*bap;
793:
795: for (k=0; k<m; k++) { /* loop over added rows */
796: row = im[k]; /* row number */
797: brow = row/bs; /* block row number */
798: if (row < 0) continue;
799: #if defined(PETSC_USE_DEBUG)
800: if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1);
801: #endif
802: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
803: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
804: rmax = imax[brow]; /* maximum space allocated for this row */
805: nrow = ailen[brow]; /* actual length of this row */
806: low = 0;
807:
808: for (l=0; l<n; l++) { /* loop over added columns */
809: if (in[l] < 0) continue;
810: #if defined(PETSC_USE_DEBUG)
811: if (in[l] >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap.N-1);
812: #endif
813: col = in[l];
814: bcol = col/bs; /* block col number */
815:
816: if (brow > bcol) {
817: if (a->ignore_ltriangular){
818: continue; /* ignore lower triangular values */
819: } else {
820: SETERRQ(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)");
821: }
822: }
823:
824: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
825: if ((brow==bcol && ridx<=cidx) || (brow<bcol)){
826: /* element value a(k,l) */
827: if (roworiented) {
828: value = v[l + k*n];
829: } else {
830: value = v[k + l*m];
831: }
832:
833: /* move pointer bap to a(k,l) quickly and add/insert value */
834: if (col <= lastcol) low = 0; high = nrow;
835: lastcol = col;
836: while (high-low > 7) {
837: t = (low+high)/2;
838: if (rp[t] > bcol) high = t;
839: else low = t;
840: }
841: for (i=low; i<high; i++) {
842: if (rp[i] > bcol) break;
843: if (rp[i] == bcol) {
844: bap = ap + bs2*i + bs*cidx + ridx;
845: if (is == ADD_VALUES) *bap += value;
846: else *bap = value;
847: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
848: if (brow == bcol && ridx < cidx){
849: bap = ap + bs2*i + bs*ridx + cidx;
850: if (is == ADD_VALUES) *bap += value;
851: else *bap = value;
852: }
853: goto noinsert1;
854: }
855: }
856:
857: if (nonew == 1) goto noinsert1;
858: if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
859: MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
860:
861: N = nrow++ - 1; high++;
862: /* shift up all the later entries in this row */
863: for (ii=N; ii>=i; ii--) {
864: rp[ii+1] = rp[ii];
865: PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));
866: }
867: if (N>=i) {
868: PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));
869: }
870: rp[i] = bcol;
871: ap[bs2*i + bs*cidx + ridx] = value;
872: noinsert1:;
873: low = i;
874: }
875: } /* end of loop over added columns */
876: ailen[brow] = nrow;
877: } /* end of loop over added rows */
878: return(0);
879: }
883: PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,MatFactorInfo *info)
884: {
885: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
886: Mat outA;
888: PetscTruth row_identity;
889:
891: if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
892: ISIdentity(row,&row_identity);
893: if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported");
894: if (inA->rmap.bs != 1) SETERRQ1(PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap.bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */
896: outA = inA;
897: inA->factor = FACTOR_CHOLESKY;
898:
899: MatMarkDiagonal_SeqSBAIJ(inA);
900: /*
901: Blocksize < 8 have a special faster factorization/solver
902: for ICC(0) factorization with natural ordering
903: */
904: switch (inA->rmap.bs){ /* Note: row_identity = PETSC_TRUE! */
905: case 1:
906: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering;
907: inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
908: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
909: inA->ops->solves = MatSolves_SeqSBAIJ_1;
910: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering;
911: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering;
912: PetscInfo(inA,"Using special in-place natural ordering solvetrans BS=1\n");
913: break;
914: case 2:
915: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
916: inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering;
917: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering;
918: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_2_NaturalOrdering;
919: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_2_NaturalOrdering;
920: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=2\n");
921: break;
922: case 3:
923: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
924: inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering;
925: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering;
926: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_3_NaturalOrdering;
927: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_3_NaturalOrdering;
928: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=3\n");
929: break;
930: case 4:
931: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
932: inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering;
933: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering;
934: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_4_NaturalOrdering;
935: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_4_NaturalOrdering;
936: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n");
937: break;
938: case 5:
939: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
940: inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering;
941: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering;
942: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_5_NaturalOrdering;
943: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_5_NaturalOrdering;
944: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n");
945: break;
946: case 6:
947: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
948: inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering;
949: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering;
950: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_6_NaturalOrdering;
951: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_6_NaturalOrdering;
952: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=6\n");
953: break;
954: case 7:
955: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
956: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering;
957: inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering;
958: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_7_NaturalOrdering;
959: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_7_NaturalOrdering;
960: PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=7\n");
961: break;
962: default:
963: inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
964: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering;
965: inA->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering;
966: inA->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_N_NaturalOrdering;
967: inA->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_N_NaturalOrdering;
968: break;
969: }
970:
971: a->row = row;
972: a->col = row;
973: PetscObjectReference((PetscObject)row);
974: PetscObjectReference((PetscObject)row);
975:
976: /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
977: if (a->icol) {ISInvertPermutation(row,PETSC_DECIDE, &a->icol);}
978: PetscLogObjectParent(inA,a->icol);
979:
980: if (!a->solve_work) {
981: PetscMalloc((inA->rmap.N+inA->rmap.bs)*sizeof(PetscScalar),&a->solve_work);
982: PetscLogObjectMemory(inA,(inA->rmap.N+inA->rmap.bs)*sizeof(PetscScalar));
983: }
984:
985: MatCholeskyFactorNumeric(inA,info,&outA);
986: return(0);
987: }
992: PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
993: {
994: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data;
995: PetscInt i,nz,n;
996:
998: nz = baij->maxnz;
999: n = mat->cmap.n;
1000: for (i=0; i<nz; i++) {
1001: baij->j[i] = indices[i];
1002: }
1003: baij->nz = nz;
1004: for (i=0; i<n; i++) {
1005: baij->ilen[i] = baij->imax[i];
1006: }
1007: return(0);
1008: }
1013: /*@
1014: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1015: in the matrix.
1016:
1017: Input Parameters:
1018: + mat - the SeqSBAIJ matrix
1019: - indices - the column indices
1020:
1021: Level: advanced
1022:
1023: Notes:
1024: This can be called if you have precomputed the nonzero structure of the
1025: matrix and want to provide it to the matrix object to improve the performance
1026: of the MatSetValues() operation.
1027:
1028: You MUST have set the correct numbers of nonzeros per row in the call to
1029: MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1030:
1031: MUST be called before any calls to MatSetValues()
1032:
1033: .seealso: MatCreateSeqSBAIJ
1034: @*/
1035: PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1036: {
1037: PetscErrorCode ierr,(*f)(Mat,PetscInt *);
1038:
1042: PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);
1043: if (f) {
1044: (*f)(mat,indices);
1045: } else {
1046: SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices");
1047: }
1048: return(0);
1049: }
1053: PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1054: {
1058: /* If the two matrices have the same copy implementation, use fast copy. */
1059: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
1060: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1061: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1063: if (a->i[A->rmap.N] != b->i[B->rmap.N]) {
1064: SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different");
1065: }
1066: PetscMemcpy(b->a,a->a,(a->i[A->rmap.N])*sizeof(PetscScalar));
1067: } else {
1068: MatGetRowUpperTriangular(A);
1069: MatCopy_Basic(A,B,str);
1070: MatRestoreRowUpperTriangular(A);
1071: }
1072: return(0);
1073: }
1077: PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A)
1078: {
1080:
1082: MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,PetscMax(A->rmap.bs,1),PETSC_DEFAULT,0);
1083: return(0);
1084: }
1088: PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1089: {
1090: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1092: *array = a->a;
1093: return(0);
1094: }
1098: PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1099: {
1101: return(0);
1102: }
1104: #include petscblaslapack.h
1107: PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1108: {
1109: Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data;
1111: PetscInt i,bs=Y->rmap.bs,bs2,j;
1112: PetscBLASInt bnz = (PetscBLASInt)x->nz,one = 1;
1113:
1115: if (str == SAME_NONZERO_PATTERN) {
1116: PetscScalar alpha = a;
1117: BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one);
1118: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1119: if (y->xtoy && y->XtoY != X) {
1120: PetscFree(y->xtoy);
1121: MatDestroy(y->XtoY);
1122: }
1123: if (!y->xtoy) { /* get xtoy */
1124: MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);
1125: y->XtoY = X;
1126: }
1127: bs2 = bs*bs;
1128: for (i=0; i<x->nz; i++) {
1129: j = 0;
1130: while (j < bs2){
1131: y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]);
1132: j++;
1133: }
1134: }
1135: PetscInfo3(0,"ratio of nnz_s(X)/nnz_s(Y): %D/%D = %G\n",bs2*x->nz,bs2*y->nz,(PetscReal)(bs2*x->nz)/(bs2*y->nz));
1136: } else {
1137: MatGetRowUpperTriangular(X);
1138: MatAXPY_Basic(Y,a,X,str);
1139: MatRestoreRowUpperTriangular(X);
1140: }
1141: return(0);
1142: }
1146: PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg)
1147: {
1149: *flg = PETSC_TRUE;
1150: return(0);
1151: }
1155: PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg)
1156: {
1158: *flg = PETSC_TRUE;
1159: return(0);
1160: }
1164: PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscTruth *flg)
1165: {
1167: *flg = PETSC_FALSE;
1168: return(0);
1169: }
1173: PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1174: {
1175: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1176: PetscInt i,nz = a->bs2*a->i[a->mbs];
1177: PetscScalar *aa = a->a;
1180: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1181: return(0);
1182: }
1186: PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1187: {
1188: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1189: PetscInt i,nz = a->bs2*a->i[a->mbs];
1190: PetscScalar *aa = a->a;
1193: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1194: return(0);
1195: }
1197: /* -------------------------------------------------------------------*/
1198: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1199: MatGetRow_SeqSBAIJ,
1200: MatRestoreRow_SeqSBAIJ,
1201: MatMult_SeqSBAIJ_N,
1202: /* 4*/ MatMultAdd_SeqSBAIJ_N,
1203: MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1204: MatMultAdd_SeqSBAIJ_N,
1205: MatSolve_SeqSBAIJ_N,
1206: 0,
1207: 0,
1208: /*10*/ 0,
1209: 0,
1210: MatCholeskyFactor_SeqSBAIJ,
1211: MatRelax_SeqSBAIJ,
1212: MatTranspose_SeqSBAIJ,
1213: /*15*/ MatGetInfo_SeqSBAIJ,
1214: MatEqual_SeqSBAIJ,
1215: MatGetDiagonal_SeqSBAIJ,
1216: MatDiagonalScale_SeqSBAIJ,
1217: MatNorm_SeqSBAIJ,
1218: /*20*/ 0,
1219: MatAssemblyEnd_SeqSBAIJ,
1220: 0,
1221: MatSetOption_SeqSBAIJ,
1222: MatZeroEntries_SeqSBAIJ,
1223: /*25*/ 0,
1224: 0,
1225: 0,
1226: MatCholeskyFactorSymbolic_SeqSBAIJ,
1227: MatCholeskyFactorNumeric_SeqSBAIJ_N,
1228: /*30*/ MatSetUpPreallocation_SeqSBAIJ,
1229: 0,
1230: MatICCFactorSymbolic_SeqSBAIJ,
1231: MatGetArray_SeqSBAIJ,
1232: MatRestoreArray_SeqSBAIJ,
1233: /*35*/ MatDuplicate_SeqSBAIJ,
1234: MatForwardSolve_SeqSBAIJ_N,
1235: MatBackwardSolve_SeqSBAIJ_N,
1236: 0,
1237: MatICCFactor_SeqSBAIJ,
1238: /*40*/ MatAXPY_SeqSBAIJ,
1239: MatGetSubMatrices_SeqSBAIJ,
1240: MatIncreaseOverlap_SeqSBAIJ,
1241: MatGetValues_SeqSBAIJ,
1242: MatCopy_SeqSBAIJ,
1243: /*45*/ 0,
1244: MatScale_SeqSBAIJ,
1245: 0,
1246: 0,
1247: 0,
1248: /*50*/ 0,
1249: MatGetRowIJ_SeqSBAIJ,
1250: MatRestoreRowIJ_SeqSBAIJ,
1251: 0,
1252: 0,
1253: /*55*/ 0,
1254: 0,
1255: 0,
1256: 0,
1257: MatSetValuesBlocked_SeqSBAIJ,
1258: /*60*/ MatGetSubMatrix_SeqSBAIJ,
1259: 0,
1260: 0,
1261: 0,
1262: 0,
1263: /*65*/ 0,
1264: 0,
1265: 0,
1266: 0,
1267: 0,
1268: /*70*/ MatGetRowMax_SeqSBAIJ,
1269: 0,
1270: 0,
1271: 0,
1272: 0,
1273: /*75*/ 0,
1274: 0,
1275: 0,
1276: 0,
1277: 0,
1278: /*80*/ 0,
1279: 0,
1280: 0,
1281: #if !defined(PETSC_USE_COMPLEX)
1282: MatGetInertia_SeqSBAIJ,
1283: #else
1284: 0,
1285: #endif
1286: MatLoad_SeqSBAIJ,
1287: /*85*/ MatIsSymmetric_SeqSBAIJ,
1288: MatIsHermitian_SeqSBAIJ,
1289: MatIsStructurallySymmetric_SeqSBAIJ,
1290: 0,
1291: 0,
1292: /*90*/ 0,
1293: 0,
1294: 0,
1295: 0,
1296: 0,
1297: /*95*/ 0,
1298: 0,
1299: 0,
1300: 0,
1301: 0,
1302: /*100*/0,
1303: 0,
1304: 0,
1305: 0,
1306: 0,
1307: /*105*/0,
1308: MatRealPart_SeqSBAIJ,
1309: MatImaginaryPart_SeqSBAIJ,
1310: MatGetRowUpperTriangular_SeqSBAIJ,
1311: MatRestoreRowUpperTriangular_SeqSBAIJ
1312: };
1317: PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1318: {
1319: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1320: PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2;
1322:
1324: if (aij->nonew != 1) {
1325: SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first");
1326: }
1327:
1328: /* allocate space for values if not already there */
1329: if (!aij->saved_values) {
1330: PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);
1331: }
1332:
1333: /* copy values over */
1334: PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));
1335: return(0);
1336: }
1342: PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1343: {
1344: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1346: PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2;
1347:
1349: if (aij->nonew != 1) {
1350: SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first");
1351: }
1352: if (!aij->saved_values) {
1353: SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1354: }
1355:
1356: /* copy values over */
1357: PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));
1358: return(0);
1359: }
1365: PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1366: {
1367: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1369: PetscInt i,mbs,bs2;
1370: PetscTruth skipallocation = PETSC_FALSE,flg;
1371:
1373: B->preallocated = PETSC_TRUE;
1374: PetscOptionsGetInt(B->prefix,"-mat_block_size",&bs,PETSC_NULL);
1375: B->rmap.bs = B->cmap.bs = bs;
1376: PetscMapInitialize(B->comm,&B->rmap);
1377: PetscMapInitialize(B->comm,&B->cmap);
1379: mbs = B->rmap.N/bs;
1380: bs2 = bs*bs;
1381:
1382: if (mbs*bs != B->rmap.N) {
1383: SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize");
1384: }
1385:
1386: if (nz == MAT_SKIP_ALLOCATION) {
1387: skipallocation = PETSC_TRUE;
1388: nz = 0;
1389: }
1391: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1392: if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1393: if (nnz) {
1394: for (i=0; i<mbs; i++) {
1395: if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]);
1396: if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],mbs);
1397: }
1398: }
1399:
1400: PetscOptionsHasName(B->prefix,"-mat_no_unroll",&flg);
1401: if (!flg) {
1402: switch (bs) {
1403: case 1:
1404: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1;
1405: B->ops->solve = MatSolve_SeqSBAIJ_1;
1406: B->ops->solves = MatSolves_SeqSBAIJ_1;
1407: B->ops->solvetranspose = MatSolve_SeqSBAIJ_1;
1408: B->ops->mult = MatMult_SeqSBAIJ_1;
1409: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1410: B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1411: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1412: break;
1413: case 2:
1414: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1415: B->ops->solve = MatSolve_SeqSBAIJ_2;
1416: B->ops->solvetranspose = MatSolve_SeqSBAIJ_2;
1417: B->ops->mult = MatMult_SeqSBAIJ_2;
1418: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1419: B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1420: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1421: break;
1422: case 3:
1423: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1424: B->ops->solve = MatSolve_SeqSBAIJ_3;
1425: B->ops->solvetranspose = MatSolve_SeqSBAIJ_3;
1426: B->ops->mult = MatMult_SeqSBAIJ_3;
1427: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1428: B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1429: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1430: break;
1431: case 4:
1432: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1433: B->ops->solve = MatSolve_SeqSBAIJ_4;
1434: B->ops->solvetranspose = MatSolve_SeqSBAIJ_4;
1435: B->ops->mult = MatMult_SeqSBAIJ_4;
1436: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1437: B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1438: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1439: break;
1440: case 5:
1441: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1442: B->ops->solve = MatSolve_SeqSBAIJ_5;
1443: B->ops->solvetranspose = MatSolve_SeqSBAIJ_5;
1444: B->ops->mult = MatMult_SeqSBAIJ_5;
1445: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1446: B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1447: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1448: break;
1449: case 6:
1450: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1451: B->ops->solve = MatSolve_SeqSBAIJ_6;
1452: B->ops->solvetranspose = MatSolve_SeqSBAIJ_6;
1453: B->ops->mult = MatMult_SeqSBAIJ_6;
1454: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1455: B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1456: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1457: break;
1458: case 7:
1459: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1460: B->ops->solve = MatSolve_SeqSBAIJ_7;
1461: B->ops->solvetranspose = MatSolve_SeqSBAIJ_7;
1462: B->ops->mult = MatMult_SeqSBAIJ_7;
1463: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1464: B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1465: B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1466: break;
1467: }
1468: }
1469:
1470: b->mbs = mbs;
1471: b->nbs = mbs;
1472: if (!skipallocation) {
1473: /* b->ilen will count nonzeros in each block row so far. */
1474: PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);
1475: for (i=0; i<mbs; i++) { b->ilen[i] = 0;}
1476: if (!nnz) {
1477: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1478: else if (nz <= 0) nz = 1;
1479: for (i=0; i<mbs; i++) {
1480: b->imax[i] = nz;
1481: }
1482: nz = nz*mbs; /* total nz */
1483: } else {
1484: nz = 0;
1485: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];}
1486: }
1487: /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1488:
1489: /* allocate the matrix space */
1490: PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.N+1,PetscInt,&b->i);
1491: PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));
1492: PetscMemzero(b->j,nz*sizeof(PetscInt));
1493: b->singlemalloc = PETSC_TRUE;
1494:
1495: /* pointer to beginning of each row */
1496: b->i[0] = 0;
1497: for (i=1; i<mbs+1; i++) {
1498: b->i[i] = b->i[i-1] + b->imax[i-1];
1499: }
1500: b->free_a = PETSC_TRUE;
1501: b->free_ij = PETSC_TRUE;
1502: } else {
1503: b->free_a = PETSC_FALSE;
1504: b->free_ij = PETSC_FALSE;
1505: }
1506:
1507: B->rmap.bs = bs;
1508: b->bs2 = bs2;
1509: b->nz = 0;
1510: b->maxnz = nz*bs2;
1511:
1512: b->inew = 0;
1513: b->jnew = 0;
1514: b->anew = 0;
1515: b->a2anew = 0;
1516: b->permute = PETSC_FALSE;
1517: return(0);
1518: }
1522: EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
1523: EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*);
1526: /*MC
1527: MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1528: based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1529:
1530: Options Database Keys:
1531: . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1532:
1533: Level: beginner
1534:
1535: .seealso: MatCreateSeqSBAIJ
1536: M*/
1541: PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1542: {
1543: Mat_SeqSBAIJ *b;
1545: PetscMPIInt size;
1546: PetscTruth flg;
1547:
1549: MPI_Comm_size(B->comm,&size);
1550: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1551:
1552: PetscNew(Mat_SeqSBAIJ,&b);
1553: B->data = (void*)b;
1554: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1555: B->ops->destroy = MatDestroy_SeqSBAIJ;
1556: B->ops->view = MatView_SeqSBAIJ;
1557: B->factor = 0;
1558: B->mapping = 0;
1559: b->row = 0;
1560: b->icol = 0;
1561: b->reallocs = 0;
1562: b->saved_values = 0;
1563:
1564:
1565: b->sorted = PETSC_FALSE;
1566: b->roworiented = PETSC_TRUE;
1567: b->nonew = 0;
1568: b->diag = 0;
1569: b->solve_work = 0;
1570: b->mult_work = 0;
1571: B->spptr = 0;
1572: b->keepzeroedrows = PETSC_FALSE;
1573: b->xtoy = 0;
1574: b->XtoY = 0;
1575:
1576: b->inew = 0;
1577: b->jnew = 0;
1578: b->anew = 0;
1579: b->a2anew = 0;
1580: b->permute = PETSC_FALSE;
1582: b->ignore_ltriangular = PETSC_FALSE;
1583: PetscOptionsHasName(PETSC_NULL,"-mat_ignore_lower_triangular",&flg);
1584: if (flg) b->ignore_ltriangular = PETSC_TRUE;
1586: b->getrow_utriangular = PETSC_FALSE;
1587: PetscOptionsHasName(PETSC_NULL,"-mat_getrow_uppertriangular",&flg);
1588: if (flg) b->getrow_utriangular = PETSC_TRUE;
1590: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C",
1591: "MatStoreValues_SeqSBAIJ",
1592: MatStoreValues_SeqSBAIJ);
1593: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C",
1594: "MatRetrieveValues_SeqSBAIJ",
1595: (void*)MatRetrieveValues_SeqSBAIJ);
1596: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",
1597: "MatSeqSBAIJSetColumnIndices_SeqSBAIJ",
1598: MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1599: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",
1600: "MatConvert_SeqSBAIJ_SeqAIJ",
1601: MatConvert_SeqSBAIJ_SeqAIJ);
1602: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",
1603: "MatConvert_SeqSBAIJ_SeqBAIJ",
1604: MatConvert_SeqSBAIJ_SeqBAIJ);
1605: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",
1606: "MatSeqSBAIJSetPreallocation_SeqSBAIJ",
1607: MatSeqSBAIJSetPreallocation_SeqSBAIJ);
1609: B->symmetric = PETSC_TRUE;
1610: B->structurally_symmetric = PETSC_TRUE;
1611: B->symmetric_set = PETSC_TRUE;
1612: B->structurally_symmetric_set = PETSC_TRUE;
1613: PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);
1614: return(0);
1615: }
1620: /*@C
1621: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
1622: compressed row) format. For good matrix assembly performance the
1623: user should preallocate the matrix storage by setting the parameter nz
1624: (or the array nnz). By setting these parameters accurately, performance
1625: during matrix assembly can be increased by more than a factor of 50.
1627: Collective on Mat
1629: Input Parameters:
1630: + A - the symmetric matrix
1631: . bs - size of block
1632: . nz - number of block nonzeros per block row (same for all rows)
1633: - nnz - array containing the number of block nonzeros in the upper triangular plus
1634: diagonal portion of each block (possibly different for each block row) or PETSC_NULL
1636: Options Database Keys:
1637: . -mat_no_unroll - uses code that does not unroll the loops in the
1638: block calculations (much slower)
1639: . -mat_block_size - size of the blocks to use
1641: Level: intermediate
1643: Notes:
1644: Specify the preallocated storage with either nz or nnz (not both).
1645: Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory
1646: allocation. For additional details, see the users manual chapter on
1647: matrices.
1649: If the nnz parameter is given then the nz parameter is ignored
1652: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ()
1653: @*/
1654: PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
1655: {
1656: PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]);
1659: PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);
1660: if (f) {
1661: (*f)(B,bs,nz,nnz);
1662: }
1663: return(0);
1664: }
1668: /*@C
1669: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
1670: compressed row) format. For good matrix assembly performance the
1671: user should preallocate the matrix storage by setting the parameter nz
1672: (or the array nnz). By setting these parameters accurately, performance
1673: during matrix assembly can be increased by more than a factor of 50.
1675: Collective on MPI_Comm
1677: Input Parameters:
1678: + comm - MPI communicator, set to PETSC_COMM_SELF
1679: . bs - size of block
1680: . m - number of rows, or number of columns
1681: . nz - number of block nonzeros per block row (same for all rows)
1682: - nnz - array containing the number of block nonzeros in the upper triangular plus
1683: diagonal portion of each block (possibly different for each block row) or PETSC_NULL
1685: Output Parameter:
1686: . A - the symmetric matrix
1688: Options Database Keys:
1689: . -mat_no_unroll - uses code that does not unroll the loops in the
1690: block calculations (much slower)
1691: . -mat_block_size - size of the blocks to use
1693: Level: intermediate
1695: Notes:
1697: Specify the preallocated storage with either nz or nnz (not both).
1698: Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory
1699: allocation. For additional details, see the users manual chapter on
1700: matrices.
1702: If the nnz parameter is given then the nz parameter is ignored
1704: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ()
1705: @*/
1706: PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
1707: {
1709:
1711: MatCreate(comm,A);
1712: MatSetSizes(*A,m,n,m,n);
1713: MatSetType(*A,MATSEQSBAIJ);
1714: MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);
1715: return(0);
1716: }
1720: PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
1721: {
1722: Mat C;
1723: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
1725: PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
1728: if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix");
1730: *B = 0;
1731: MatCreate(A->comm,&C);
1732: MatSetSizes(C,A->rmap.N,A->cmap.n,A->rmap.N,A->cmap.n);
1733: MatSetType(C,A->type_name);
1734: PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));
1735: c = (Mat_SeqSBAIJ*)C->data;
1737: C->preallocated = PETSC_TRUE;
1738: C->factor = A->factor;
1739: c->row = 0;
1740: c->icol = 0;
1741: c->saved_values = 0;
1742: c->keepzeroedrows = a->keepzeroedrows;
1743: C->assembled = PETSC_TRUE;
1745: PetscMapCopy(A->comm,&A->rmap,&C->rmap);
1746: PetscMapCopy(A->comm,&A->cmap,&C->cmap);
1747: c->bs2 = a->bs2;
1748: c->mbs = a->mbs;
1749: c->nbs = a->nbs;
1751: PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);
1752: for (i=0; i<mbs; i++) {
1753: c->imax[i] = a->imax[i];
1754: c->ilen[i] = a->ilen[i];
1755: }
1757: /* allocate the matrix space */
1758: PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);
1759: c->singlemalloc = PETSC_TRUE;
1760: PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));
1761: PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));
1762: if (mbs > 0) {
1763: PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));
1764: if (cpvalues == MAT_COPY_VALUES) {
1765: PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));
1766: } else {
1767: PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));
1768: }
1769: }
1771: c->sorted = a->sorted;
1772: c->roworiented = a->roworiented;
1773: c->nonew = a->nonew;
1775: if (a->diag) {
1776: PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);
1777: PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));
1778: for (i=0; i<mbs; i++) {
1779: c->diag[i] = a->diag[i];
1780: }
1781: } else c->diag = 0;
1782: c->nz = a->nz;
1783: c->maxnz = a->maxnz;
1784: c->solve_work = 0;
1785: c->mult_work = 0;
1786: c->free_a = PETSC_TRUE;
1787: c->free_ij = PETSC_TRUE;
1788: *B = C;
1789: PetscFListDuplicate(A->qlist,&C->qlist);
1790: return(0);
1791: }
1795: PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, MatType type,Mat *A)
1796: {
1797: Mat_SeqSBAIJ *a;
1798: Mat B;
1800: int fd;
1801: PetscMPIInt size;
1802: PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1;
1803: PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount;
1804: PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows;
1805: PetscInt *masked,nmask,tmp,bs2,ishift;
1806: PetscScalar *aa;
1807: MPI_Comm comm = ((PetscObject)viewer)->comm;
1810: PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);
1811: bs2 = bs*bs;
1813: MPI_Comm_size(comm,&size);
1814: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor");
1815: PetscViewerBinaryGetDescriptor(viewer,&fd);
1816: PetscBinaryRead(fd,header,4,PETSC_INT);
1817: if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object");
1818: M = header[1]; N = header[2]; nz = header[3];
1820: if (header[3] < 0) {
1821: SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ");
1822: }
1824: if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices");
1826: /*
1827: This code adds extra rows to make sure the number of rows is
1828: divisible by the blocksize
1829: */
1830: mbs = M/bs;
1831: extra_rows = bs - M + bs*(mbs);
1832: if (extra_rows == bs) extra_rows = 0;
1833: else mbs++;
1834: if (extra_rows) {
1835: PetscInfo(0,"Padding loaded matrix to match blocksize\n");
1836: }
1838: /* read in row lengths */
1839: PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);
1840: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1841: for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;
1843: /* read in column indices */
1844: PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);
1845: PetscBinaryRead(fd,jj,nz,PETSC_INT);
1846: for (i=0; i<extra_rows; i++) jj[nz+i] = M+i;
1848: /* loop over row lengths determining block row lengths */
1849: PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);
1850: PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));
1851: PetscMalloc(2*mbs*sizeof(PetscInt),&mask);
1852: PetscMemzero(mask,mbs*sizeof(PetscInt));
1853: masked = mask + mbs;
1854: rowcount = 0; nzcount = 0;
1855: for (i=0; i<mbs; i++) {
1856: nmask = 0;
1857: for (j=0; j<bs; j++) {
1858: kmax = rowlengths[rowcount];
1859: for (k=0; k<kmax; k++) {
1860: tmp = jj[nzcount++]/bs; /* block col. index */
1861: if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;}
1862: }
1863: rowcount++;
1864: }
1865: s_browlengths[i] += nmask;
1866:
1867: /* zero out the mask elements we set */
1868: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
1869: }
1871: /* create our matrix */
1872: MatCreate(comm,&B);
1873: MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);
1874: MatSetType(B,type);
1875: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);
1876: a = (Mat_SeqSBAIJ*)B->data;
1878: /* set matrix "i" values */
1879: a->i[0] = 0;
1880: for (i=1; i<= mbs; i++) {
1881: a->i[i] = a->i[i-1] + s_browlengths[i-1];
1882: a->ilen[i-1] = s_browlengths[i-1];
1883: }
1884: a->nz = a->i[mbs];
1886: /* read in nonzero values */
1887: PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);
1888: PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);
1889: for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0;
1891: /* set "a" and "j" values into matrix */
1892: nzcount = 0; jcount = 0;
1893: for (i=0; i<mbs; i++) {
1894: nzcountb = nzcount;
1895: nmask = 0;
1896: for (j=0; j<bs; j++) {
1897: kmax = rowlengths[i*bs+j];
1898: for (k=0; k<kmax; k++) {
1899: tmp = jj[nzcount++]/bs; /* block col. index */
1900: if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;}
1901: }
1902: }
1903: /* sort the masked values */
1904: PetscSortInt(nmask,masked);
1906: /* set "j" values into matrix */
1907: maskcount = 1;
1908: for (j=0; j<nmask; j++) {
1909: a->j[jcount++] = masked[j];
1910: mask[masked[j]] = maskcount++;
1911: }
1913: /* set "a" values into matrix */
1914: ishift = bs2*a->i[i];
1915: for (j=0; j<bs; j++) {
1916: kmax = rowlengths[i*bs+j];
1917: for (k=0; k<kmax; k++) {
1918: tmp = jj[nzcountb]/bs ; /* block col. index */
1919: if (tmp >= i){
1920: block = mask[tmp] - 1;
1921: point = jj[nzcountb] - bs*tmp;
1922: idx = ishift + bs2*block + j + bs*point;
1923: a->a[idx] = aa[nzcountb];
1924: }
1925: nzcountb++;
1926: }
1927: }
1928: /* zero out the mask elements we set */
1929: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
1930: }
1931: if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix");
1933: PetscFree(rowlengths);
1934: PetscFree(s_browlengths);
1935: PetscFree(aa);
1936: PetscFree(jj);
1937: PetscFree(mask);
1939: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1940: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1941: MatView_Private(B);
1942: *A = B;
1943: return(0);
1944: }
1948: PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1949: {
1950: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1951: MatScalar *aa=a->a,*v,*v1;
1952: PetscScalar *x,*b,*t,sum,d;
1954: PetscInt m=a->mbs,bs=A->rmap.bs,*ai=a->i,*aj=a->j;
1955: PetscInt nz,nz1,*vj,*vj1,i;
1958: if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat");
1960: its = its*lits;
1961: if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
1963: if (bs > 1)
1964: SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented");
1966: VecGetArray(xx,&x);
1967: if (xx != bb) {
1968: VecGetArray(bb,&b);
1969: } else {
1970: b = x;
1971: }
1973: if (!a->relax_work) {
1974: PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);
1975: }
1976: t = a->relax_work;
1977:
1978: if (flag & SOR_ZERO_INITIAL_GUESS) {
1979: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
1980: PetscMemcpy(t,b,m*sizeof(PetscScalar));
1982: for (i=0; i<m; i++){
1983: d = *(aa + ai[i]); /* diag[i] */
1984: v = aa + ai[i] + 1;
1985: vj = aj + ai[i] + 1;
1986: nz = ai[i+1] - ai[i] - 1;
1987: PetscLogFlops(2*nz-1);
1988: x[i] = omega*t[i]/d;
1989: while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */
1990: }
1991: }
1993: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
1994: if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){
1995: t = b;
1996: }
1997:
1998: for (i=m-1; i>=0; i--){
1999: d = *(aa + ai[i]);
2000: v = aa + ai[i] + 1;
2001: vj = aj + ai[i] + 1;
2002: nz = ai[i+1] - ai[i] - 1;
2003: PetscLogFlops(2*nz-1);
2004: sum = t[i];
2005: while (nz--) sum -= x[*vj++]*(*v++);
2006: x[i] = (1-omega)*x[i] + omega*sum/d;
2007: }
2008: t = a->relax_work;
2009: }
2010: its--;
2011: }
2013: while (its--) {
2014: /*
2015: forward sweep:
2016: for i=0,...,m-1:
2017: sum[i] = (b[i] - U(i,:)x )/d[i];
2018: x[i] = (1-omega)x[i] + omega*sum[i];
2019: b = b - x[i]*U^T(i,:);
2020:
2021: */
2022: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
2023: PetscMemcpy(t,b,m*sizeof(PetscScalar));
2025: for (i=0; i<m; i++){
2026: d = *(aa + ai[i]); /* diag[i] */
2027: v = aa + ai[i] + 1; v1=v;
2028: vj = aj + ai[i] + 1; vj1=vj;
2029: nz = ai[i+1] - ai[i] - 1; nz1=nz;
2030: sum = t[i];
2031: PetscLogFlops(4*nz-2);
2032: while (nz1--) sum -= (*v1++)*x[*vj1++];
2033: x[i] = (1-omega)*x[i] + omega*sum/d;
2034: while (nz--) t[*vj++] -= x[i]*(*v++);
2035: }
2036: }
2037:
2038: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
2039: /*
2040: backward sweep:
2041: b = b - x[i]*U^T(i,:), i=0,...,n-2
2042: for i=m-1,...,0:
2043: sum[i] = (b[i] - U(i,:)x )/d[i];
2044: x[i] = (1-omega)x[i] + omega*sum[i];
2045: */
2046: /* if there was a forward sweep done above then I thing the next two for loops are not needed */
2047: PetscMemcpy(t,b,m*sizeof(PetscScalar));
2048:
2049: for (i=0; i<m-1; i++){ /* update rhs */
2050: v = aa + ai[i] + 1;
2051: vj = aj + ai[i] + 1;
2052: nz = ai[i+1] - ai[i] - 1;
2053: PetscLogFlops(2*nz-1);
2054: while (nz--) t[*vj++] -= x[i]*(*v++);
2055: }
2056: for (i=m-1; i>=0; i--){
2057: d = *(aa + ai[i]);
2058: v = aa + ai[i] + 1;
2059: vj = aj + ai[i] + 1;
2060: nz = ai[i+1] - ai[i] - 1;
2061: PetscLogFlops(2*nz-1);
2062: sum = t[i];
2063: while (nz--) sum -= x[*vj++]*(*v++);
2064: x[i] = (1-omega)*x[i] + omega*sum/d;
2065: }
2066: }
2067: }
2069: VecRestoreArray(xx,&x);
2070: if (bb != xx) {
2071: VecRestoreArray(bb,&b);
2072: }
2073: return(0);
2074: }
2078: /*@
2079: MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2080: (upper triangular entries in CSR format) provided by the user.
2082: Collective on MPI_Comm
2084: Input Parameters:
2085: + comm - must be an MPI communicator of size 1
2086: . bs - size of block
2087: . m - number of rows
2088: . n - number of columns
2089: . i - row indices
2090: . j - column indices
2091: - a - matrix values
2093: Output Parameter:
2094: . mat - the matrix
2096: Level: intermediate
2098: Notes:
2099: The i, j, and a arrays are not copied by this routine, the user must free these arrays
2100: once the matrix is destroyed
2102: You cannot set new nonzero locations into this matrix, that will generate an error.
2104: The i and j indices are 0 based
2106: .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ()
2108: @*/
2109: PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat)
2110: {
2112: PetscInt ii;
2113: Mat_SeqSBAIJ *sbaij;
2116: if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2117: if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2118:
2119: MatCreate(comm,mat);
2120: MatSetSizes(*mat,m,n,m,n);
2121: MatSetType(*mat,MATSEQSBAIJ);
2122: MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);
2123: sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2124: PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);
2126: sbaij->i = i;
2127: sbaij->j = j;
2128: sbaij->a = a;
2129: sbaij->singlemalloc = PETSC_FALSE;
2130: sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2131: sbaij->free_a = PETSC_FALSE;
2132: sbaij->free_ij = PETSC_FALSE;
2134: for (ii=0; ii<m; ii++) {
2135: sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2136: #if defined(PETSC_USE_DEBUG)
2137: if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]);
2138: #endif
2139: }
2140: #if defined(PETSC_USE_DEBUG)
2141: for (ii=0; ii<sbaij->i[m]; ii++) {
2142: if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2143: if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]);
2144: }
2145: #endif
2147: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
2148: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
2149: return(0);
2150: }