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: }