Actual source code: sbaij.c
1: /*$Id: sbaij.c,v 1.62 2001/08/07 03:03:01 balay Exp $*/
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/vec/vecimpl.h
9: #include src/inline/spops.h
10: #include src/mat/impls/sbaij/seq/sbaij.h
12: #define CHUNKSIZE 10
14: /*
15: Checks for missing diagonals
16: */
17: #undef __FUNCT__
19: int MatMissingDiagonal_SeqSBAIJ(Mat A)
20: {
21: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
22: int *diag,*jj = a->j,i,ierr;
25: MatMarkDiagonal_SeqSBAIJ(A);
26: diag = a->diag;
27: for (i=0; i<a->mbs; i++) {
28: if (jj[diag[i]] != i) SETERRQ1(1,"Matrix is missing diagonal number %d",i);
29: }
30: return(0);
31: }
33: #undef __FUNCT__
35: int MatMarkDiagonal_SeqSBAIJ(Mat A)
36: {
37: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
38: int i,mbs = a->mbs,ierr;
41: if (a->diag) return(0);
43: PetscMalloc((mbs+1)*sizeof(int),&a->diag);
44: PetscLogObjectMemory(A,(mbs+1)*sizeof(int));
45: for (i=0; i<mbs; i++) a->diag[i] = a->i[i];
46: return(0);
47: }
49: #undef __FUNCT__
51: static int MatGetRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done)
52: {
53: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
54: int 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: int s_nz = a->i[n];
63: for (i=0; i<s_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: }
72: #undef __FUNCT__
74: static int MatRestoreRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done)
75: {
76: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
77: int i,n = a->mbs;
80: if (!ia) return(0);
82: if (oshift == 1) {
83: int s_nz = a->i[n]-1;
84: for (i=0; i<s_nz; i++) a->j[i]--;
85: for (i=0; i<n+1; i++) a->i[i]--;
86: }
87: return(0);
88: }
90: #undef __FUNCT__
92: int MatGetBlockSize_SeqSBAIJ(Mat mat,int *bs)
93: {
94: Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)mat->data;
97: *bs = sbaij->bs;
98: return(0);
99: }
101: #undef __FUNCT__
103: int MatDestroy_SeqSBAIJ(Mat A)
104: {
105: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
106: int ierr;
109: #if defined(PETSC_USE_LOG)
110: PetscLogObjectState((PetscObject)A,"Rows=%d, s_NZ=%d",A->m,a->s_nz);
111: #endif
112: PetscFree(a->a);
113: if (!a->singlemalloc) {
114: PetscFree(a->i);
115: PetscFree(a->j);
116: }
117: if (a->row) {
118: ISDestroy(a->row);
119: }
120: if (a->diag) {PetscFree(a->diag);}
121: if (a->ilen) {PetscFree(a->ilen);}
122: if (a->imax) {PetscFree(a->imax);}
123: if (a->icol) {ISDestroy(a->icol);}
124: if (a->solve_work) {PetscFree(a->solve_work);}
125: if (a->solves_work) {PetscFree(a->solves_work);}
126: if (a->mult_work) {PetscFree(a->mult_work);}
127: if (a->saved_values) {PetscFree(a->saved_values);}
129: if (a->inew){
130: PetscFree(a->inew);
131: a->inew = 0;
132: }
133: PetscFree(a);
134: return(0);
135: }
137: #undef __FUNCT__
139: int MatSetOption_SeqSBAIJ(Mat A,MatOption op)
140: {
141: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
144: switch (op) {
145: case MAT_ROW_ORIENTED:
146: a->roworiented = PETSC_TRUE;
147: break;
148: case MAT_COLUMN_ORIENTED:
149: a->roworiented = PETSC_FALSE;
150: break;
151: case MAT_COLUMNS_SORTED:
152: a->sorted = PETSC_TRUE;
153: break;
154: case MAT_COLUMNS_UNSORTED:
155: a->sorted = PETSC_FALSE;
156: break;
157: case MAT_KEEP_ZEROED_ROWS:
158: a->keepzeroedrows = PETSC_TRUE;
159: break;
160: case MAT_NO_NEW_NONZERO_LOCATIONS:
161: a->nonew = 1;
162: break;
163: case MAT_NEW_NONZERO_LOCATION_ERR:
164: a->nonew = -1;
165: break;
166: case MAT_NEW_NONZERO_ALLOCATION_ERR:
167: a->nonew = -2;
168: break;
169: case MAT_YES_NEW_NONZERO_LOCATIONS:
170: a->nonew = 0;
171: break;
172: case MAT_ROWS_SORTED:
173: case MAT_ROWS_UNSORTED:
174: case MAT_YES_NEW_DIAGONALS:
175: case MAT_IGNORE_OFF_PROC_ENTRIES:
176: case MAT_USE_HASH_TABLE:
177: PetscLogInfo(A,"MatSetOption_SeqSBAIJ:Option ignoredn");
178: break;
179: case MAT_NO_NEW_DIAGONALS:
180: SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS");
181: default:
182: SETERRQ(PETSC_ERR_SUP,"unknown option");
183: }
184: return(0);
185: }
187: #undef __FUNCT__
189: int MatGetRow_SeqSBAIJ(Mat A,int row,int *ncols,int **cols,PetscScalar **v)
190: {
191: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
192: int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2,ierr;
193: MatScalar *aa,*aa_i;
194: PetscScalar *v_i;
197: bs = a->bs;
198: ai = a->i;
199: aj = a->j;
200: aa = a->a;
201: bs2 = a->bs2;
202:
203: if (row < 0 || row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range");
204:
205: bn = row/bs; /* Block number */
206: bp = row % bs; /* Block position */
207: M = ai[bn+1] - ai[bn];
208: *ncols = bs*M;
209:
210: if (v) {
211: *v = 0;
212: if (*ncols) {
213: PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);
214: for (i=0; i<M; i++) { /* for each block in the block row */
215: v_i = *v + i*bs;
216: aa_i = aa + bs2*(ai[bn] + i);
217: for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];}
218: }
219: }
220: }
222: if (cols) {
223: *cols = 0;
224: if (*ncols) {
225: PetscMalloc((*ncols+row)*sizeof(int),cols);
226: for (i=0; i<M; i++) { /* for each block in the block row */
227: cols_i = *cols + i*bs;
228: itmp = bs*aj[ai[bn] + i];
229: for (j=0; j<bs; j++) {cols_i[j] = itmp++;}
230: }
231: }
232: }
234: /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */
235: /* this segment is currently removed, so only entries in the upper triangle are obtained */
236: #ifdef column_search
237: v_i = *v + M*bs;
238: cols_i = *cols + M*bs;
239: for (i=0; i<bn; i++){ /* for each block row */
240: M = ai[i+1] - ai[i];
241: for (j=0; j<M; j++){
242: itmp = aj[ai[i] + j]; /* block column value */
243: if (itmp == bn){
244: aa_i = aa + bs2*(ai[i] + j) + bs*bp;
245: for (k=0; k<bs; k++) {
246: *cols_i++ = i*bs+k;
247: *v_i++ = aa_i[k];
248: }
249: *ncols += bs;
250: break;
251: }
252: }
253: }
254: #endif
256: return(0);
257: }
259: #undef __FUNCT__
261: int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v)
262: {
266: if (idx) {if (*idx) {PetscFree(*idx);}}
267: if (v) {if (*v) {PetscFree(*v);}}
268: return(0);
269: }
271: #undef __FUNCT__
273: int MatTranspose_SeqSBAIJ(Mat A,Mat *B)
274: {
277: MatDuplicate(A,MAT_COPY_VALUES,B);
278: return(0);
279: }
281: #undef __FUNCT__
283: static int MatView_SeqSBAIJ_Binary(Mat A,PetscViewer viewer)
284: {
285: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
286: int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2;
287: PetscScalar *aa;
288: FILE *file;
291: ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);
292: ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);
293: col_lens[0] = MAT_FILE_COOKIE;
295: col_lens[1] = A->m;
296: col_lens[2] = A->m;
297: col_lens[3] = a->s_nz*bs2;
299: /* store lengths of each row and write (including header) to file */
300: count = 0;
301: for (i=0; i<a->mbs; i++) {
302: for (j=0; j<bs; j++) {
303: col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]);
304: }
305: }
306: PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);
307: PetscFree(col_lens);
309: /* store column indices (zero start index) */
310: ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(int),&jj);
311: count = 0;
312: for (i=0; i<a->mbs; i++) {
313: for (j=0; j<bs; j++) {
314: for (k=a->i[i]; k<a->i[i+1]; k++) {
315: for (l=0; l<bs; l++) {
316: jj[count++] = bs*a->j[k] + l;
317: }
318: }
319: }
320: }
321: PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);
322: PetscFree(jj);
324: /* store nonzero values */
325: ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(PetscScalar),&aa);
326: count = 0;
327: for (i=0; i<a->mbs; i++) {
328: for (j=0; j<bs; j++) {
329: for (k=a->i[i]; k<a->i[i+1]; k++) {
330: for (l=0; l<bs; l++) {
331: aa[count++] = a->a[bs2*k + l*bs + j];
332: }
333: }
334: }
335: }
336: PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);
337: PetscFree(aa);
339: PetscViewerBinaryGetInfoPointer(viewer,&file);
340: if (file) {
341: fprintf(file,"-matload_block_size %dn",a->bs);
342: }
343: return(0);
344: }
346: #undef __FUNCT__
348: static int MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
349: {
350: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
351: int ierr,i,j,bs = a->bs,k,l,bs2=a->bs2;
352: char *name;
353: PetscViewerFormat format;
356: PetscObjectGetName((PetscObject)A,&name);
357: PetscViewerGetFormat(viewer,&format);
358: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
359: PetscViewerASCIIPrintf(viewer," block size is %dn",bs);
360: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
361: SETERRQ(PETSC_ERR_SUP,"Matlab format not supported");
362: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
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 && PetscRealPart(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 && PetscRealPart(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 if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
377: PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));
378: }
379: #else
380: if (a->a[bs2*k + l*bs + j] != 0.0) {
381: PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);
382: }
383: #endif
384: }
385: }
386: PetscViewerASCIIPrintf(viewer,"n");
387: }
388: }
389: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
390: } else {
391: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
392: for (i=0; i<a->mbs; i++) {
393: for (j=0; j<bs; j++) {
394: PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);
395: for (k=a->i[i]; k<a->i[i+1]; k++) {
396: for (l=0; l<bs; l++) {
397: #if defined(PETSC_USE_COMPLEX)
398: if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
399: PetscViewerASCIIPrintf(viewer," (%d, %g + %g i) ",bs*a->j[k]+l,
400: PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
401: } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
402: PetscViewerASCIIPrintf(viewer," (%d, %g - %g i) ",bs*a->j[k]+l,
403: PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
404: } else {
405: PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));
406: }
407: #else
408: PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);
409: #endif
410: }
411: }
412: PetscViewerASCIIPrintf(viewer,"n");
413: }
414: }
415: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
416: }
417: PetscViewerFlush(viewer);
418: return(0);
419: }
421: #undef __FUNCT__
423: static int MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
424: {
425: Mat A = (Mat) Aa;
426: Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
427: int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank;
428: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
429: MatScalar *aa;
430: MPI_Comm comm;
431: PetscViewer viewer;
434: /*
435: This is nasty. If this is called from an originally parallel matrix
436: then all processes call this,but only the first has the matrix so the
437: rest should return immediately.
438: */
439: PetscObjectGetComm((PetscObject)draw,&comm);
440: MPI_Comm_rank(comm,&rank);
441: if (rank) return(0);
443: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
445: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
446: PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");
448: /* loop over matrix elements drawing boxes */
449: color = PETSC_DRAW_BLUE;
450: for (i=0,row=0; i<mbs; i++,row+=bs) {
451: for (j=a->i[i]; j<a->i[i+1]; j++) {
452: y_l = A->m - row - 1.0; y_r = y_l + 1.0;
453: x_l = a->j[j]*bs; x_r = x_l + 1.0;
454: aa = a->a + j*bs2;
455: for (k=0; k<bs; k++) {
456: for (l=0; l<bs; l++) {
457: if (PetscRealPart(*aa++) >= 0.) continue;
458: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
459: }
460: }
461: }
462: }
463: color = PETSC_DRAW_CYAN;
464: for (i=0,row=0; i<mbs; i++,row+=bs) {
465: for (j=a->i[i]; j<a->i[i+1]; j++) {
466: y_l = A->m - row - 1.0; y_r = y_l + 1.0;
467: x_l = a->j[j]*bs; x_r = x_l + 1.0;
468: aa = a->a + j*bs2;
469: for (k=0; k<bs; k++) {
470: for (l=0; l<bs; l++) {
471: if (PetscRealPart(*aa++) != 0.) continue;
472: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
473: }
474: }
475: }
476: }
478: color = PETSC_DRAW_RED;
479: for (i=0,row=0; i<mbs; i++,row+=bs) {
480: for (j=a->i[i]; j<a->i[i+1]; j++) {
481: y_l = A->m - row - 1.0; y_r = y_l + 1.0;
482: x_l = a->j[j]*bs; x_r = x_l + 1.0;
483: aa = a->a + j*bs2;
484: for (k=0; k<bs; k++) {
485: for (l=0; l<bs; l++) {
486: if (PetscRealPart(*aa++) <= 0.) continue;
487: PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
488: }
489: }
490: }
491: }
492: return(0);
493: }
495: #undef __FUNCT__
497: static int MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
498: {
499: int ierr;
500: PetscReal xl,yl,xr,yr,w,h;
501: PetscDraw draw;
502: PetscTruth isnull;
506: PetscViewerDrawGetDraw(viewer,0,&draw);
507: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
509: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
510: xr = A->m; yr = A->m; h = yr/10.0; w = xr/10.0;
511: xr += w; yr += h; xl = -w; yl = -h;
512: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
513: PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);
514: PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
515: return(0);
516: }
518: #undef __FUNCT__
520: int MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
521: {
522: int ierr;
523: PetscTruth issocket,isascii,isbinary,isdraw;
526: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);
527: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
528: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
529: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
530: if (issocket) {
531: SETERRQ(PETSC_ERR_SUP,"Socket viewer not supported");
532: } else if (isascii){
533: MatView_SeqSBAIJ_ASCII(A,viewer);
534: } else if (isbinary) {
535: MatView_SeqSBAIJ_Binary(A,viewer);
536: } else if (isdraw) {
537: MatView_SeqSBAIJ_Draw(A,viewer);
538: } else {
539: SETERRQ1(1,"Viewer type %s not supported by SeqSBAIJ matrices",((PetscObject)viewer)->type_name);
540: }
541: return(0);
542: }
545: #undef __FUNCT__
547: int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v)
548: {
549: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
550: int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
551: int *ai = a->i,*ailen = a->ilen;
552: int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2;
553: MatScalar *ap,*aa = a->a,zero = 0.0;
556: for (k=0; k<m; k++) { /* loop over rows */
557: row = im[k]; brow = row/bs;
558: if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row");
559: if (row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
560: rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ;
561: nrow = ailen[brow];
562: for (l=0; l<n; l++) { /* loop over columns */
563: if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative column");
564: if (in[l] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
565: col = in[l] ;
566: bcol = col/bs;
567: cidx = col%bs;
568: ridx = row%bs;
569: high = nrow;
570: low = 0; /* assume unsorted */
571: while (high-low > 5) {
572: t = (low+high)/2;
573: if (rp[t] > bcol) high = t;
574: else low = t;
575: }
576: for (i=low; i<high; i++) {
577: if (rp[i] > bcol) break;
578: if (rp[i] == bcol) {
579: *v++ = ap[bs2*i+bs*cidx+ridx];
580: goto finished;
581: }
582: }
583: *v++ = zero;
584: finished:;
585: }
586: }
587: return(0);
588: }
591: #undef __FUNCT__
593: int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is)
594: {
595: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
596: int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted;
597: int *imax=a->imax,*ai=a->i,*ailen=a->ilen;
598: int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval,ierr;
599: PetscTruth roworiented=a->roworiented;
600: MatScalar *value = v,*ap,*aa = a->a,*bap;
601:
603: if (roworiented) {
604: stepval = (n-1)*bs;
605: } else {
606: stepval = (m-1)*bs;
607: }
608: for (k=0; k<m; k++) { /* loop over added rows */
609: row = im[k];
610: if (row < 0) continue;
611: #if defined(PETSC_USE_BOPT_g)
612: if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
613: #endif
614: rp = aj + ai[row];
615: ap = aa + bs2*ai[row];
616: rmax = imax[row];
617: nrow = ailen[row];
618: low = 0;
619: for (l=0; l<n; l++) { /* loop over added columns */
620: if (in[l] < 0) continue;
621: #if defined(PETSC_USE_BOPT_g)
622: if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
623: #endif
624: col = in[l];
625: if (col < row) continue; /* ignore lower triangular block */
626: if (roworiented) {
627: value = v + k*(stepval+bs)*bs + l*bs;
628: } else {
629: value = v + l*(stepval+bs)*bs + k*bs;
630: }
631: if (!sorted) low = 0; high = nrow;
632: while (high-low > 7) {
633: t = (low+high)/2;
634: if (rp[t] > col) high = t;
635: else low = t;
636: }
637: for (i=low; i<high; i++) {
638: if (rp[i] > col) break;
639: if (rp[i] == col) {
640: bap = ap + bs2*i;
641: if (roworiented) {
642: if (is == ADD_VALUES) {
643: for (ii=0; ii<bs; ii++,value+=stepval) {
644: for (jj=ii; jj<bs2; jj+=bs) {
645: bap[jj] += *value++;
646: }
647: }
648: } else {
649: for (ii=0; ii<bs; ii++,value+=stepval) {
650: for (jj=ii; jj<bs2; jj+=bs) {
651: bap[jj] = *value++;
652: }
653: }
654: }
655: } else {
656: if (is == ADD_VALUES) {
657: for (ii=0; ii<bs; ii++,value+=stepval) {
658: for (jj=0; jj<bs; jj++) {
659: *bap++ += *value++;
660: }
661: }
662: } else {
663: for (ii=0; ii<bs; ii++,value+=stepval) {
664: for (jj=0; jj<bs; jj++) {
665: *bap++ = *value++;
666: }
667: }
668: }
669: }
670: goto noinsert2;
671: }
672: }
673: if (nonew == 1) goto noinsert2;
674: else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
675: if (nrow >= rmax) {
676: /* there is no extra room in row, therefore enlarge */
677: int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j;
678: MatScalar *new_a;
680: if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
682: /* malloc new storage space */
683: len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int);
684: ierr = PetscMalloc(len,&new_a);
685: new_j = (int*)(new_a + bs2*new_nz);
686: new_i = new_j + new_nz;
688: /* copy over old data into new slots */
689: for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];}
690: for (ii=row+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;}
691: PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));
692: len = (new_nz - CHUNKSIZE - ai[row] - nrow);
693: PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));
694: PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar));
695: PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));
696: PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar));
697: /* free up old matrix storage */
698: PetscFree(a->a);
699: if (!a->singlemalloc) {
700: PetscFree(a->i);
701: PetscFree(a->j);
702: }
703: aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j;
704: a->singlemalloc = PETSC_TRUE;
706: rp = aj + ai[row]; ap = aa + bs2*ai[row];
707: rmax = imax[row] = imax[row] + CHUNKSIZE;
708: PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar)));
709: a->s_maxnz += bs2*CHUNKSIZE;
710: a->reallocs++;
711: a->s_nz++;
712: }
713: N = nrow++ - 1;
714: /* shift up all the later entries in this row */
715: for (ii=N; ii>=i; ii--) {
716: rp[ii+1] = rp[ii];
717: PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));
718: }
719: if (N >= i) {
720: PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));
721: }
722: rp[i] = col;
723: bap = ap + bs2*i;
724: if (roworiented) {
725: for (ii=0; ii<bs; ii++,value+=stepval) {
726: for (jj=ii; jj<bs2; jj+=bs) {
727: bap[jj] = *value++;
728: }
729: }
730: } else {
731: for (ii=0; ii<bs; ii++,value+=stepval) {
732: for (jj=0; jj<bs; jj++) {
733: *bap++ = *value++;
734: }
735: }
736: }
737: noinsert2:;
738: low = i;
739: }
740: ailen[row] = nrow;
741: }
742: return(0);
743: }
745: #undef __FUNCT__
747: int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
748: {
749: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
750: int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
751: int m = A->m,*ip,N,*ailen = a->ilen;
752: int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr;
753: MatScalar *aa = a->a,*ap;
754: #if defined(PETSC_HAVE_SPOOLES)
755: PetscTruth flag;
756: #endif
759: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
761: if (m) rmax = ailen[0];
762: for (i=1; i<mbs; i++) {
763: /* move each row back by the amount of empty slots (fshift) before it*/
764: fshift += imax[i-1] - ailen[i-1];
765: rmax = PetscMax(rmax,ailen[i]);
766: if (fshift) {
767: ip = aj + ai[i]; ap = aa + bs2*ai[i];
768: N = ailen[i];
769: for (j=0; j<N; j++) {
770: ip[j-fshift] = ip[j];
771: PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));
772: }
773: }
774: ai[i] = ai[i-1] + ailen[i-1];
775: }
776: if (mbs) {
777: fshift += imax[mbs-1] - ailen[mbs-1];
778: ai[mbs] = ai[mbs-1] + ailen[mbs-1];
779: }
780: /* reset ilen and imax for each row */
781: for (i=0; i<mbs; i++) {
782: ailen[i] = imax[i] = ai[i+1] - ai[i];
783: }
784: a->s_nz = ai[mbs];
786: /* diagonals may have moved, reset it */
787: if (a->diag) {
788: PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(int));
789: }
790: PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d usedn",
791: m,A->m,a->bs,fshift*bs2,a->s_nz*bs2);
792: PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %dn",
793: a->reallocs);
794: PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %dn",rmax);
795: a->reallocs = 0;
796: A->info.nz_unneeded = (PetscReal)fshift*bs2;
797:
798: #if defined(PETSC_HAVE_SPOOLES)
799: PetscOptionsHasName(A->prefix,"-mat_sbaij_spooles",&flag);
800: if (flag) { MatUseSpooles_SeqSBAIJ(A); }
801: #endif
803: return(0);
804: }
806: /*
807: This function returns an array of flags which indicate the locations of contiguous
808: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
809: then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
810: Assume: sizes should be long enough to hold all the values.
811: */
812: #undef __FUNCT__
814: int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max)
815: {
816: int i,j,k,row;
817: PetscTruth flg;
820: for (i=0,j=0; i<n; j++) {
821: row = idx[i];
822: if (row%bs!=0) { /* Not the begining of a block */
823: sizes[j] = 1;
824: i++;
825: } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
826: sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
827: i++;
828: } else { /* Begining of the block, so check if the complete block exists */
829: flg = PETSC_TRUE;
830: for (k=1; k<bs; k++) {
831: if (row+k != idx[i+k]) { /* break in the block */
832: flg = PETSC_FALSE;
833: break;
834: }
835: }
836: if (flg == PETSC_TRUE) { /* No break in the bs */
837: sizes[j] = bs;
838: i+= bs;
839: } else {
840: sizes[j] = 1;
841: i++;
842: }
843: }
844: }
845: *bs_max = j;
846: return(0);
847: }
848:
849: #undef __FUNCT__
851: int MatZeroRows_SeqSBAIJ(Mat A,IS is,PetscScalar *diag)
852: {
854: SETERRQ(PETSC_ERR_SUP,"No support for this function yet");
855: }
857: /* Only add/insert a(i,j) with i<=j (blocks).
858: Any a(i,j) with i>j input by user is ingored.
859: */
861: #undef __FUNCT__
863: int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is)
864: {
865: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
866: int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted;
867: int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
868: int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol;
869: int ridx,cidx,bs2=a->bs2,ierr;
870: MatScalar *ap,value,*aa=a->a,*bap;
874: for (k=0; k<m; k++) { /* loop over added rows */
875: row = im[k]; /* row number */
876: brow = row/bs; /* block row number */
877: if (row < 0) continue;
878: #if defined(PETSC_USE_BOPT_g)
879: if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m);
880: #endif
881: rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
882: ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
883: rmax = imax[brow]; /* maximum space allocated for this row */
884: nrow = ailen[brow]; /* actual length of this row */
885: low = 0;
887: for (l=0; l<n; l++) { /* loop over added columns */
888: if (in[l] < 0) continue;
889: #if defined(PETSC_USE_BOPT_g)
890: if (in[l] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->m);
891: #endif
892: col = in[l];
893: bcol = col/bs; /* block col number */
895: if (brow <= bcol){
896: ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
897: if ((brow==bcol && ridx<=cidx) || (brow<bcol)){
898: /* element value a(k,l) */
899: if (roworiented) {
900: value = v[l + k*n];
901: } else {
902: value = v[k + l*m];
903: }
905: /* move pointer bap to a(k,l) quickly and add/insert value */
906: if (!sorted) low = 0; high = nrow;
907: while (high-low > 7) {
908: t = (low+high)/2;
909: if (rp[t] > bcol) high = t;
910: else low = t;
911: }
912: for (i=low; i<high; i++) {
913: /* printf("The loop of i=low.., rp[%d]=%dn",i,rp[i]); */
914: if (rp[i] > bcol) break;
915: if (rp[i] == bcol) {
916: bap = ap + bs2*i + bs*cidx + ridx;
917: if (is == ADD_VALUES) *bap += value;
918: else *bap = value;
919: /* for diag block, add/insert its symmetric element a(cidx,ridx) */
920: if (brow == bcol && ridx < cidx){
921: bap = ap + bs2*i + bs*ridx + cidx;
922: if (is == ADD_VALUES) *bap += value;
923: else *bap = value;
924: }
925: goto noinsert1;
926: }
927: }
928:
929: if (nonew == 1) goto noinsert1;
930: else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
931: if (nrow >= rmax) {
932: /* there is no extra room in row, therefore enlarge */
933: int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j;
934: MatScalar *new_a;
936: if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
938: /* Malloc new storage space */
939: len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int);
940: ierr = PetscMalloc(len,&new_a);
941: new_j = (int*)(new_a + bs2*new_nz);
942: new_i = new_j + new_nz;
944: /* copy over old data into new slots */
945: for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];}
946: for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;}
947: PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));
948: len = (new_nz - CHUNKSIZE - ai[brow] - nrow);
949: PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));
950: PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));
951: PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));
952: PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));
953: /* free up old matrix storage */
954: PetscFree(a->a);
955: if (!a->singlemalloc) {
956: PetscFree(a->i);
957: PetscFree(a->j);
958: }
959: aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j;
960: a->singlemalloc = PETSC_TRUE;
962: rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
963: rmax = imax[brow] = imax[brow] + CHUNKSIZE;
964: PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar)));
965: a->s_maxnz += bs2*CHUNKSIZE;
966: a->reallocs++;
967: a->s_nz++;
968: }
970: N = nrow++ - 1;
971: /* shift up all the later entries in this row */
972: for (ii=N; ii>=i; ii--) {
973: rp[ii+1] = rp[ii];
974: ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));
975: }
976: if (N>=i) {
977: PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));
978: }
979: rp[i] = bcol;
980: ap[bs2*i + bs*cidx + ridx] = value;
981: noinsert1:;
982: low = i;
983: /* } */
984: }
985: } /* end of if .. if.. */
986: } /* end of loop over added columns */
987: ailen[brow] = nrow;
988: } /* end of loop over added rows */
990: return(0);
991: }
993: extern int MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*);
994: extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal);
995: extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int);
996: extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*);
997: extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**);
998: extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec);
999: extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec);
1000: extern int MatScale_SeqSBAIJ(PetscScalar*,Mat);
1001: extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *);
1002: extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*);
1003: extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec);
1004: extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec);
1005: extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *);
1006: extern int MatZeroEntries_SeqSBAIJ(Mat);
1007: extern int MatGetRowMax_SeqSBAIJ(Mat,Vec);
1009: extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec);
1010: extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec);
1011: extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec);
1012: extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec);
1013: extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec);
1014: extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec);
1015: extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec);
1016: extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec);
1017: extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec);
1018: extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec);
1019: extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec);
1020: extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec);
1021: extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec);
1022: extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec);
1023: extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec);
1025: extern int MatSolves_SeqSBAIJ_1(Mat,Vecs,Vecs);
1027: extern int MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec);
1028: extern int MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec);
1029: extern int MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec);
1030: extern int MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec);
1031: extern int MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec);
1032: extern int MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec);
1033: extern int MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec);
1034: extern int MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec);
1036: extern int MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*);
1037: extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*);
1038: extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*);
1039: extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*);
1040: extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*);
1041: extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*);
1042: extern int MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,Mat*);
1043: extern int MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,Mat*);
1044: extern int MatGetInertia_SeqSBAIJ(Mat,int*,int*,int*);
1046: extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec);
1047: extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec);
1048: extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec);
1049: extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec);
1050: extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec);
1051: extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec);
1052: extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec);
1053: extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec);
1055: extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec);
1056: extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec);
1057: extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec);
1058: extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec);
1059: extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec);
1060: extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec);
1061: extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec);
1062: extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec);
1064: #ifdef HAVE_MatICCFactor
1065: /* modefied from MatILUFactor_SeqSBAIJ, needs further work! */
1066: #undef __FUNCT__
1068: int MatICCFactor_SeqSBAIJ(Mat inA,IS row,PetscReal fill,int level)
1069: {
1070: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1071: Mat outA;
1072: int ierr;
1073: PetscTruth row_identity,col_identity;
1076: /*
1077: if (level != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU");
1078: ISIdentity(row,&row_identity);
1079: ISIdentity(col,&col_identity);
1080: if (!row_identity || !col_identity) {
1081: SETERRQ(1,"Row and column permutations must be identity for in-place ICC");
1082: }
1083: */
1085: outA = inA;
1086: inA->factor = FACTOR_CHOLESKY;
1088: if (!a->diag) {
1089: MatMarkDiagonal_SeqSBAIJ(inA);
1090: }
1091: /*
1092: Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver
1093: for ILU(0) factorization with natural ordering
1094: */
1095: switch (a->bs) {
1096: case 1:
1097: inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1098: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering;
1099: inA->ops->solves = MatSolves_SeqSBAIJ_1;
1100: PetscLoginfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering solvetrans BS=1n");
1101: case 2:
1102: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1103: inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering;
1104: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering;
1105: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=2n");
1106: break;
1107: case 3:
1108: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1109: inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering;
1110: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering;
1111: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=3n");
1112: break;
1113: case 4:
1114: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1115: inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering;
1116: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering;
1117: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=4n");
1118: break;
1119: case 5:
1120: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1121: inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering;
1122: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering;
1123: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=5n");
1124: break;
1125: case 6:
1126: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1127: inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering;
1128: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering;
1129: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=6n");
1130: break;
1131: case 7:
1132: inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1133: inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering;
1134: inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering;
1135: PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=7n");
1136: break;
1137: default:
1138: a->row = row;
1139: a->icol = col;
1140: ierr = PetscObjectReference((PetscObject)row);
1141: ierr = PetscObjectReference((PetscObject)col);
1143: /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */
1144: ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));
1145: PetscLogObjectParent(inA,a->icol);
1147: if (!a->solve_work) {
1148: PetscMalloc((A->m+a->bs)*sizeof(PetscScalar),&a->solve_work);
1149: PetscLogObjectMemory(inA,(A->m+a->bs)*sizeof(PetscScalar));
1150: }
1151: }
1153: MatCholeskyFactorNumeric(inA,&outA);
1155: return(0);
1156: }
1157: #endif
1159: #undef __FUNCT__
1161: int MatPrintHelp_SeqSBAIJ(Mat A)
1162: {
1163: static PetscTruth called = PETSC_FALSE;
1164: MPI_Comm comm = A->comm;
1165: int ierr;
1168: if (called) {return(0);} else called = PETSC_TRUE;
1169: (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):n");
1170: (*PetscHelpPrintf)(comm," -mat_block_size <block_size>n");
1171: return(0);
1172: }
1174: EXTERN_C_BEGIN
1175: #undef __FUNCT__
1177: int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices)
1178: {
1179: Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data;
1180: int i,nz,n;
1183: nz = baij->s_maxnz;
1184: n = mat->n;
1185: for (i=0; i<nz; i++) {
1186: baij->j[i] = indices[i];
1187: }
1188: baij->s_nz = nz;
1189: for (i=0; i<n; i++) {
1190: baij->ilen[i] = baij->imax[i];
1191: }
1193: return(0);
1194: }
1195: EXTERN_C_END
1197: #undef __FUNCT__
1199: /*@
1200: MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1201: in the matrix.
1203: Input Parameters:
1204: + mat - the SeqSBAIJ matrix
1205: - indices - the column indices
1207: Level: advanced
1209: Notes:
1210: This can be called if you have precomputed the nonzero structure of the
1211: matrix and want to provide it to the matrix object to improve the performance
1212: of the MatSetValues() operation.
1214: You MUST have set the correct numbers of nonzeros per row in the call to
1215: MatCreateSeqSBAIJ().
1217: MUST be called before any calls to MatSetValues()
1219: .seealso: MatCreateSeqSBAIJ
1220: @*/
1221: int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices)
1222: {
1223: int ierr,(*f)(Mat,int *);
1227: PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);
1228: if (f) {
1229: (*f)(mat,indices);
1230: } else {
1231: SETERRQ(1,"Wrong type of matrix to set column indices");
1232: }
1233: return(0);
1234: }
1236: #undef __FUNCT__
1238: int MatSetUpPreallocation_SeqSBAIJ(Mat A)
1239: {
1240: int ierr;
1243: MatSeqSBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);
1244: return(0);
1245: }
1247: #undef __FUNCT__
1249: int MatGetArray_SeqSBAIJ(Mat A,PetscScalar **array)
1250: {
1251: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1253: *array = a->a;
1254: return(0);
1255: }
1257: #undef __FUNCT__
1259: int MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar **array)
1260: {
1262: return(0);
1263: }
1265: #include petscblaslapack.h
1266: #undef __FUNCT__
1268: int MatAXPY_SeqSBAIJ(PetscScalar *a,Mat X,Mat Y,MatStructure str)
1269: {
1270: int ierr,one=1;
1271: Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data,*y = (Mat_SeqSBAIJ *)Y->data;
1274: if (str == SAME_NONZERO_PATTERN) {
1275: BLaxpy_(&x->s_nz,a,x->a,&one,y->a,&one);
1276: } else {
1277: MatAXPY_Basic(a,X,Y,str);
1278: }
1279: return(0);
1280: }
1282: /* -------------------------------------------------------------------*/
1283: static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1284: MatGetRow_SeqSBAIJ,
1285: MatRestoreRow_SeqSBAIJ,
1286: MatMult_SeqSBAIJ_N,
1287: MatMultAdd_SeqSBAIJ_N,
1288: MatMultTranspose_SeqSBAIJ,
1289: MatMultTransposeAdd_SeqSBAIJ,
1290: MatSolve_SeqSBAIJ_N,
1291: 0,
1292: 0,
1293: 0,
1294: 0,
1295: MatCholeskyFactor_SeqSBAIJ,
1296: MatRelax_SeqSBAIJ,
1297: MatTranspose_SeqSBAIJ,
1298: MatGetInfo_SeqSBAIJ,
1299: MatEqual_SeqSBAIJ,
1300: MatGetDiagonal_SeqSBAIJ,
1301: MatDiagonalScale_SeqSBAIJ,
1302: MatNorm_SeqSBAIJ,
1303: 0,
1304: MatAssemblyEnd_SeqSBAIJ,
1305: 0,
1306: MatSetOption_SeqSBAIJ,
1307: MatZeroEntries_SeqSBAIJ,
1308: MatZeroRows_SeqSBAIJ,
1309: 0,
1310: 0,
1311: MatCholeskyFactorSymbolic_SeqSBAIJ,
1312: MatCholeskyFactorNumeric_SeqSBAIJ_N,
1313: MatSetUpPreallocation_SeqSBAIJ,
1314: 0,
1315: MatICCFactorSymbolic_SeqSBAIJ,
1316: MatGetArray_SeqSBAIJ,
1317: MatRestoreArray_SeqSBAIJ,
1318: MatDuplicate_SeqSBAIJ,
1319: 0,
1320: 0,
1321: 0,
1322: 0,
1323: MatAXPY_SeqSBAIJ,
1324: MatGetSubMatrices_SeqSBAIJ,
1325: MatIncreaseOverlap_SeqSBAIJ,
1326: MatGetValues_SeqSBAIJ,
1327: 0,
1328: MatPrintHelp_SeqSBAIJ,
1329: MatScale_SeqSBAIJ,
1330: 0,
1331: 0,
1332: 0,
1333: MatGetBlockSize_SeqSBAIJ,
1334: MatGetRowIJ_SeqSBAIJ,
1335: MatRestoreRowIJ_SeqSBAIJ,
1336: 0,
1337: 0,
1338: 0,
1339: 0,
1340: 0,
1341: 0,
1342: MatSetValuesBlocked_SeqSBAIJ,
1343: MatGetSubMatrix_SeqSBAIJ,
1344: 0,
1345: 0,
1346: MatGetPetscMaps_Petsc,
1347: 0,
1348: 0,
1349: 0,
1350: 0,
1351: 0,
1352: 0,
1353: MatGetRowMax_SeqSBAIJ,
1354: 0,
1355: 0,
1356: 0,
1357: 0,
1358: 0,
1359: 0,
1360: 0,
1361: 0,
1362: 0,
1363: 0,
1364: 0,
1365: 0,
1366: 0,
1367: #if !defined(PETSC_USE_COMPLEX)
1368: MatGetInertia_SeqSBAIJ
1369: #else
1370: 0
1371: #endif
1372: };
1374: EXTERN_C_BEGIN
1375: #undef __FUNCT__
1377: int MatStoreValues_SeqSBAIJ(Mat mat)
1378: {
1379: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1380: int nz = aij->i[mat->m]*aij->bs*aij->bs2;
1381: int ierr;
1384: if (aij->nonew != 1) {
1385: SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first");
1386: }
1388: /* allocate space for values if not already there */
1389: if (!aij->saved_values) {
1390: PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);
1391: }
1393: /* copy values over */
1394: PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));
1395: return(0);
1396: }
1397: EXTERN_C_END
1399: EXTERN_C_BEGIN
1400: #undef __FUNCT__
1402: int MatRetrieveValues_SeqSBAIJ(Mat mat)
1403: {
1404: Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data;
1405: int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr;
1408: if (aij->nonew != 1) {
1409: SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first");
1410: }
1411: if (!aij->saved_values) {
1412: SETERRQ(1,"Must call MatStoreValues(A);first");
1413: }
1415: /* copy values over */
1416: PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));
1417: return(0);
1418: }
1419: EXTERN_C_END
1421: EXTERN_C_BEGIN
1422: #undef __FUNCT__
1424: int MatCreate_SeqSBAIJ(Mat B)
1425: {
1426: Mat_SeqSBAIJ *b;
1427: int ierr,size;
1430: MPI_Comm_size(B->comm,&size);
1431: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1432: B->m = B->M = PetscMax(B->m,B->M);
1433: B->n = B->N = PetscMax(B->n,B->N);
1435: ierr = PetscNew(Mat_SeqSBAIJ,&b);
1436: B->data = (void*)b;
1437: ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));
1438: ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1439: B->ops->destroy = MatDestroy_SeqSBAIJ;
1440: B->ops->view = MatView_SeqSBAIJ;
1441: B->factor = 0;
1442: B->lupivotthreshold = 1.0;
1443: B->mapping = 0;
1444: b->row = 0;
1445: b->icol = 0;
1446: b->reallocs = 0;
1447: b->saved_values = 0;
1448:
1449: PetscMapCreateMPI(B->comm,B->m,B->M,&B->rmap);
1450: PetscMapCreateMPI(B->comm,B->n,B->N,&B->cmap);
1452: b->sorted = PETSC_FALSE;
1453: b->roworiented = PETSC_TRUE;
1454: b->nonew = 0;
1455: b->diag = 0;
1456: b->solve_work = 0;
1457: b->mult_work = 0;
1458: B->spptr = 0;
1459: b->keepzeroedrows = PETSC_FALSE;
1460:
1461: b->inew = 0;
1462: b->jnew = 0;
1463: b->anew = 0;
1464: b->a2anew = 0;
1465: b->permute = PETSC_FALSE;
1467: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C",
1468: "MatStoreValues_SeqSBAIJ",
1469: (void*)MatStoreValues_SeqSBAIJ);
1470: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C",
1471: "MatRetrieveValues_SeqSBAIJ",
1472: (void*)MatRetrieveValues_SeqSBAIJ);
1473: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",
1474: "MatSeqSBAIJSetColumnIndices_SeqSBAIJ",
1475: (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);
1476: return(0);
1477: }
1478: EXTERN_C_END
1480: #undef __FUNCT__
1482: /*@C
1483: MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
1484: compressed row) format. For good matrix assembly performance the
1485: user should preallocate the matrix storage by setting the parameter nz
1486: (or the array nnz). By setting these parameters accurately, performance
1487: during matrix assembly can be increased by more than a factor of 50.
1489: Collective on Mat
1491: Input Parameters:
1492: + A - the symmetric matrix
1493: . bs - size of block
1494: . nz - number of block nonzeros per block row (same for all rows)
1495: - nnz - array containing the number of block nonzeros in the upper triangular plus
1496: diagonal portion of each block (possibly different for each block row) or PETSC_NULL
1498: Options Database Keys:
1499: . -mat_no_unroll - uses code that does not unroll the loops in the
1500: block calculations (much slower)
1501: . -mat_block_size - size of the blocks to use
1503: Level: intermediate
1505: Notes:
1506: Specify the preallocated storage with either nz or nnz (not both).
1507: Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory
1508: allocation. For additional details, see the users manual chapter on
1509: matrices.
1511: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ()
1512: @*/
1513: int MatSeqSBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz)
1514: {
1515: Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1516: int i,len,ierr,mbs,bs2;
1517: PetscTruth flg;
1518: int s_nz;
1521: B->preallocated = PETSC_TRUE;
1522: PetscOptionsGetInt(B->prefix,"-mat_block_size",&bs,PETSC_NULL);
1523: mbs = B->m/bs;
1524: bs2 = bs*bs;
1526: if (mbs*bs != B->m) {
1527: SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize");
1528: }
1530: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1531: if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz);
1532: if (nnz) {
1533: for (i=0; i<mbs; i++) {
1534: if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]);
1535: 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);
1536: }
1537: }
1539: ierr = PetscOptionsHasName(B->prefix,"-mat_no_unroll",&flg);
1540: if (!flg) {
1541: switch (bs) {
1542: case 1:
1543: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1;
1544: B->ops->solve = MatSolve_SeqSBAIJ_1;
1545: B->ops->solves = MatSolves_SeqSBAIJ_1;
1546: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1;
1547: B->ops->mult = MatMult_SeqSBAIJ_1;
1548: B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1549: break;
1550: case 2:
1551: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1552: B->ops->solve = MatSolve_SeqSBAIJ_2;
1553: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2;
1554: B->ops->mult = MatMult_SeqSBAIJ_2;
1555: B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1556: break;
1557: case 3:
1558: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1559: B->ops->solve = MatSolve_SeqSBAIJ_3;
1560: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3;
1561: B->ops->mult = MatMult_SeqSBAIJ_3;
1562: B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1563: break;
1564: case 4:
1565: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1566: B->ops->solve = MatSolve_SeqSBAIJ_4;
1567: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4;
1568: B->ops->mult = MatMult_SeqSBAIJ_4;
1569: B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1570: break;
1571: case 5:
1572: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1573: B->ops->solve = MatSolve_SeqSBAIJ_5;
1574: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5;
1575: B->ops->mult = MatMult_SeqSBAIJ_5;
1576: B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1577: break;
1578: case 6:
1579: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1580: B->ops->solve = MatSolve_SeqSBAIJ_6;
1581: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6;
1582: B->ops->mult = MatMult_SeqSBAIJ_6;
1583: B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1584: break;
1585: case 7:
1586: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1587: B->ops->solve = MatSolve_SeqSBAIJ_7;
1588: B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7;
1589: B->ops->mult = MatMult_SeqSBAIJ_7;
1590: B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1591: break;
1592: }
1593: }
1595: b->mbs = mbs;
1596: b->nbs = mbs;
1597: ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);
1598: if (!nnz) {
1599: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1600: else if (nz <= 0) nz = 1;
1601: for (i=0; i<mbs; i++) {
1602: b->imax[i] = nz;
1603: }
1604: nz = nz*mbs; /* total nz */
1605: } else {
1606: nz = 0;
1607: for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];}
1608: }
1609: /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1610: s_nz = nz;
1612: /* allocate the matrix space */
1613: len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int);
1614: PetscMalloc(len,&b->a);
1615: PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));
1616: b->j = (int*)(b->a + s_nz*bs2);
1617: PetscMemzero(b->j,s_nz*sizeof(int));
1618: b->i = b->j + s_nz;
1619: b->singlemalloc = PETSC_TRUE;
1621: /* pointer to beginning of each row */
1622: b->i[0] = 0;
1623: for (i=1; i<mbs+1; i++) {
1624: b->i[i] = b->i[i-1] + b->imax[i-1];
1625: }
1627: /* b->ilen will count nonzeros in each block row so far. */
1628: PetscMalloc((mbs+1)*sizeof(int),&b->ilen);
1629: PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ));
1630: for (i=0; i<mbs; i++) { b->ilen[i] = 0;}
1631:
1632: b->bs = bs;
1633: b->bs2 = bs2;
1634: b->s_nz = 0;
1635: b->s_maxnz = s_nz*bs2;
1636:
1637: b->inew = 0;
1638: b->jnew = 0;
1639: b->anew = 0;
1640: b->a2anew = 0;
1641: b->permute = PETSC_FALSE;
1642: return(0);
1643: }
1646: #undef __FUNCT__
1648: /*@C
1649: MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
1650: compressed row) format. For good matrix assembly performance the
1651: user should preallocate the matrix storage by setting the parameter nz
1652: (or the array nnz). By setting these parameters accurately, performance
1653: during matrix assembly can be increased by more than a factor of 50.
1655: Collective on MPI_Comm
1657: Input Parameters:
1658: + comm - MPI communicator, set to PETSC_COMM_SELF
1659: . bs - size of block
1660: . m - number of rows, or number of columns
1661: . nz - number of block nonzeros per block row (same for all rows)
1662: - nnz - array containing the number of block nonzeros in the upper triangular plus
1663: diagonal portion of each block (possibly different for each block row) or PETSC_NULL
1665: Output Parameter:
1666: . A - the symmetric matrix
1668: Options Database Keys:
1669: . -mat_no_unroll - uses code that does not unroll the loops in the
1670: block calculations (much slower)
1671: . -mat_block_size - size of the blocks to use
1673: Level: intermediate
1675: Notes:
1677: Specify the preallocated storage with either nz or nnz (not both).
1678: Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory
1679: allocation. For additional details, see the users manual chapter on
1680: matrices.
1682: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ()
1683: @*/
1684: int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A)
1685: {
1687:
1689: MatCreate(comm,m,n,m,n,A);
1690: MatSetType(*A,MATSEQSBAIJ);
1691: MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);
1692: return(0);
1693: }
1695: #undef __FUNCT__
1697: int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
1698: {
1699: Mat C;
1700: Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
1701: int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr;
1704: if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix");
1706: *B = 0;
1707: MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);
1708: MatSetType(C,MATSEQSBAIJ);
1709: c = (Mat_SeqSBAIJ*)C->data;
1711: ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));
1712: C->preallocated = PETSC_TRUE;
1713: C->factor = A->factor;
1714: c->row = 0;
1715: c->icol = 0;
1716: c->saved_values = 0;
1717: c->keepzeroedrows = a->keepzeroedrows;
1718: C->assembled = PETSC_TRUE;
1720: c->bs = a->bs;
1721: c->bs2 = a->bs2;
1722: c->mbs = a->mbs;
1723: c->nbs = a->nbs;
1725: PetscMalloc((mbs+1)*sizeof(int),&c->imax);
1726: PetscMalloc((mbs+1)*sizeof(int),&c->ilen);
1727: for (i=0; i<mbs; i++) {
1728: c->imax[i] = a->imax[i];
1729: c->ilen[i] = a->ilen[i];
1730: }
1732: /* allocate the matrix space */
1733: c->singlemalloc = PETSC_TRUE;
1734: len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int));
1735: PetscMalloc(len,&c->a);
1736: c->j = (int*)(c->a + nz*bs2);
1737: c->i = c->j + nz;
1738: PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));
1739: if (mbs > 0) {
1740: PetscMemcpy(c->j,a->j,nz*sizeof(int));
1741: if (cpvalues == MAT_COPY_VALUES) {
1742: PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));
1743: } else {
1744: PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));
1745: }
1746: }
1748: PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ));
1749: c->sorted = a->sorted;
1750: c->roworiented = a->roworiented;
1751: c->nonew = a->nonew;
1753: if (a->diag) {
1754: PetscMalloc((mbs+1)*sizeof(int),&c->diag);
1755: PetscLogObjectMemory(C,(mbs+1)*sizeof(int));
1756: for (i=0; i<mbs; i++) {
1757: c->diag[i] = a->diag[i];
1758: }
1759: } else c->diag = 0;
1760: c->s_nz = a->s_nz;
1761: c->s_maxnz = a->s_maxnz;
1762: c->solve_work = 0;
1763: C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */
1764: c->mult_work = 0;
1765: *B = C;
1766: PetscFListDuplicate(A->qlist,&C->qlist);
1767: return(0);
1768: }
1770: EXTERN_C_BEGIN
1771: #undef __FUNCT__
1773: int MatLoad_SeqSBAIJ(PetscViewer viewer,MatType type,Mat *A)
1774: {
1775: Mat_SeqSBAIJ *a;
1776: Mat B;
1777: int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1;
1778: int *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount;
1779: int kmax,jcount,block,idx,point,nzcountb,extra_rows;
1780: int *masked,nmask,tmp,bs2,ishift;
1781: PetscScalar *aa;
1782: MPI_Comm comm = ((PetscObject)viewer)->comm;
1785: PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);
1786: bs2 = bs*bs;
1788: MPI_Comm_size(comm,&size);
1789: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor");
1790: PetscViewerBinaryGetDescriptor(viewer,&fd);
1791: PetscBinaryRead(fd,header,4,PETSC_INT);
1792: if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object");
1793: M = header[1]; N = header[2]; nz = header[3];
1795: if (header[3] < 0) {
1796: SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ");
1797: }
1799: if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices");
1801: /*
1802: This code adds extra rows to make sure the number of rows is
1803: divisible by the blocksize
1804: */
1805: mbs = M/bs;
1806: extra_rows = bs - M + bs*(mbs);
1807: if (extra_rows == bs) extra_rows = 0;
1808: else mbs++;
1809: if (extra_rows) {
1810: PetscLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksizen");
1811: }
1813: /* read in row lengths */
1814: PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);
1815: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1816: for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;
1818: /* read in column indices */
1819: PetscMalloc((nz+extra_rows)*sizeof(int),&jj);
1820: PetscBinaryRead(fd,jj,nz,PETSC_INT);
1821: for (i=0; i<extra_rows; i++) jj[nz+i] = M+i;
1823: /* loop over row lengths determining block row lengths */
1824: ierr = PetscMalloc(mbs*sizeof(int),&s_browlengths);
1825: ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));
1826: ierr = PetscMalloc(2*mbs*sizeof(int),&mask);
1827: ierr = PetscMemzero(mask,mbs*sizeof(int));
1828: masked = mask + mbs;
1829: rowcount = 0; nzcount = 0;
1830: for (i=0; i<mbs; i++) {
1831: nmask = 0;
1832: for (j=0; j<bs; j++) {
1833: kmax = rowlengths[rowcount];
1834: for (k=0; k<kmax; k++) {
1835: tmp = jj[nzcount++]/bs; /* block col. index */
1836: if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;}
1837: }
1838: rowcount++;
1839: }
1840: s_browlengths[i] += nmask;
1841:
1842: /* zero out the mask elements we set */
1843: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
1844: }
1846: /* create our matrix */
1847: MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,s_browlengths,A);
1848: B = *A;
1849: a = (Mat_SeqSBAIJ*)B->data;
1851: /* set matrix "i" values */
1852: a->i[0] = 0;
1853: for (i=1; i<= mbs; i++) {
1854: a->i[i] = a->i[i-1] + s_browlengths[i-1];
1855: a->ilen[i-1] = s_browlengths[i-1];
1856: }
1857: a->s_nz = a->i[mbs];
1859: /* read in nonzero values */
1860: PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);
1861: PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);
1862: for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0;
1864: /* set "a" and "j" values into matrix */
1865: nzcount = 0; jcount = 0;
1866: for (i=0; i<mbs; i++) {
1867: nzcountb = nzcount;
1868: nmask = 0;
1869: for (j=0; j<bs; j++) {
1870: kmax = rowlengths[i*bs+j];
1871: for (k=0; k<kmax; k++) {
1872: tmp = jj[nzcount++]/bs; /* block col. index */
1873: if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;}
1874: }
1875: }
1876: /* sort the masked values */
1877: PetscSortInt(nmask,masked);
1879: /* set "j" values into matrix */
1880: maskcount = 1;
1881: for (j=0; j<nmask; j++) {
1882: a->j[jcount++] = masked[j];
1883: mask[masked[j]] = maskcount++;
1884: }
1886: /* set "a" values into matrix */
1887: ishift = bs2*a->i[i];
1888: for (j=0; j<bs; j++) {
1889: kmax = rowlengths[i*bs+j];
1890: for (k=0; k<kmax; k++) {
1891: tmp = jj[nzcountb]/bs ; /* block col. index */
1892: if (tmp >= i){
1893: block = mask[tmp] - 1;
1894: point = jj[nzcountb] - bs*tmp;
1895: idx = ishift + bs2*block + j + bs*point;
1896: a->a[idx] = aa[nzcountb];
1897: }
1898: nzcountb++;
1899: }
1900: }
1901: /* zero out the mask elements we set */
1902: for (j=0; j<nmask; j++) mask[masked[j]] = 0;
1903: }
1904: if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix");
1906: PetscFree(rowlengths);
1907: PetscFree(s_browlengths);
1908: PetscFree(aa);
1909: PetscFree(jj);
1910: PetscFree(mask);
1912: B->assembled = PETSC_TRUE;
1913: MatView_Private(B);
1914: return(0);
1915: }
1916: EXTERN_C_END
1918: #undef __FUNCT__
1920: int MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx)
1921: {
1922: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1923: MatScalar *aa=a->a,*v,*v1;
1924: PetscScalar *x,*b,*t,sum,d;
1925: int m=a->mbs,bs=a->bs,*ai=a->i,*aj=a->j,ierr;
1926: int nz,nz1,*vj,*vj1,i;
1929: its = its*lits;
1930: if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits);
1932: if (bs > 1)
1933: SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented");
1935: VecGetArray(xx,&x);
1936: if (xx != bb) {
1937: VecGetArray(bb,&b);
1938: } else {
1939: b = x;
1940: }
1942: PetscMalloc(m*sizeof(PetscScalar),&t);
1943:
1944: if (flag & SOR_ZERO_INITIAL_GUESS) {
1945: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
1946: for (i=0; i<m; i++)
1947: t[i] = b[i];
1949: for (i=0; i<m; i++){
1950: d = *(aa + ai[i]); /* diag[i] */
1951: v = aa + ai[i] + 1;
1952: vj = aj + ai[i] + 1;
1953: nz = ai[i+1] - ai[i] - 1;
1954: x[i] = omega*t[i]/d;
1955: while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */
1956: PetscLogFlops(2*nz-1);
1957: }
1958: }
1960: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
1961: for (i=0; i<m; i++)
1962: t[i] = b[i];
1963:
1964: for (i=0; i<m-1; i++){ /* update rhs */
1965: v = aa + ai[i] + 1;
1966: vj = aj + ai[i] + 1;
1967: nz = ai[i+1] - ai[i] - 1;
1968: while (nz--) t[*vj++] -= x[i]*(*v++);
1969: PetscLogFlops(2*nz-1);
1970: }
1971: for (i=m-1; i>=0; i--){
1972: d = *(aa + ai[i]);
1973: v = aa + ai[i] + 1;
1974: vj = aj + ai[i] + 1;
1975: nz = ai[i+1] - ai[i] - 1;
1976: sum = t[i];
1977: while (nz--) sum -= x[*vj++]*(*v++);
1978: PetscLogFlops(2*nz-1);
1979: x[i] = (1-omega)*x[i] + omega*sum/d;
1980: }
1981: }
1982: its--;
1983: }
1985: while (its--) {
1986: /*
1987: forward sweep:
1988: for i=0,...,m-1:
1989: sum[i] = (b[i] - U(i,:)x )/d[i];
1990: x[i] = (1-omega)x[i] + omega*sum[i];
1991: b = b - x[i]*U^T(i,:);
1992:
1993: */
1994: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
1995: for (i=0; i<m; i++)
1996: t[i] = b[i];
1998: for (i=0; i<m; i++){
1999: d = *(aa + ai[i]); /* diag[i] */
2000: v = aa + ai[i] + 1; v1=v;
2001: vj = aj + ai[i] + 1; vj1=vj;
2002: nz = ai[i+1] - ai[i] - 1; nz1=nz;
2003: sum = t[i];
2004: while (nz1--) sum -= (*v1++)*x[*vj1++];
2005: x[i] = (1-omega)*x[i] + omega*sum/d;
2006: while (nz--) t[*vj++] -= x[i]*(*v++);
2007: PetscLogFlops(4*nz-2);
2008: }
2009: }
2010:
2011: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
2012: /*
2013: backward sweep:
2014: b = b - x[i]*U^T(i,:), i=0,...,n-2
2015: for i=m-1,...,0:
2016: sum[i] = (b[i] - U(i,:)x )/d[i];
2017: x[i] = (1-omega)x[i] + omega*sum[i];
2018: */
2019: for (i=0; i<m; i++)
2020: t[i] = b[i];
2021:
2022: for (i=0; i<m-1; i++){ /* update rhs */
2023: v = aa + ai[i] + 1;
2024: vj = aj + ai[i] + 1;
2025: nz = ai[i+1] - ai[i] - 1;
2026: while (nz--) t[*vj++] -= x[i]*(*v++);
2027: PetscLogFlops(2*nz-1);
2028: }
2029: for (i=m-1; i>=0; i--){
2030: d = *(aa + ai[i]);
2031: v = aa + ai[i] + 1;
2032: vj = aj + ai[i] + 1;
2033: nz = ai[i+1] - ai[i] - 1;
2034: sum = t[i];
2035: while (nz--) sum -= x[*vj++]*(*v++);
2036: PetscLogFlops(2*nz-1);
2037: x[i] = (1-omega)*x[i] + omega*sum/d;
2038: }
2039: }
2040: }
2042: PetscFree(t);
2043: VecRestoreArray(xx,&x);
2044: if (bb != xx) {
2045: VecRestoreArray(bb,&b);
2046: }
2048: return(0);
2049: }