Actual source code: bdiag3.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
     7:  #include petscsys.h

     11: PetscErrorCode MatGetInfo_SeqBDiag(Mat A,MatInfoType flag,MatInfo *info)
     12: {
     13:   Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;

     16:   info->rows_global       = (double)A->m;
     17:   info->columns_global    = (double)A->n;
     18:   info->rows_local        = (double)A->m;
     19:   info->columns_local     = (double)A->n;
     20:   info->block_size        = A->bs;
     21:   info->nz_allocated      = (double)a->maxnz;
     22:   info->nz_used           = (double)a->nz;
     23:   info->nz_unneeded       = (double)(a->maxnz - a->nz);
     24:   info->assemblies        = (double)A->num_ass;
     25:   info->mallocs           = (double)a->reallocs;
     26:   info->memory            = A->mem;
     27:   info->fill_ratio_given  = 0; /* supports ILU(0) only */
     28:   info->fill_ratio_needed = 0;
     29:   info->factor_mallocs    = 0;
     30:   return(0);
     31: }

     33: /*
     34:      Note: this currently will generate different answers if you request
     35:  all items or a subset. If you request all items it checks if the value is
     36:  nonzero and only includes it if it is nonzero; if you check a subset of items
     37:  it returns a list of all active columns in the row (some which may contain
     38:  a zero)
     39: */
     42: PetscErrorCode MatGetRow_SeqBDiag(Mat A,PetscInt row,PetscInt *nz,PetscInt **col,PetscScalar **v)
     43: {
     44:   Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
     45:   PetscInt     nd = a->nd,bs = A->bs;
     46:   PetscInt     nc = A->n,*diag = a->diag,pcol,shift,i,j,k;

     49:   /* For efficiency, if ((nz) && (col) && (v)) then do all at once */
     50:   if ((nz) && (col) && (v)) {
     51:     *col = a->colloc;
     52:     *v   = a->dvalue;
     53:     k    = 0;
     54:     if (bs == 1) {
     55:       for (j=0; j<nd; j++) {
     56:         pcol = row - diag[j];
     57: #if defined(PETSC_USE_COMPLEX)
     58:         if (pcol > -1 && pcol < nc && PetscAbsScalar((a->diagv[j])[row])) {
     59: #else
     60:         if (pcol > -1 && pcol < nc && (a->diagv[j])[row]) {
     61: #endif
     62:           (*v)[k]   = (a->diagv[j])[row];
     63:           (*col)[k] = pcol;
     64:           k++;
     65:         }
     66:       }
     67:       *nz = k;
     68:     } else {
     69:       shift = (row/bs)*bs*bs + row%bs;
     70:       for (j=0; j<nd; j++) {
     71:         pcol = bs * (row/bs - diag[j]);
     72:         if (pcol > -1 && pcol < nc) {
     73:           for (i=0; i<bs; i++) {
     74: #if defined(PETSC_USE_COMPLEX)
     75:             if (PetscAbsScalar((a->diagv[j])[shift + i*bs])) {
     76: #else
     77:             if ((a->diagv[j])[shift + i*bs]) {
     78: #endif
     79:               (*v)[k]   = (a->diagv[j])[shift + i*bs];
     80:               (*col)[k] = pcol + i;
     81:               k++;
     82:             }
     83:           }
     84:         }
     85:       }
     86:       *nz = k;
     87:     }
     88:   } else {
     89:     if (bs == 1) {
     90:       if (nz) {
     91:         k = 0;
     92:         for (j=0; j<nd; j++) {
     93:           pcol = row - diag[j];
     94:           if (pcol > -1 && pcol < nc) k++;
     95:         }
     96:         *nz = k;
     97:       }
     98:       if (col) {
     99:         *col = a->colloc;
    100:         k = 0;
    101:         for (j=0; j<nd; j++) {
    102:           pcol = row - diag[j];
    103:           if (pcol > -1 && pcol < nc) {
    104:             (*col)[k] = pcol;  k++;
    105:           }
    106:         }
    107:       }
    108:       if (v) {
    109:         *v = a->dvalue;
    110:         k = 0;
    111:         for (j=0; j<nd; j++) {
    112:           pcol = row - diag[j];
    113:           if (pcol > -1 && pcol < nc) {
    114:             (*v)[k] = (a->diagv[j])[row]; k++;
    115:           }
    116:         }
    117:       }
    118:     } else {
    119:       if (nz) {
    120:         k = 0;
    121:         for (j=0; j<nd; j++) {
    122:           pcol = bs * (row/bs- diag[j]);
    123:           if (pcol > -1 && pcol < nc) k += bs;
    124:         }
    125:         *nz = k;
    126:       }
    127:       if (col) {
    128:         *col = a->colloc;
    129:         k = 0;
    130:         for (j=0; j<nd; j++) {
    131:           pcol = bs * (row/bs - diag[j]);
    132:           if (pcol > -1 && pcol < nc) {
    133:             for (i=0; i<bs; i++) {
    134:               (*col)[k+i] = pcol + i;
    135:             }
    136:             k += bs;
    137:           }
    138:         }
    139:       }
    140:       if (v) {
    141:         shift = (row/bs)*bs*bs + row%bs;
    142:         *v = a->dvalue;
    143:         k = 0;
    144:         for (j=0; j<nd; j++) {
    145:           pcol = bs * (row/bs - diag[j]);
    146:           if (pcol > -1 && pcol < nc) {
    147:             for (i=0; i<bs; i++) {
    148:              (*v)[k+i] = (a->diagv[j])[shift + i*bs];
    149:             }
    150:             k += bs;
    151:           }
    152:         }
    153:       }
    154:     }
    155:   }
    156:   return(0);
    157: }

    161: PetscErrorCode MatRestoreRow_SeqBDiag(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
    162: {
    164:   /* Work space is allocated during matrix creation and freed
    165:      when matrix is destroyed */
    166:   return(0);
    167: }

    169: /*
    170:    MatNorm_SeqBDiag_Columns - Computes the column norms of a block diagonal
    171:    matrix.  We code this separately from MatNorm_SeqBDiag() so that the
    172:    routine can be used for the parallel version as well.
    173:  */
    176: PetscErrorCode MatNorm_SeqBDiag_Columns(Mat A,PetscReal *tmp,PetscInt n)
    177: {
    178:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    180:   PetscInt       d,i,j,k,nd = a->nd,bs = A->bs,diag,kshift,kloc,len;
    181:   PetscScalar    *dv;

    184:   PetscMemzero(tmp,A->n*sizeof(PetscReal));
    185:   if (bs == 1) {
    186:     for (d=0; d<nd; d++) {
    187:       dv   = a->diagv[d];
    188:       diag = a->diag[d];
    189:       len  = a->bdlen[d];
    190:       if (diag > 0) {        /* lower triangle */
    191:         for (i=0; i<len; i++) {
    192:           tmp[i] += PetscAbsScalar(dv[i+diag]);
    193:         }
    194:       } else {        /* upper triangle */
    195:         for (i=0; i<len; i++) {
    196:           tmp[i-diag] += PetscAbsScalar(dv[i]);
    197:         }
    198:       }
    199:     }
    200:   } else {
    201:     for (d=0; d<nd; d++) {
    202:       dv   = a->diagv[d];
    203:       diag = a->diag[d];
    204:       len  = a->bdlen[d];

    206:       if (diag > 0) {        /* lower triangle */
    207:         for (k=0; k<len; k++) {
    208:           kloc = k*bs; kshift = kloc*bs + diag*bs;
    209:           for (i=0; i<bs; i++) {        /* i = local row */
    210:             for (j=0; j<bs; j++) {        /* j = local column */
    211:               tmp[kloc + j] += PetscAbsScalar(dv[kshift + j*bs + i]);
    212:             }
    213:           }
    214:         }
    215:       } else {        /* upper triangle */
    216:         for (k=0; k<len; k++) {
    217:           kloc = k*bs; kshift = kloc*bs;
    218:           for (i=0; i<bs; i++) {        /* i = local row */
    219:             for (j=0; j<bs; j++) {        /* j = local column */
    220:               tmp[kloc + j - bs*diag] += PetscAbsScalar(dv[kshift + j*bs + i]);
    221:             }
    222:           }
    223:         }
    224:       }
    225:     }
    226:   }
    227:   return(0);
    228: }

    232: PetscErrorCode MatNorm_SeqBDiag(Mat A,NormType type,PetscReal *nrm)
    233: {
    234:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    235:   PetscReal      sum = 0.0,*tmp;
    237:   PetscInt       d,i,j,k,nd = a->nd,bs = A->bs,diag,kshift,kloc,len;
    238:   PetscScalar    *dv;

    241:   if (type == NORM_FROBENIUS) {
    242:     for (d=0; d<nd; d++) {
    243:       dv   = a->diagv[d];
    244:       len  = a->bdlen[d]*bs*bs;
    245:       diag = a->diag[d];
    246:       if (diag > 0) {
    247:         for (i=0; i<len; i++) {
    248: #if defined(PETSC_USE_COMPLEX)
    249:           sum += PetscRealPart(PetscConj(dv[i+diag])*dv[i+diag]);
    250: #else
    251:           sum += dv[i+diag]*dv[i+diag];
    252: #endif
    253:         }
    254:       } else {
    255:         for (i=0; i<len; i++) {
    256: #if defined(PETSC_USE_COMPLEX)
    257:           sum += PetscRealPart(PetscConj(dv[i])*dv[i]);
    258: #else
    259:           sum += dv[i]*dv[i];
    260: #endif
    261:         }
    262:       }
    263:     }
    264:     *nrm = sqrt(sum);
    265:   } else if (type == NORM_1) { /* max column norm */
    266:     PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);
    267:     MatNorm_SeqBDiag_Columns(A,tmp,A->n);
    268:     *nrm = 0.0;
    269:     for (j=0; j<A->n; j++) {
    270:       if (tmp[j] > *nrm) *nrm = tmp[j];
    271:     }
    272:     PetscFree(tmp);
    273:   } else if (type == NORM_INFINITY) { /* max row norm */
    274:     PetscMalloc((A->m+1)*sizeof(PetscReal),&tmp);
    275:     PetscMemzero(tmp,A->m*sizeof(PetscReal));
    276:     *nrm = 0.0;
    277:     if (bs == 1) {
    278:       for (d=0; d<nd; d++) {
    279:         dv   = a->diagv[d];
    280:         diag = a->diag[d];
    281:         len  = a->bdlen[d];
    282:         if (diag > 0) {        /* lower triangle */
    283:           for (i=0; i<len; i++) {
    284:             tmp[i+diag] += PetscAbsScalar(dv[i+diag]);
    285:           }
    286:         } else {        /* upper triangle */
    287:           for (i=0; i<len; i++) {
    288:             tmp[i] += PetscAbsScalar(dv[i]);
    289:           }
    290:         }
    291:       }
    292:     } else {
    293:       for (d=0; d<nd; d++) {
    294:         dv   = a->diagv[d];
    295:         diag = a->diag[d];
    296:         len  = a->bdlen[d];
    297:         if (diag > 0) {
    298:           for (k=0; k<len; k++) {
    299:             kloc = k*bs; kshift = kloc*bs + bs*diag;
    300:             for (i=0; i<bs; i++) {        /* i = local row */
    301:               for (j=0; j<bs; j++) {        /* j = local column */
    302:                 tmp[kloc + i + bs*diag] += PetscAbsScalar(dv[kshift+j*bs+i]);
    303:               }
    304:             }
    305:           }
    306:         } else {
    307:           for (k=0; k<len; k++) {
    308:             kloc = k*bs; kshift = kloc*bs;
    309:             for (i=0; i<bs; i++) {        /* i = local row */
    310:               for (j=0; j<bs; j++) {        /* j = local column */
    311:                 tmp[kloc + i] += PetscAbsScalar(dv[kshift + j*bs + i]);
    312:               }
    313:             }
    314:           }
    315:         }
    316:       }
    317:     }
    318:     for (j=0; j<A->m; j++) {
    319:       if (tmp[j] > *nrm) *nrm = tmp[j];
    320:     }
    321:     PetscFree(tmp);
    322:   } else {
    323:     SETERRQ(PETSC_ERR_SUP,"No support for two norm");
    324:   }
    325:   return(0);
    326: }

    330: PetscErrorCode MatTranspose_SeqBDiag(Mat A,Mat *matout)
    331: {
    332:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data,*anew;
    333:   Mat            tmat;
    335:   PetscInt       i,j,k,d,nd = a->nd,*diag = a->diag,*diagnew;
    336:   PetscInt       bs = A->bs,kshift,shifto,shiftn;
    337:   PetscScalar    *dwork,*dvnew;

    340:   PetscMalloc((nd+1)*sizeof(PetscInt),&diagnew);
    341:   for (i=0; i<nd; i++) {
    342:     diagnew[i] = -diag[nd-i-1]; /* assume sorted in descending order */
    343:   }
    344:   MatCreate(A->comm,&tmat);
    345:   MatSetSizes(tmat,A->n,A->m,A->n,A->m);
    346:   MatSetType(tmat,A->type_name);
    347:   MatSeqBDiagSetPreallocation(tmat,nd,bs,diagnew,PETSC_NULL);
    348:   PetscFree(diagnew);
    349:   anew = (Mat_SeqBDiag*)tmat->data;
    350:   for (d=0; d<nd; d++) {
    351:     dvnew = anew->diagv[d];
    352:     dwork = a->diagv[nd-d-1];
    353:     if (anew->bdlen[d] != a->bdlen[nd-d-1]) SETERRQ(PETSC_ERR_ARG_SIZ,"Incompatible diagonal lengths");
    354:     shifto = a->diag[nd-d-1];
    355:     shiftn = anew->diag[d];
    356:     if (shifto > 0)  shifto = bs*bs*shifto; else shifto = 0;
    357:     if (shiftn > 0)  shiftn = bs*bs*shiftn; else shiftn = 0;
    358:     if (bs == 1) {
    359:       for (k=0; k<anew->bdlen[d]; k++) dvnew[shiftn+k] = dwork[shifto+k];
    360:     } else {
    361:       for (k=0; k<anew->bdlen[d]; k++) {
    362:         kshift = k*bs*bs;
    363:         for (i=0; i<bs; i++) {        /* i = local row */
    364:           for (j=0; j<bs; j++) {        /* j = local column */
    365:             dvnew[shiftn + kshift + j + i*bs] = dwork[shifto + kshift + j*bs + i];
    366:           }
    367:         }
    368:       }
    369:     }
    370:   }
    371:   MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
    372:   MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
    373:   if (matout) {
    374:     *matout = tmat;
    375:   } else {
    376:     /* This isn't really an in-place transpose ... but free data
    377:        structures from a.  We should fix this. */
    378:     if (!a->user_alloc) { /* Free the actual diagonals */
    379:       for (i=0; i<a->nd; i++) {
    380:         if (a->diag[i] > 0) {
    381:           PetscScalar *dummy = a->diagv[i] + bs*bs*a->diag[i];
    382:           PetscFree(dummy);
    383:         } else {
    384:           PetscFree(a->diagv[i]);
    385:         }
    386:       }
    387:     }
    388:     if (a->pivot) {PetscFree(a->pivot);}
    389:     PetscFree(a->diagv);
    390:     PetscFree(a->diag);
    391:     PetscFree(a->colloc);
    392:     PetscFree(a->dvalue);
    393:     PetscFree(a);
    394:     PetscMemcpy(A,tmat,sizeof(struct _p_Mat));
    395:     PetscHeaderDestroy(tmat);
    396:   }
    397:   return(0);
    398: }

    400: /* ----------------------------------------------------------------*/

    405: PetscErrorCode MatView_SeqBDiag_Binary(Mat A,PetscViewer viewer)
    406: {
    407:   Mat_SeqBDiag   *a = (Mat_SeqBDiag*)A->data;
    409:   PetscInt       i,ict,*col_lens,*cval,*col,nz;
    410:   PetscScalar    *anonz,*val;
    411:   int            fd;

    414:   PetscViewerBinaryGetDescriptor(viewer,&fd);

    416:   /* For MATSEQBDIAG format,maxnz = nz */
    417:   PetscMalloc((4+A->m)*sizeof(PetscInt),&col_lens);
    418:   col_lens[0] = MAT_FILE_COOKIE;
    419:   col_lens[1] = A->m;
    420:   col_lens[2] = A->n;
    421:   col_lens[3] = a->maxnz;

    423:   /* Should do translation using less memory; this is just a quick initial version */
    424:   PetscMalloc((a->maxnz)*sizeof(PetscInt),&cval);
    425:   PetscMalloc((a->maxnz)*sizeof(PetscScalar),&anonz);

    427:   ict = 0;
    428:   for (i=0; i<A->m; i++) {
    429:     MatGetRow_SeqBDiag(A,i,&nz,&col,&val);
    430:     col_lens[4+i] = nz;
    431:     PetscMemcpy(&cval[ict],col,nz*sizeof(PetscInt));
    432:     PetscMemcpy(&anonz[ict],anonz,nz*sizeof(PetscScalar));
    433:     MatRestoreRow_SeqBDiag(A,i,&nz,&col,&val);
    434:     ict += nz;
    435:   }
    436:   if (ict != a->maxnz) SETERRQ(PETSC_ERR_PLIB,"Error in nonzero count");

    438:   /* Store lengths of each row and write (including header) to file */
    439:   PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,PETSC_TRUE);
    440:   PetscFree(col_lens);

    442:   /* Store column indices (zero start index) */
    443:   PetscBinaryWrite(fd,cval,a->maxnz,PETSC_INT,PETSC_FALSE);

    445:   /* Store nonzero values */
    446:   PetscBinaryWrite(fd,anonz,a->maxnz,PETSC_SCALAR,PETSC_FALSE);
    447:   return(0);
    448: }

    452: PetscErrorCode MatView_SeqBDiag_ASCII(Mat A,PetscViewer viewer)
    453: {
    454:   Mat_SeqBDiag      *a = (Mat_SeqBDiag*)A->data;
    455:   const char        *name;
    456:   PetscErrorCode    ierr;
    457:   PetscInt          *col,i,j,len,diag,nr = A->m,bs = A->bs,iprint,nz;
    458:   PetscScalar       *val,*dv,zero = 0.0;
    459:   PetscViewerFormat format;

    462:   PetscObjectGetName((PetscObject)A,&name);
    463:   PetscViewerGetFormat(viewer,&format);
    464:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    465:     PetscInt nline = PetscMin(10,a->nd),k,nk,np;
    466:     if (a->user_alloc) {
    467:       PetscViewerASCIIPrintf(viewer,"block size=%D, number of diagonals=%D, user-allocated storage\n",bs,a->nd);
    468:     } else {
    469:       PetscViewerASCIIPrintf(viewer,"block size=%D, number of diagonals=%D, PETSc-allocated storage\n",bs,a->nd);
    470:     }
    471:     nk = (a->nd-1)/nline + 1;
    472:     for (k=0; k<nk; k++) {
    473:       PetscViewerASCIIPrintf(viewer,"diag numbers:");
    474:       np = PetscMin(nline,a->nd - nline*k);
    475:       PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    476:       for (i=0; i<np; i++) {
    477:         PetscViewerASCIIPrintf(viewer,"  %D",a->diag[i+nline*k]);
    478:       }
    479:       PetscViewerASCIIPrintf(viewer,"\n");
    480:       PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    481:     }
    482:   } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
    483:     PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    484:     PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",nr,A->n);
    485:     PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);
    486:     PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz);
    487:     PetscViewerASCIIPrintf(viewer,"zzz = [\n");
    488:     for (i=0; i<A->m; i++) {
    489:       MatGetRow_SeqBDiag(A,i,&nz,&col,&val);
    490:       for (j=0; j<nz; j++) {
    491:         if (val[j] != zero) {
    492: #if defined(PETSC_USE_COMPLEX)
    493:           PetscViewerASCIIPrintf(viewer,"%D %D  %18.16e  %18.16ei \n",
    494:              i+1,col[j]+1,PetscRealPart(val[j]),PetscImaginaryPart(val[j]));
    495: #else
    496:           PetscViewerASCIIPrintf(viewer,"%D %D  %18.16ei \n",i+1,col[j]+1,val[j]);
    497: #endif
    498:         }
    499:       }
    500:       MatRestoreRow_SeqBDiag(A,i,&nz,&col,&val);
    501:     }
    502:     PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);
    503:     PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    504:   } else if (format == PETSC_VIEWER_ASCII_IMPL) {
    505:     PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    506:     if (bs == 1) { /* diagonal format */
    507:       for (i=0; i<a->nd; i++) {
    508:         dv   = a->diagv[i];
    509:         diag = a->diag[i];
    510:         PetscViewerASCIIPrintf(viewer,"\n<diagonal %D>\n",diag);
    511:         /* diag[i] is (row-col)/bs */
    512:         if (diag > 0) {  /* lower triangle */
    513:           len  = a->bdlen[i];
    514:           for (j=0; j<len; j++) {
    515:             if (dv[diag+j] != zero) {
    516: #if defined(PETSC_USE_COMPLEX)
    517:               if (PetscImaginaryPart(dv[diag+j]) != 0.0) {
    518:                 PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %e + %e i\n",
    519:                           j+diag,j,PetscRealPart(dv[diag+j]),PetscImaginaryPart(dv[diag+j]));
    520:               } else {
    521:                 PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %e\n",j+diag,j,PetscRealPart(dv[diag+j]));
    522:               }
    523: #else
    524:               PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %e\n",j+diag,j,dv[diag+j]);

    526: #endif
    527:             }
    528:           }
    529:         } else {         /* upper triangle, including main diagonal */
    530:           len  = a->bdlen[i];
    531:           for (j=0; j<len; j++) {
    532:             if (dv[j] != zero) {
    533: #if defined(PETSC_USE_COMPLEX)
    534:               if (PetscImaginaryPart(dv[j]) != 0.0) {
    535:                 PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %g + %g i\n",
    536:                                          j,j-diag,PetscRealPart(dv[j]),PetscImaginaryPart(dv[j]));
    537:               } else {
    538:                 PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %g\n",j,j-diag,PetscRealPart(dv[j]));
    539:               }
    540: #else
    541:               PetscViewerASCIIPrintf(viewer,"A[ %D , %D ] = %g\n",j,j-diag,dv[j]);
    542: #endif
    543:             }
    544:           }
    545:         }
    546:       }
    547:     } else {  /* Block diagonals */
    548:       PetscInt d,k,kshift;
    549:       for (d=0; d< a->nd; d++) {
    550:         dv   = a->diagv[d];
    551:         diag = a->diag[d];
    552:         len  = a->bdlen[d];
    553:         PetscViewerASCIIPrintf(viewer,"\n<diagonal %D>\n", diag);
    554:         if (diag > 0) {                /* lower triangle */
    555:           for (k=0; k<len; k++) {
    556:             kshift = (diag+k)*bs*bs;
    557:             for (i=0; i<bs; i++) {
    558:               iprint = 0;
    559:               for (j=0; j<bs; j++) {
    560:                 if (dv[kshift + j*bs + i] != zero) {
    561:                   iprint = 1;
    562: #if defined(PETSC_USE_COMPLEX)
    563:                   if (PetscImaginaryPart(dv[kshift + j*bs + i])){
    564:                     PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e + %5.2e i  ",(k+diag)*bs+i,k*bs+j,
    565:                       PetscRealPart(dv[kshift + j*bs + i]),PetscImaginaryPart(dv[kshift + j*bs + i]));
    566:                   } else {
    567:                     PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e   ",(k+diag)*bs+i,k*bs+j,
    568:                       PetscRealPart(dv[kshift + j*bs + i]));
    569:                   }
    570: #else
    571:                   PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e   ",(k+diag)*bs+i,k*bs+j,
    572:                       dv[kshift + j*bs + i]);
    573: #endif
    574:                 }
    575:               }
    576:               if (iprint) {PetscViewerASCIIPrintf(viewer,"\n");}
    577:             }
    578:           }
    579:         } else {                /* upper triangle, including main diagonal */
    580:           for (k=0; k<len; k++) {
    581:             kshift = k*bs*bs;
    582:             for (i=0; i<bs; i++) {
    583:               iprint = 0;
    584:               for (j=0; j<bs; j++) {
    585:                 if (dv[kshift + j*bs + i] != zero) {
    586:                   iprint = 1;
    587: #if defined(PETSC_USE_COMPLEX)
    588:                   if (PetscImaginaryPart(dv[kshift + j*bs + i])){
    589:                     PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e + %5.2e i  ",k*bs+i,(k-diag)*bs+j,
    590:                        PetscRealPart(dv[kshift + j*bs + i]),PetscImaginaryPart(dv[kshift + j*bs + i]));
    591:                   } else {
    592:                     PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e   ",k*bs+i,(k-diag)*bs+j,
    593:                        PetscRealPart(dv[kshift + j*bs + i]));
    594:                   }
    595: #else
    596:                   PetscViewerASCIIPrintf(viewer,"A[%D,%D]=%5.2e   ",k*bs+i,(k-diag)*bs+j,
    597:                      dv[kshift + j*bs + i]);
    598: #endif
    599:                 }
    600:               }
    601:               if (iprint) {PetscViewerASCIIPrintf(viewer,"\n");}
    602:             }
    603:           }
    604:         }
    605:       }
    606:     }
    607:     PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    608:   } else {
    609:     PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    610:     /* the usual row format (PETSC_VIEWER_ASCII_NONZERO_ONLY) */
    611:     for (i=0; i<A->m; i++) {
    612:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
    613:       MatGetRow_SeqBDiag(A,i,&nz,&col,&val);
    614:       for (j=0; j<nz; j++) {
    615: #if defined(PETSC_USE_COMPLEX)
    616:         if (PetscImaginaryPart(val[j]) != 0.0 && PetscRealPart(val[j]) != 0.0) {
    617:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",col[j],PetscRealPart(val[j]),PetscImaginaryPart(val[j]));
    618:         } else if (PetscRealPart(val[j]) != 0.0) {
    619:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",col[j],PetscRealPart(val[j]));
    620:         }
    621: #else
    622:         if (val[j] != 0.0) {PetscViewerASCIIPrintf(viewer," (%D, %g) ",col[j],val[j]);}
    623: #endif
    624:       }
    625:       PetscViewerASCIIPrintf(viewer,"\n");
    626:       MatRestoreRow_SeqBDiag(A,i,&nz,&col,&val);
    627:     }
    628:     PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    629:   }
    630:   PetscViewerFlush(viewer);
    631:   return(0);
    632: }

    636: static PetscErrorCode MatView_SeqBDiag_Draw(Mat A,PetscViewer viewer)
    637: {
    638:   PetscDraw      draw;
    639:   PetscReal      xl,yl,xr,yr,w,h;
    641:   PetscInt       nz,*col,i,j,nr = A->m;
    642:   PetscTruth     isnull;

    645:   PetscViewerDrawGetDraw(viewer,0,&draw);
    646:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

    648:   xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0;
    649:   xr += w; yr += h; xl = -w; yl = -h;
    650:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);

    652:   /* loop over matrix elements drawing boxes; we really should do this
    653:      by diagonals.  What do we really want to draw here: nonzeros,
    654:      allocated space? */
    655:   for (i=0; i<nr; i++) {
    656:     yl = nr - i - 1.0; yr = yl + 1.0;
    657:     MatGetRow_SeqBDiag(A,i,&nz,&col,0);
    658:     for (j=0; j<nz; j++) {
    659:       xl = col[j]; xr = xl + 1.0;
    660:       PetscDrawRectangle(draw,xl,yl,xr,yr,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,
    661:                            PETSC_DRAW_BLACK,PETSC_DRAW_BLACK);
    662:     }
    663:     MatRestoreRow_SeqBDiag(A,i,&nz,&col,0);
    664:   }
    665:   PetscDrawFlush(draw);
    666:   PetscDrawPause(draw);
    667:   return(0);
    668: }

    672: PetscErrorCode MatView_SeqBDiag(Mat A,PetscViewer viewer)
    673: {
    675:   PetscTruth     iascii,isbinary,isdraw;

    678:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    679:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
    680:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
    681:   if (iascii) {
    682:     MatView_SeqBDiag_ASCII(A,viewer);
    683:   } else if (isbinary) {
    684:     MatView_SeqBDiag_Binary(A,viewer);
    685:   } else if (isdraw) {
    686:     MatView_SeqBDiag_Draw(A,viewer);
    687:   } else {
    688:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by BDiag matrices",((PetscObject)viewer)->type_name);
    689:   }
    690:   return(0);
    691: }