Actual source code: bdiag.c

      1: #define PETSCMAT_DLL

      3: /* Block diagonal matrix format */

     5:  #include src/mat/impls/bdiag/seq/bdiag.h
     6:  #include src/inline/ilu.h

     10: PetscErrorCode MatDestroy_SeqBDiag(Mat A)
     11: {
     12:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
     14:   PetscInt       i,bs = A->bs;

     17: #if defined(PETSC_USE_LOG)
     18:   PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D, BSize=%D, NDiag=%D",A->m,A->n,a->nz,A->bs,a->nd);
     19: #endif
     20:   if (!a->user_alloc) { /* Free the actual diagonals */
     21:     for (i=0; i<a->nd; i++) {
     22:       if (a->diag[i] > 0) {
     23:         PetscScalar *dummy = a->diagv[i] + bs*bs*a->diag[i];
     24:         PetscFree(dummy);
     25:       } else {
     26:         PetscFree(a->diagv[i]);
     27:       }
     28:     }
     29:   }
     30:   if (a->pivot) {PetscFree(a->pivot);}
     31:   PetscFree(a->diagv);
     32:   PetscFree(a->diag);
     33:   PetscFree(a->colloc);
     34:   PetscFree(a->dvalue);
     35:   if (a->solvework) {PetscFree(a->solvework);}
     36:   PetscFree(a);
     37:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBDiagSetPreallocation_C","",PETSC_NULL);
     38:   return(0);
     39: }

     43: PetscErrorCode MatAssemblyEnd_SeqBDiag(Mat A,MatAssemblyType mode)
     44: {
     45:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
     46:   PetscInt       i,k,temp,*diag = a->diag,*bdlen = a->bdlen;
     47:   PetscScalar    *dtemp,**dv = a->diagv;

     51:   if (mode == MAT_FLUSH_ASSEMBLY) return(0);

     53:   /* Sort diagonals */
     54:   for (i=0; i<a->nd; i++) {
     55:     for (k=i+1; k<a->nd; k++) {
     56:       if (diag[i] < diag[k]) {
     57:         temp     = diag[i];
     58:         diag[i]  = diag[k];
     59:         diag[k]  = temp;
     60:         temp     = bdlen[i];
     61:         bdlen[i] = bdlen[k];
     62:         bdlen[k] = temp;
     63:         dtemp    = dv[i];
     64:         dv[i]    = dv[k];
     65:         dv[k]    = dtemp;
     66:       }
     67:     }
     68:   }

     70:   /* Set location of main diagonal */
     71:   for (i=0; i<a->nd; i++) {
     72:     if (!a->diag[i]) {a->mainbd = i; break;}
     73:   }
     74:   PetscLogInfo((A,"MatAssemblyEnd_SeqBDiag:Number diagonals %D,memory used %D, block size %D\n",a->nd,a->maxnz,A->bs));
     75:   return(0);
     76: }

     80: PetscErrorCode MatSetOption_SeqBDiag(Mat A,MatOption op)
     81: {
     82:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;

     86:   switch (op) {
     87:   case MAT_NO_NEW_NONZERO_LOCATIONS:
     88:     a->nonew       = 1;
     89:     break;
     90:   case MAT_YES_NEW_NONZERO_LOCATIONS:
     91:     a->nonew       = 0;
     92:     break;
     93:   case MAT_NO_NEW_DIAGONALS:
     94:     a->nonew_diag  = 1;
     95:     break;
     96:   case MAT_YES_NEW_DIAGONALS:
     97:     a->nonew_diag  = 0;
     98:     break;
     99:   case MAT_COLUMN_ORIENTED:
    100:     a->roworiented = PETSC_FALSE;
    101:     break;
    102:   case MAT_ROW_ORIENTED:
    103:     a->roworiented = PETSC_TRUE;
    104:     break;
    105:   case MAT_ROWS_SORTED:
    106:   case MAT_ROWS_UNSORTED:
    107:   case MAT_COLUMNS_SORTED:
    108:   case MAT_COLUMNS_UNSORTED:
    109:   case MAT_IGNORE_OFF_PROC_ENTRIES:
    110:   case MAT_NEW_NONZERO_LOCATION_ERR:
    111:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
    112:   case MAT_USE_HASH_TABLE:
    113:     PetscLogInfo((A,"MatSetOption_SeqBDiag:Option ignored\n"));
    114:     break;
    115:   case MAT_SYMMETRIC:
    116:   case MAT_STRUCTURALLY_SYMMETRIC:
    117:   case MAT_NOT_SYMMETRIC:
    118:   case MAT_NOT_STRUCTURALLY_SYMMETRIC:
    119:   case MAT_HERMITIAN:
    120:   case MAT_NOT_HERMITIAN:
    121:   case MAT_SYMMETRY_ETERNAL:
    122:   case MAT_NOT_SYMMETRY_ETERNAL:
    123:     break;
    124:   default:
    125:     SETERRQ(PETSC_ERR_SUP,"unknown option");
    126:   }
    127:   return(0);
    128: }

    132: PetscErrorCode MatPrintHelp_SeqBDiag(Mat A)
    133: {
    134:   static PetscTruth called = PETSC_FALSE;
    135:   MPI_Comm          comm = A->comm;
    136:   PetscErrorCode    ierr;

    139:   if (called) {return(0);} else called = PETSC_TRUE;
    140:   (*PetscHelpPrintf)(comm," Options for MATSEQBDIAG and MATMPIBDIAG matrix formats:\n");
    141:   (*PetscHelpPrintf)(comm,"  -mat_block_size <block_size>\n");
    142:   (*PetscHelpPrintf)(comm,"  -mat_bdiag_diags <d1,d2,d3,...> (diagonal numbers)\n");
    143:   (*PetscHelpPrintf)(comm,"   (for example) -mat_bdiag_diags -5,-1,0,1,5\n");
    144:   return(0);
    145: }

    149: static PetscErrorCode MatGetDiagonal_SeqBDiag_N(Mat A,Vec v)
    150: {
    151:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    153:   PetscInt       i,j,n,len,ibase,bs = A->bs,iloc;
    154:   PetscScalar    *x,*dd,zero = 0.0;

    157:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    158:   VecSet(v,zero);
    159:   VecGetLocalSize(v,&n);
    160:   if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
    161:   if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal not set");
    162:   len = PetscMin(a->mblock,a->nblock);
    163:   dd = a->diagv[a->mainbd];
    164:   VecGetArray(v,&x);
    165:   for (i=0; i<len; i++) {
    166:     ibase = i*bs*bs;  iloc = i*bs;
    167:     for (j=0; j<bs; j++) x[j + iloc] = dd[ibase + j*(bs+1)];
    168:   }
    169:   VecRestoreArray(v,&x);
    170:   return(0);
    171: }

    175: static PetscErrorCode MatGetDiagonal_SeqBDiag_1(Mat A,Vec v)
    176: {
    177:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    179:   PetscInt       i,n,len;
    180:   PetscScalar    *x,*dd,zero = 0.0;

    183:   VecSet(v,zero);
    184:   VecGetLocalSize(v,&n);
    185:   if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
    186:   if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal not set");
    187:   dd = a->diagv[a->mainbd];
    188:   len = PetscMin(A->m,A->n);
    189:   VecGetArray(v,&x);
    190:   for (i=0; i<len; i++) x[i] = dd[i];
    191:   VecRestoreArray(v,&x);
    192:   return(0);
    193: }

    197: PetscErrorCode MatZeroEntries_SeqBDiag(Mat A)
    198: {
    199:   Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
    200:   PetscInt     d,i,len,bs = A->bs;
    201:   PetscScalar  *dv;

    204:   for (d=0; d<a->nd; d++) {
    205:     dv  = a->diagv[d];
    206:     if (a->diag[d] > 0) {
    207:       dv += bs*bs*a->diag[d];
    208:     }
    209:     len = a->bdlen[d]*bs*bs;
    210:     for (i=0; i<len; i++) dv[i] = 0.0;
    211:   }
    212:   return(0);
    213: }

    217: PetscErrorCode MatZeroRows_SeqBDiag(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag)
    218: {
    219:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    221:   PetscInt       i,m = A->m - 1,nz;
    222:   PetscScalar    *dd;
    223:   PetscScalar    *val;

    226:   for (i=0; i<N; i++) {
    227:     if (rows[i]<0 || rows[i]>m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range");
    228:     MatGetRow_SeqBDiag(A,rows[i],&nz,PETSC_NULL,&val);
    229:     PetscMemzero((void*)val,nz*sizeof(PetscScalar));
    230:     MatRestoreRow_SeqBDiag(A,rows[i],&nz,PETSC_NULL,&val);
    231:   }
    232:   if (diag != 0.0) {
    233:     if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal does not exist");
    234:     dd = a->diagv[a->mainbd];
    235:     for (i=0; i<N; i++) dd[rows[i]] = diag;
    236:   }
    237:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    238:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    239:   return(0);
    240: }

    244: PetscErrorCode MatGetSubMatrix_SeqBDiag(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *submat)
    245: {
    247:   PetscInt       nznew,*smap,i,j,oldcols = A->n;
    248:   PetscInt       *irow,*icol,newr,newc,*cwork,nz,bs;
    249:   PetscInt       *col;
    250:   PetscScalar    *vwork;
    251:   PetscScalar    *val;
    252:   Mat            newmat;

    255:   if (scall == MAT_REUSE_MATRIX) { /* no support for reuse so simply destroy all */
    256:     MatDestroy(*submat);
    257:   }

    259:   ISGetIndices(isrow,&irow);
    260:   ISGetIndices(iscol,&icol);
    261:   ISGetLocalSize(isrow,&newr);
    262:   ISGetLocalSize(iscol,&newc);

    264:   PetscMalloc((oldcols+1)*sizeof(PetscInt),&smap);
    265:   PetscMalloc((newc+1)*sizeof(PetscInt),&cwork);
    266:   PetscMalloc((newc+1)*sizeof(PetscScalar),&vwork);
    267:   PetscMemzero((char*)smap,oldcols*sizeof(PetscInt));
    268:   for (i=0; i<newc; i++) smap[icol[i]] = i+1;

    270:   /* Determine diagonals; then create submatrix */
    271:   bs = A->bs; /* Default block size remains the same */
    272:   MatCreate(A->comm,&newmat);
    273:   MatSetSizes(newmat,newr,newc,newr,newc);
    274:   MatSetType(newmat,A->type_name);
    275:   MatSeqBDiagSetPreallocation(newmat,0,bs,PETSC_NULL,PETSC_NULL);

    277:   /* Fill new matrix */
    278:   for (i=0; i<newr; i++) {
    279:     MatGetRow_SeqBDiag(A,irow[i],&nz,&col,&val);
    280:     nznew = 0;
    281:     for (j=0; j<nz; j++) {
    282:       if (smap[col[j]]) {
    283:         cwork[nznew]   = smap[col[j]] - 1;
    284:         vwork[nznew++] = val[j];
    285:       }
    286:     }
    287:     MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES);
    288:     MatRestoreRow_SeqBDiag(A,i,&nz,&col,&val);
    289:   }
    290:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
    291:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

    293:   /* Free work space */
    294:   PetscFree(smap);
    295:   PetscFree(cwork);
    296:   PetscFree(vwork);
    297:   ISRestoreIndices(isrow,&irow);
    298:   ISRestoreIndices(iscol,&icol);
    299:   *submat = newmat;
    300:   return(0);
    301: }

    305: PetscErrorCode MatGetSubMatrices_SeqBDiag(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
    306: {
    308:   PetscInt       i;

    311:   if (scall == MAT_INITIAL_MATRIX) {
    312:     PetscMalloc((n+1)*sizeof(Mat),B);
    313:   }

    315:   for (i=0; i<n; i++) {
    316:     MatGetSubMatrix_SeqBDiag(A,irow[i],icol[i],scall,&(*B)[i]);
    317:   }
    318:   return(0);
    319: }

    323: PetscErrorCode MatScale_SeqBDiag(Mat inA,PetscScalar alpha)
    324: {
    325:   Mat_SeqBDiag *a = (Mat_SeqBDiag*)inA->data;
    326:   PetscInt          i,bs = inA->bs;
    327:   PetscScalar  oalpha = alpha;
    328:   PetscBLASInt one = 1,len;

    332:   for (i=0; i<a->nd; i++) {
    333:     len = (PetscBLASInt)bs*bs*a->bdlen[i];
    334:     if (a->diag[i] > 0) {
    335:       BLASscal_(&len,&oalpha,a->diagv[i] + bs*bs*a->diag[i],&one);
    336:     } else {
    337:       BLASscal_(&len,&oalpha,a->diagv[i],&one);
    338:     }
    339:   }
    340:   PetscLogFlops(a->nz);
    341:   return(0);
    342: }

    346: PetscErrorCode MatDiagonalScale_SeqBDiag(Mat A,Vec ll,Vec rr)
    347: {
    348:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    349:   PetscScalar    *l,*r,*dv;
    351:   PetscInt       d,j,len;
    352:   PetscInt       nd = a->nd,bs = A->bs,diag,m,n;

    355:   if (ll) {
    356:     VecGetSize(ll,&m);
    357:     if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length");
    358:     if (bs == 1) {
    359:       VecGetArray(ll,&l);
    360:       for (d=0; d<nd; d++) {
    361:         dv   = a->diagv[d];
    362:         diag = a->diag[d];
    363:         len  = a->bdlen[d];
    364:         if (diag > 0) for (j=0; j<len; j++) dv[j+diag] *= l[j+diag];
    365:         else          for (j=0; j<len; j++) dv[j]      *= l[j];
    366:       }
    367:       VecRestoreArray(ll,&l);
    368:       PetscLogFlops(a->nz);
    369:     } else SETERRQ(PETSC_ERR_SUP,"Not yet done for bs>1");
    370:   }
    371:   if (rr) {
    372:     VecGetSize(rr,&n);
    373:     if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length");
    374:     if (bs == 1) {
    375:       VecGetArray(rr,&r);
    376:       for (d=0; d<nd; d++) {
    377:         dv   = a->diagv[d];
    378:         diag = a->diag[d];
    379:         len  = a->bdlen[d];
    380:         if (diag > 0) for (j=0; j<len; j++) dv[j+diag] *= r[j];
    381:         else          for (j=0; j<len; j++) dv[j]      *= r[j-diag];
    382:       }
    383:       VecRestoreArray(rr,&r);
    384:       PetscLogFlops(a->nz);
    385:     } else SETERRQ(PETSC_ERR_SUP,"Not yet done for bs>1");
    386:   }
    387:   return(0);
    388: }

    390: static PetscErrorCode MatDuplicate_SeqBDiag(Mat,MatDuplicateOption,Mat *);

    394: PetscErrorCode MatSetUpPreallocation_SeqBDiag(Mat A)
    395: {

    399:    MatSeqBDiagSetPreallocation(A,PETSC_DEFAULT,PETSC_DEFAULT,0,0);
    400:   return(0);
    401: }

    403: /* -------------------------------------------------------------------*/
    404: static struct _MatOps MatOps_Values = {MatSetValues_SeqBDiag_N,
    405:        MatGetRow_SeqBDiag,
    406:        MatRestoreRow_SeqBDiag,
    407:        MatMult_SeqBDiag_N,
    408: /* 4*/ MatMultAdd_SeqBDiag_N,
    409:        MatMultTranspose_SeqBDiag_N,
    410:        MatMultTransposeAdd_SeqBDiag_N,
    411:        MatSolve_SeqBDiag_N,
    412:        0,
    413:        0,
    414: /*10*/ 0,
    415:        0,
    416:        0,
    417:        MatRelax_SeqBDiag_N,
    418:        MatTranspose_SeqBDiag,
    419: /*15*/ MatGetInfo_SeqBDiag,
    420:        0,
    421:        MatGetDiagonal_SeqBDiag_N,
    422:        MatDiagonalScale_SeqBDiag,
    423:        MatNorm_SeqBDiag,
    424: /*20*/ 0,
    425:        MatAssemblyEnd_SeqBDiag,
    426:        0,
    427:        MatSetOption_SeqBDiag,
    428:        MatZeroEntries_SeqBDiag,
    429: /*25*/ MatZeroRows_SeqBDiag,
    430:        0,
    431:        MatLUFactorNumeric_SeqBDiag_N,
    432:        0,
    433:        0,
    434: /*30*/ MatSetUpPreallocation_SeqBDiag,
    435:        MatILUFactorSymbolic_SeqBDiag,
    436:        0,
    437:        0,
    438:        0,
    439: /*35*/ MatDuplicate_SeqBDiag,
    440:        0,
    441:        0,
    442:        MatILUFactor_SeqBDiag,
    443:        0,
    444: /*40*/ 0,
    445:        MatGetSubMatrices_SeqBDiag,
    446:        0,
    447:        MatGetValues_SeqBDiag_N,
    448:        0,
    449: /*45*/ MatPrintHelp_SeqBDiag,
    450:        MatScale_SeqBDiag,
    451:        0,
    452:        0,
    453:        0,
    454: /*50*/ 0,
    455:        0,
    456:        0,
    457:        0,
    458:        0,
    459: /*55*/ 0,
    460:        0,
    461:        0,
    462:        0,
    463:        0,
    464: /*60*/ 0,
    465:        MatDestroy_SeqBDiag,
    466:        MatView_SeqBDiag,
    467:        MatGetPetscMaps_Petsc,
    468:        0,
    469: /*65*/ 0,
    470:        0,
    471:        0,
    472:        0,
    473:        0,
    474: /*70*/ 0,
    475:        0,
    476:        0,
    477:        0,
    478:        0,
    479: /*75*/ 0,
    480:        0,
    481:        0,
    482:        0,
    483:        0,
    484: /*80*/ 0,
    485:        0,
    486:        0,
    487:        0,
    488:        MatLoad_SeqBDiag,
    489: /*85*/ 0,
    490:        0,
    491:        0,
    492:        0,
    493:        0,
    494: /*90*/ 0,
    495:        0,
    496:        0,
    497:        0,
    498:        0,
    499: /*95*/ 0,
    500:        0,
    501:        0,
    502:        0};

    506: /*@C
    507:    MatSeqBDiagSetPreallocation - Sets the nonzero structure and (optionally) arrays.

    509:    Collective on MPI_Comm

    511:    Input Parameters:
    512: +  B - the matrix
    513: .  nd - number of block diagonals (optional)
    514: .  bs - each element of a diagonal is an bs x bs dense matrix
    515: .  diag - optional array of block diagonal numbers (length nd).
    516:    For a matrix element A[i,j], where i=row and j=column, the
    517:    diagonal number is
    518: $     diag = i/bs - j/bs  (integer division)
    519:    Set diag=PETSC_NULL on input for PETSc to dynamically allocate memory as
    520:    needed (expensive).
    521: -  diagv - pointer to actual diagonals (in same order as diag array),
    522:    if allocated by user.  Otherwise, set diagv=PETSC_NULL on input for PETSc
    523:    to control memory allocation.

    525:    Options Database Keys:
    526: .  -mat_block_size <bs> - Sets blocksize
    527: .  -mat_bdiag_diags <s1,s2,s3,...> - Sets diagonal numbers

    529:    Notes:
    530:    See the users manual for further details regarding this storage format.

    532:    Fortran Note:
    533:    Fortran programmers cannot set diagv; this value is ignored.

    535:    Level: intermediate

    537: .keywords: matrix, block, diagonal, sparse

    539: .seealso: MatCreate(), MatCreateMPIBDiag(), MatSetValues()
    540: @*/
    541: PetscErrorCode PETSCMAT_DLLEXPORT MatSeqBDiagSetPreallocation(Mat B,PetscInt nd,PetscInt bs,const PetscInt diag[],PetscScalar *diagv[])
    542: {
    543:   PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[],PetscScalar*[]);

    546:   PetscObjectQueryFunction((PetscObject)B,"MatSeqBDiagSetPreallocation_C",(void (**)(void))&f);
    547:   if (f) {
    548:     (*f)(B,nd,bs,diag,diagv);
    549:   }
    550:   return(0);
    551: }

    556: PetscErrorCode PETSCMAT_DLLEXPORT MatSeqBDiagSetPreallocation_SeqBDiag(Mat B,PetscInt nd,PetscInt bs,PetscInt *diag,PetscScalar **diagv)
    557: {
    558:   Mat_SeqBDiag   *b;
    560:   PetscInt       i,nda,sizetot, nd2 = 128,idiag[128];
    561:   PetscTruth     flg1;

    565:   B->preallocated = PETSC_TRUE;
    566:   if (bs == PETSC_DEFAULT) bs = 1;
    567:   if (!bs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize cannot be zero");
    568:   if (nd == PETSC_DEFAULT) nd = 0;
    569:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
    570:   PetscOptionsGetIntArray(PETSC_NULL,"-mat_bdiag_diags",idiag,&nd2,&flg1);
    571:   if (flg1) {
    572:     diag = idiag;
    573:     nd   = nd2;
    574:   }

    576:   if ((B->n%bs) || (B->m%bs)) SETERRQ(PETSC_ERR_ARG_SIZ,"Invalid block size");
    577:   if (!nd) nda = nd + 1;
    578:   else     nda = nd;
    579:   b            = (Mat_SeqBDiag*)B->data;

    581:   PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg1);
    582:   if (!flg1) {
    583:     switch (bs) {
    584:       case 1:
    585:         B->ops->setvalues       = MatSetValues_SeqBDiag_1;
    586:         B->ops->getvalues       = MatGetValues_SeqBDiag_1;
    587:         B->ops->getdiagonal     = MatGetDiagonal_SeqBDiag_1;
    588:         B->ops->mult            = MatMult_SeqBDiag_1;
    589:         B->ops->multadd         = MatMultAdd_SeqBDiag_1;
    590:         B->ops->multtranspose   = MatMultTranspose_SeqBDiag_1;
    591:         B->ops->multtransposeadd= MatMultTransposeAdd_SeqBDiag_1;
    592:         B->ops->relax           = MatRelax_SeqBDiag_1;
    593:         B->ops->solve           = MatSolve_SeqBDiag_1;
    594:         B->ops->lufactornumeric = MatLUFactorNumeric_SeqBDiag_1;
    595:         break;
    596:       case 2:
    597:         B->ops->mult            = MatMult_SeqBDiag_2;
    598:         B->ops->multadd         = MatMultAdd_SeqBDiag_2;
    599:         B->ops->solve           = MatSolve_SeqBDiag_2;
    600:         break;
    601:       case 3:
    602:         B->ops->mult            = MatMult_SeqBDiag_3;
    603:         B->ops->multadd         = MatMultAdd_SeqBDiag_3;
    604:         B->ops->solve           = MatSolve_SeqBDiag_3;
    605:         break;
    606:       case 4:
    607:         B->ops->mult            = MatMult_SeqBDiag_4;
    608:         B->ops->multadd         = MatMultAdd_SeqBDiag_4;
    609:         B->ops->solve           = MatSolve_SeqBDiag_4;
    610:         break;
    611:       case 5:
    612:         B->ops->mult            = MatMult_SeqBDiag_5;
    613:         B->ops->multadd         = MatMultAdd_SeqBDiag_5;
    614:         B->ops->solve           = MatSolve_SeqBDiag_5;
    615:         break;
    616:    }
    617:   }

    619:   b->mblock = B->m/bs;
    620:   b->nblock = B->n/bs;
    621:   b->nd     = nd;
    622:   B->bs     = bs;
    623:   b->ndim   = 0;
    624:   b->mainbd = -1;
    625:   b->pivot  = 0;

    627:   PetscMalloc(2*nda*sizeof(PetscInt),&b->diag);
    628:   b->bdlen  = b->diag + nda;
    629:   PetscMalloc((B->n+1)*sizeof(PetscInt),&b->colloc);
    630:   PetscMalloc(nda*sizeof(PetscScalar*),&b->diagv);
    631:   sizetot   = 0;

    633:   if (diagv) { /* user allocated space */
    634:     b->user_alloc = PETSC_TRUE;
    635:     for (i=0; i<nd; i++) b->diagv[i] = diagv[i];
    636:   } else b->user_alloc = PETSC_FALSE;

    638:   for (i=0; i<nd; i++) {
    639:     b->diag[i] = diag[i];
    640:     if (diag[i] > 0) { /* lower triangular */
    641:       b->bdlen[i] = PetscMin(b->nblock,b->mblock - diag[i]);
    642:     } else {           /* upper triangular */
    643:       b->bdlen[i] = PetscMin(b->mblock,b->nblock + diag[i]);
    644:     }
    645:     sizetot += b->bdlen[i];
    646:   }
    647:   sizetot   *= bs*bs;
    648:   b->maxnz  =  sizetot;
    649:   PetscMalloc((B->n+1)*sizeof(PetscScalar),&b->dvalue);
    650:   PetscLogObjectMemory(B,(nda*(bs+2))*sizeof(PetscInt) + bs*nda*sizeof(PetscScalar)
    651:                     + nda*sizeof(PetscScalar*) + sizeof(Mat_SeqBDiag)
    652:                     + sizeof(struct _p_Mat) + sizetot*sizeof(PetscScalar));

    654:   if (!b->user_alloc) {
    655:     for (i=0; i<nd; i++) {
    656:       PetscMalloc(bs*bs*b->bdlen[i]*sizeof(PetscScalar),&b->diagv[i]);
    657:       PetscMemzero(b->diagv[i],bs*bs*b->bdlen[i]*sizeof(PetscScalar));
    658:     }
    659:     b->nonew = 0; b->nonew_diag = 0;
    660:   } else { /* diagonals are set on input; don't allow dynamic allocation */
    661:     b->nonew = 1; b->nonew_diag = 1;
    662:   }

    664:   /* adjust diagv so one may access rows with diagv[diag][row] for all rows */
    665:   for (i=0; i<nd; i++) {
    666:     if (diag[i] > 0) {
    667:       b->diagv[i] -= bs*bs*diag[i];
    668:     }
    669:   }

    671:   b->nz          = b->maxnz; /* Currently not keeping track of exact count */
    672:   b->roworiented = PETSC_TRUE;
    673:   B->info.nz_unneeded = (double)b->maxnz;
    674:   return(0);
    675: }

    680: static PetscErrorCode MatDuplicate_SeqBDiag(Mat A,MatDuplicateOption cpvalues,Mat *matout)
    681: {
    682:   Mat_SeqBDiag   *newmat,*a = (Mat_SeqBDiag*)A->data;
    684:   PetscInt       i,len,diag,bs = A->bs;
    685:   Mat            mat;

    688:   MatCreate(A->comm,matout);
    689:   MatSetSizes(*matout,A->m,A->n,A->m,A->n);
    690:   MatSetType(*matout,A->type_name);
    691:   MatSeqBDiagSetPreallocation(*matout,a->nd,bs,a->diag,PETSC_NULL);

    693:   /* Copy contents of diagonals */
    694:   mat = *matout;
    695:   newmat = (Mat_SeqBDiag*)mat->data;
    696:   if (cpvalues == MAT_COPY_VALUES) {
    697:     for (i=0; i<a->nd; i++) {
    698:       len = a->bdlen[i] * bs * bs * sizeof(PetscScalar);
    699:       diag = a->diag[i];
    700:       if (diag > 0) {
    701:         PetscMemcpy(newmat->diagv[i]+bs*bs*diag,a->diagv[i]+bs*bs*diag,len);
    702:       } else {
    703:         PetscMemcpy(newmat->diagv[i],a->diagv[i],len);
    704:       }
    705:     }
    706:   }
    707:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    708:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    709:   return(0);
    710: }

    714: PetscErrorCode MatLoad_SeqBDiag(PetscViewer viewer, MatType type,Mat *A)
    715: {
    716:   Mat            B;
    718:   PetscMPIInt    size;
    719:   int            fd;
    720:   PetscInt       *scols,i,nz,header[4],nd = 128;
    721:   PetscInt       bs,*rowlengths = 0,M,N,*cols,extra_rows,*diag = 0;
    722:   PetscInt       idiag[128];
    723:   PetscScalar    *vals,*svals;
    724:   MPI_Comm       comm;
    725:   PetscTruth     flg;
    726:
    728:   PetscObjectGetComm((PetscObject)viewer,&comm);
    729:   MPI_Comm_size(comm,&size);
    730:   if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor");
    731:   PetscViewerBinaryGetDescriptor(viewer,&fd);
    732:   PetscBinaryRead(fd,header,4,PETSC_INT);
    733:   if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
    734:   M = header[1]; N = header[2]; nz = header[3];
    735:   if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only load square matrices");
    736:   if (header[3] < 0) {
    737:     SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBDiag");
    738:   }

    740:   /*
    741:      This code adds extra rows to make sure the number of rows is
    742:     divisible by the blocksize
    743:   */
    744:   bs = 1;
    745:   PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);
    746:   extra_rows = bs - M + bs*(M/bs);
    747:   if (extra_rows == bs) extra_rows = 0;
    748:   if (extra_rows) {
    749:     PetscLogInfo((0,"MatLoad_SeqBDiag:Padding loaded matrix to match blocksize\n"));
    750:   }

    752:   /* read row lengths */
    753:   PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);
    754:   PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
    755:   for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;

    757:   /* load information about diagonals */
    758:   PetscOptionsGetIntArray(PETSC_NULL,"-matload_bdiag_diags",idiag,&nd,&flg);
    759:   if (flg) {
    760:     diag = idiag;
    761:   }

    763:   /* create our matrix */
    764:   MatCreate(comm,A);
    765:   MatSetSizes(*A,M+extra_rows,M+extra_rows,M+extra_rows,M+extra_rows);
    766:   MatSetType(*A,type);
    767:   MatSeqBDiagSetPreallocation(*A,nd,bs,diag,PETSC_NULL);
    768:   B = *A;

    770:   /* read column indices and nonzeros */
    771:   PetscMalloc(nz*sizeof(PetscInt),&scols);
    772:   cols = scols;
    773:   PetscBinaryRead(fd,cols,nz,PETSC_INT);
    774:   PetscMalloc(nz*sizeof(PetscScalar),&svals);
    775:   vals = svals;
    776:   PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
    777:   /* insert into matrix */

    779:   for (i=0; i<M; i++) {
    780:     MatSetValues(B,1,&i,rowlengths[i],scols,svals,INSERT_VALUES);
    781:     scols += rowlengths[i]; svals += rowlengths[i];
    782:   }
    783:   vals[0] = 1.0;
    784:   for (i=M; i<M+extra_rows; i++) {
    785:     MatSetValues(B,1,&i,1,&i,vals,INSERT_VALUES);
    786:   }

    788:   PetscFree(cols);
    789:   PetscFree(vals);
    790:   PetscFree(rowlengths);

    792:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    793:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    794:   return(0);
    795: }

    797: /*MC
    798:    MATSEQBDIAG - MATSEQBDIAG = "seqbdiag" - A matrix type to be used for sequential block diagonal matrices.

    800:    Options Database Keys:
    801: . -mat_type seqbdiag - sets the matrix type to "seqbdiag" during a call to MatSetFromOptions()

    803:   Level: beginner

    805: .seealso: MatCreateSeqBDiag
    806: M*/

    811: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqBDiag(Mat B)
    812: {
    813:   Mat_SeqBDiag   *b;
    815:   PetscMPIInt    size;

    818:   MPI_Comm_size(B->comm,&size);
    819:   if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");

    821:   B->m = B->M = PetscMax(B->m,B->M);
    822:   B->n = B->N = PetscMax(B->n,B->N);

    824:   PetscNew(Mat_SeqBDiag,&b);
    825:   B->data         = (void*)b;
    826:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
    827:   B->factor       = 0;
    828:   B->mapping      = 0;

    830:   PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);
    831:   PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);

    833:   b->ndim   = 0;
    834:   b->mainbd = -1;
    835:   b->pivot  = 0;

    837:   b->roworiented = PETSC_TRUE;
    838:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqBDiagSetPreallocation_C",
    839:                                     "MatSeqBDiagSetPreallocation_SeqBDiag",
    840:                                      MatSeqBDiagSetPreallocation_SeqBDiag);

    842:   return(0);
    843: }

    848: /*@C
    849:    MatCreateSeqBDiag - Creates a sequential block diagonal matrix.

    851:    Collective on MPI_Comm

    853:    Input Parameters:
    854: +  comm - MPI communicator, set to PETSC_COMM_SELF
    855: .  m - number of rows
    856: .  n - number of columns
    857: .  nd - number of block diagonals (optional)
    858: .  bs - each element of a diagonal is an bs x bs dense matrix
    859: .  diag - optional array of block diagonal numbers (length nd).
    860:    For a matrix element A[i,j], where i=row and j=column, the
    861:    diagonal number is
    862: $     diag = i/bs - j/bs  (integer division)
    863:    Set diag=PETSC_NULL on input for PETSc to dynamically allocate memory as
    864:    needed (expensive).
    865: -  diagv - pointer to actual diagonals (in same order as diag array),
    866:    if allocated by user.  Otherwise, set diagv=PETSC_NULL on input for PETSc
    867:    to control memory allocation.

    869:    Output Parameters:
    870: .  A - the matrix

    872:    Options Database Keys:
    873: .  -mat_block_size <bs> - Sets blocksize
    874: .  -mat_bdiag_diags <s1,s2,s3,...> - Sets diagonal numbers

    876:    Notes:
    877:    See the users manual for further details regarding this storage format.

    879:    Fortran Note:
    880:    Fortran programmers cannot set diagv; this value is ignored.

    882:    Level: intermediate

    884: .keywords: matrix, block, diagonal, sparse

    886: .seealso: MatCreate(), MatCreateMPIBDiag(), MatSetValues()
    887: @*/
    888: PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqBDiag(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nd,PetscInt bs,const PetscInt diag[],PetscScalar *diagv[],Mat *A)
    889: {

    893:   MatCreate(comm,A);
    894:   MatSetSizes(*A,m,n,m,n);
    895:   MatSetType(*A,MATSEQBDIAG);
    896:   MatSeqBDiagSetPreallocation(*A,nd,bs,diag,diagv);
    897:   return(0);
    898: }