Actual source code: matimpl.h

     5:  #include petscmat.h

      7: /*
      8:   This file defines the parts of the matrix data structure that are
      9:   shared by all matrix types.
     10: */

     12: /*
     13:     If you add entries here also add them to the MATOP enum
     14:     in include/petscmat.h and include/finclude/petscmat.h
     15: */
     16: typedef struct _MatOps *MatOps;
     17: struct _MatOps {
     18:   /* 0*/
     19:   PetscErrorCode (*setvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const MatScalar[],InsertMode);
     20:   PetscErrorCode (*getrow)(Mat,PetscInt,PetscInt *,PetscInt*[],PetscScalar*[]);
     21:   PetscErrorCode (*restorerow)(Mat,PetscInt,PetscInt *,PetscInt *[],PetscScalar *[]);
     22:   PetscErrorCode (*mult)(Mat,Vec,Vec);
     23:   PetscErrorCode (*multadd)(Mat,Vec,Vec,Vec);
     24:   /* 5*/
     25:   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
     26:   PetscErrorCode (*multtransposeadd)(Mat,Vec,Vec,Vec);
     27:   PetscErrorCode (*solve)(Mat,Vec,Vec);
     28:   PetscErrorCode (*solveadd)(Mat,Vec,Vec,Vec);
     29:   PetscErrorCode (*solvetranspose)(Mat,Vec,Vec);
     30:   /*10*/
     31:   PetscErrorCode (*solvetransposeadd)(Mat,Vec,Vec,Vec);
     32:   PetscErrorCode (*lufactor)(Mat,IS,IS,MatFactorInfo*);
     33:   PetscErrorCode (*choleskyfactor)(Mat,IS,MatFactorInfo*);
     34:   PetscErrorCode (*relax)(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
     35:   PetscErrorCode (*transpose)(Mat,Mat *);
     36:   /*15*/
     37:   PetscErrorCode (*getinfo)(Mat,MatInfoType,MatInfo*);
     38:   PetscErrorCode (*equal)(Mat,Mat,PetscTruth *);
     39:   PetscErrorCode (*getdiagonal)(Mat,Vec);
     40:   PetscErrorCode (*diagonalscale)(Mat,Vec,Vec);
     41:   PetscErrorCode (*norm)(Mat,NormType,PetscReal*);
     42:   /*20*/
     43:   PetscErrorCode (*assemblybegin)(Mat,MatAssemblyType);
     44:   PetscErrorCode (*assemblyend)(Mat,MatAssemblyType);
     45:   PetscErrorCode (*compress)(Mat);
     46:   PetscErrorCode (*setoption)(Mat,MatOption);
     47:   PetscErrorCode (*zeroentries)(Mat);
     48:   /*25*/
     49:   PetscErrorCode (*zerorows)(Mat,PetscInt,const PetscInt[],PetscScalar);
     50:   PetscErrorCode (*lufactorsymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
     51:   PetscErrorCode (*lufactornumeric)(Mat,MatFactorInfo*,Mat*);
     52:   PetscErrorCode (*choleskyfactorsymbolic)(Mat,IS,MatFactorInfo*,Mat*);
     53:   PetscErrorCode (*choleskyfactornumeric)(Mat,MatFactorInfo*,Mat*);
     54:   /*30*/
     55:   PetscErrorCode (*setuppreallocation)(Mat);
     56:   PetscErrorCode (*ilufactorsymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
     57:   PetscErrorCode (*iccfactorsymbolic)(Mat,IS,MatFactorInfo*,Mat*);
     58:   PetscErrorCode (*getarray)(Mat,PetscScalar**);
     59:   PetscErrorCode (*restorearray)(Mat,PetscScalar**);
     60:   /*35*/
     61:   PetscErrorCode (*duplicate)(Mat,MatDuplicateOption,Mat*);
     62:   PetscErrorCode (*forwardsolve)(Mat,Vec,Vec);
     63:   PetscErrorCode (*backwardsolve)(Mat,Vec,Vec);
     64:   PetscErrorCode (*ilufactor)(Mat,IS,IS,MatFactorInfo*);
     65:   PetscErrorCode (*iccfactor)(Mat,IS,MatFactorInfo*);
     66:   /*40*/
     67:   PetscErrorCode (*axpy)(Mat,PetscScalar,Mat,MatStructure);
     68:   PetscErrorCode (*getsubmatrices)(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *[]);
     69:   PetscErrorCode (*increaseoverlap)(Mat,PetscInt,IS[],PetscInt);
     70:   PetscErrorCode (*getvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscScalar []);
     71:   PetscErrorCode (*copy)(Mat,Mat,MatStructure);
     72:   /*45*/
     73:   PetscErrorCode (*printhelp)(Mat);
     74:   PetscErrorCode (*scale)(Mat,PetscScalar);
     75:   PetscErrorCode (*shift)(Mat,PetscScalar);
     76:   PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
     77:   PetscErrorCode (*iludtfactor)(Mat,IS,IS,MatFactorInfo*,Mat *);
     78:   /*50*/
     79:   PetscErrorCode (*setblocksize)(Mat,PetscInt);
     80:   PetscErrorCode (*getrowij)(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
     81:   PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscTruth,PetscInt *,PetscInt *[],PetscInt *[],PetscTruth *);
     82:   PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
     83:   PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
     84:   /*55*/
     85:   PetscErrorCode (*fdcoloringcreate)(Mat,ISColoring,MatFDColoring);
     86:   PetscErrorCode (*coloringpatch)(Mat,PetscInt,PetscInt,ISColoringValue[],ISColoring*);
     87:   PetscErrorCode (*setunfactored)(Mat);
     88:   PetscErrorCode (*permute)(Mat,IS,IS,Mat*);
     89:   PetscErrorCode (*setvaluesblocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
     90:   /*60*/
     91:   PetscErrorCode (*getsubmatrix)(Mat,IS,IS,PetscInt,MatReuse,Mat*);
     92:   PetscErrorCode (*destroy)(Mat);
     93:   PetscErrorCode (*view)(Mat,PetscViewer);
     94:   PetscErrorCode (*getmaps)(Mat,PetscMap*,PetscMap*);
     95:   PetscErrorCode (*usescaledform)(Mat,PetscTruth);
     96:   /*65*/
     97:   PetscErrorCode (*scalesystem)(Mat,Vec,Vec);
     98:   PetscErrorCode (*unscalesystem)(Mat,Vec,Vec);
     99:   PetscErrorCode (*setlocaltoglobalmapping)(Mat,ISLocalToGlobalMapping);
    100:   PetscErrorCode (*setvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
    101:   PetscErrorCode (*zerorowslocal)(Mat,PetscInt,const PetscInt[],PetscScalar);
    102:   /*70*/
    103:   PetscErrorCode (*getrowmax)(Mat,Vec);
    104:   PetscErrorCode (*convert)(Mat, MatType,MatReuse,Mat*);
    105:   PetscErrorCode (*setcoloring)(Mat,ISColoring);
    106:   PetscErrorCode (*setvaluesadic)(Mat,void*);
    107:   PetscErrorCode (*setvaluesadifor)(Mat,PetscInt,void*);
    108:   /*75*/
    109:   PetscErrorCode (*fdcoloringapply)(Mat,MatFDColoring,Vec,MatStructure*,void*);
    110:   PetscErrorCode (*setfromoptions)(Mat);
    111:   PetscErrorCode (*multconstrained)(Mat,Vec,Vec);
    112:   PetscErrorCode (*multtransposeconstrained)(Mat,Vec,Vec);
    113:   PetscErrorCode (*ilufactorsymbolicconstrained)(Mat,IS,IS,double,PetscInt,PetscInt,Mat *);
    114:   /*80*/
    115:   PetscErrorCode (*permutesparsify)(Mat, PetscInt, double, double, IS, IS, Mat *);
    116:   PetscErrorCode (*mults)(Mat, Vecs, Vecs);
    117:   PetscErrorCode (*solves)(Mat, Vecs, Vecs);
    118:   PetscErrorCode (*getinertia)(Mat,PetscInt*,PetscInt*,PetscInt*);
    119:   PetscErrorCode (*load)(PetscViewer, MatType,Mat*);
    120:   /*85*/
    121:   PetscErrorCode (*issymmetric)(Mat,PetscReal,PetscTruth*);
    122:   PetscErrorCode (*ishermitian)(Mat,PetscTruth*);
    123:   PetscErrorCode (*isstructurallysymmetric)(Mat,PetscTruth*);
    124:   PetscErrorCode (*pbrelax)(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
    125:   PetscErrorCode (*getvecs)(Mat,Vec*,Vec*);
    126:   /*90*/
    127:   PetscErrorCode (*matmult)(Mat,Mat,MatReuse,PetscReal,Mat*);
    128:   PetscErrorCode (*matmultsymbolic)(Mat,Mat,PetscReal,Mat*);
    129:   PetscErrorCode (*matmultnumeric)(Mat,Mat,Mat);
    130:   PetscErrorCode (*ptap)(Mat,Mat,MatReuse,PetscReal,Mat*);
    131:   PetscErrorCode (*ptapsymbolic)(Mat,Mat,PetscReal,Mat*); /* double dispatch wrapper routine */
    132:   /*95*/
    133:   PetscErrorCode (*ptapnumeric)(Mat,Mat,Mat);             /* double dispatch wrapper routine */
    134:   PetscErrorCode (*matmulttranspose)(Mat,Mat,MatReuse,PetscReal,Mat*);
    135:   PetscErrorCode (*matmulttransposesymbolic)(Mat,Mat,PetscReal,Mat*);
    136:   PetscErrorCode (*matmulttransposenumeric)(Mat,Mat,Mat);
    137:   PetscErrorCode (*ptapsymbolic_seqaij)(Mat,Mat,PetscReal,Mat*); /* actual implememtation, A=seqaij */
    138:   /*100*/
    139:   PetscErrorCode (*ptapnumeric_seqaij)(Mat,Mat,Mat);             /* actual implememtation, A=seqaij */
    140:   PetscErrorCode (*ptapsymbolic_mpiaij)(Mat,Mat,PetscReal,Mat*); /* actual implememtation, A=mpiaij */
    141:   PetscErrorCode (*ptapnumeric_mpiaij)(Mat,Mat,Mat);             /* actual implememtation, A=mpiaij */
    142:   PetscErrorCode (*conjugate)(Mat);                              /* complex conjugate */
    143: };
    144: /*
    145:     If you add MatOps entries above also add them to the MATOP enum
    146:     in include/petscmat.h and include/finclude/petscmat.h
    147: */

    149: /*
    150:    Utility private matrix routines
    151: */
    152: EXTERN PetscErrorCode MatConvert_Basic(Mat, MatType,MatReuse,Mat*);
    153: EXTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
    154: EXTERN PetscErrorCode MatView_Private(Mat);
    155: EXTERN PetscErrorCode MatGetPetscMaps_Petsc(Mat,PetscMap *,PetscMap *);
    156: EXTERN PetscErrorCode MatHeaderCopy(Mat,Mat);
    157: EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
    158: EXTERN PetscErrorCode MatAXPYGetxtoy_Private(PetscInt,PetscInt*,PetscInt*,PetscInt*, PetscInt*,PetscInt*,PetscInt*, PetscInt**);
    159: EXTERN PetscErrorCode MatPtAP_Basic(Mat,Mat,MatReuse,PetscReal,Mat*);

    161: /*
    162:   The stash is used to temporarily store inserted matrix values that
    163:   belong to another processor. During the assembly phase the stashed
    164:   values are moved to the correct processor and
    165: */

    167: typedef struct {
    168:   PetscInt      nmax;                   /* maximum stash size */
    169:   PetscInt      umax;                   /* user specified max-size */
    170:   PetscInt      oldnmax;                /* the nmax value used previously */
    171:   PetscInt      n;                      /* stash size */
    172:   PetscInt      bs;                     /* block size of the stash */
    173:   PetscInt      reallocs;               /* preserve the no of mallocs invoked */
    174:   PetscInt      *idx;                   /* global row numbers in stash */
    175:   PetscInt      *idy;                   /* global column numbers in stash */
    176:   MatScalar     *array;                 /* array to hold stashed values */
    177:   /* The following variables are used for communication */
    178:   MPI_Comm      comm;
    179:   PetscMPIInt   size,rank;
    180:   PetscMPIInt   tag1,tag2;
    181:   MPI_Request   *send_waits;            /* array of send requests */
    182:   MPI_Request   *recv_waits;            /* array of receive requests */
    183:   MPI_Status    *send_status;           /* array of send status */
    184:   PetscInt      nsends,nrecvs;          /* numbers of sends and receives */
    185:   MatScalar     *svalues;               /* sending data */
    186:   MatScalar     **rvalues;              /* receiving data (values) */
    187:   PetscInt      **rindices;             /* receiving data (indices) */
    188:   PetscMPIInt   *nprocs;                /* tmp data used both during scatterbegin and end */
    189:   PetscInt      nprocessed;             /* number of messages already processed */
    190: } MatStash;

    192: EXTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
    193: EXTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
    194: EXTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
    195: EXTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
    196: EXTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
    197: EXTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[]);
    198: EXTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt);
    199: EXTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
    200: EXTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
    201: EXTERN PetscErrorCode MatStashScatterBegin_Private(MatStash*,PetscInt*);
    202: EXTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,MatScalar**,PetscInt*);

    204: #define FACTOR_LU       1
    205: #define FACTOR_CHOLESKY 2

    207: typedef struct {
    208:   PetscInt   dim;
    209:   PetscInt   dims[4];
    210:   PetscInt   starts[4];
    211:   PetscTruth noc;        /* this is a single component problem, hence user will not set MatStencil.c */
    212: } MatStencilInfo;

    214: /* Info about using compressed row format */
    215: typedef struct {
    216:   PetscTruth use;
    217:   PetscInt   nrows;                         /* number of non-zero rows */
    218:   PetscInt   *i;                            /* compressed row pointer  */
    219:   PetscInt   *rindex;                       /* compressed row index               */
    220:   PetscTruth checked;                       /* if compressed row format have been checked for */
    221: } Mat_CompressedRow;
    222: EXTERN PetscErrorCode Mat_CheckCompressedRow(Mat,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);

    224: struct _p_Mat {
    225:   PETSCHEADER(struct _MatOps);
    226:   PetscMap               rmap,cmap;
    227:   void                   *data;            /* implementation-specific data */
    228:   PetscInt               factor;           /* 0, FACTOR_LU, or FACTOR_CHOLESKY */
    229:   PetscTruth             assembled;        /* is the matrix assembled? */
    230:   PetscTruth             was_assembled;    /* new values inserted into assembled mat */
    231:   PetscInt               num_ass;          /* number of times matrix has been assembled */
    232:   PetscTruth             same_nonzero;     /* matrix has same nonzero pattern as previous */
    233:   PetscInt               M,N;             /* global numbers of rows, columns */
    234:   PetscInt               m,n;             /* local numbers of rows, columns */
    235:   MatInfo                info;             /* matrix information */
    236:   ISLocalToGlobalMapping mapping;          /* mapping used in MatSetValuesLocal() */
    237:   ISLocalToGlobalMapping bmapping;         /* mapping used in MatSetValuesBlockedLocal() */
    238:   InsertMode             insertmode;       /* have values been inserted in matrix or added? */
    239:   MatStash               stash,bstash;     /* used for assembling off-proc mat emements */
    240:   MatNullSpace           nullsp;
    241:   PetscTruth             preallocated;
    242:   MatStencilInfo         stencil;          /* information for structured grid */
    243:   PetscTruth             symmetric,hermitian,structurally_symmetric;
    244:   PetscTruth             symmetric_set,hermitian_set,structurally_symmetric_set; /* if true, then corresponding flag is correct*/
    245:   PetscTruth             symmetric_eternal;
    246:   PetscInt               bs;
    247:   void                   *spptr;          /* pointer for special library like SuperLU */
    248: };

    250: #define MatPreallocated(A)  ((!(A)->preallocated) ? MatSetUpPreallocation(A) : 0)
    252: /*
    253:     Frees the a, i, and j arrays from the XAIJ (AIJ, BAIJ, and SBAIJ) matrix types
    254: */
    257: PETSC_STATIC_INLINE PetscErrorCode MatSeqXAIJFreeAIJ(PetscTruth singlemalloc,PetscScalar **a,PetscInt **j,PetscInt **i) {
    259:                                      if (singlemalloc) {
    260:                                        PetscFree3(*a,*j,*i);
    261:                                      } else {
    262:                                        if (*a) {PetscFree(*a);}
    263:                                        if (*j) {PetscFree(*j);}
    264:                                        if (*i) {PetscFree(*i);}
    265:                                      }
    266:                                      *a = 0; *j = 0; *i = 0;
    267:                                      return 0;
    268:                                    }

    270: /*
    271:     Allocates larger a, i, and j arrays for the XAIJ (AIJ, BAIJ, and SBAIJ) matrix types
    272: */
    273: #define MatSeqXAIJReallocateAIJ(A,BS2,NROW,ROW,COL,RMAX,AA,AI,AJ,AM,RP,AP,AIMAX,NONEW) \
    274:       if (NROW >= RMAX) { \
    275:         /* there is no extra room in row, therefore enlarge */ \
    276:         PetscInt    new_nz = AI[AM] + CHUNKSIZE,len,*new_i,*new_j; \
    277:         PetscScalar *new_a; \
    278:  \
    279:         if (NONEW == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc",ROW,COL); \
    280:         /* malloc new storage space */ \
    281:         PetscMalloc3(BS2*new_nz,PetscScalar,&new_a,new_nz,PetscInt,&new_j,AM+1,PetscInt,&new_i);\
    282:  \
    283:         /* copy over old data into new slots */ \
    284:         for (ii=0; ii<ROW+1; ii++) {new_i[ii] = AI[ii];} \
    285:         for (ii=ROW+1; ii<AM+1; ii++) {new_i[ii] = AI[ii]+CHUNKSIZE;} \
    286:         PetscMemcpy(new_j,AJ,(AI[ROW]+NROW)*sizeof(PetscInt)); \
    287:         len = (new_nz - CHUNKSIZE - AI[ROW] - NROW); \
    288:         PetscMemcpy(new_j+AI[ROW]+NROW+CHUNKSIZE,AJ+AI[ROW]+NROW,len*sizeof(PetscInt)); \
    289:         PetscMemcpy(new_a,AA,BS2*(AI[ROW]+NROW)*sizeof(PetscScalar)); \
    290:         PetscMemzero(new_a+BS2*(AI[ROW]+NROW),BS2*CHUNKSIZE*sizeof(MatScalar));\
    291:         PetscMemcpy(new_a+BS2*(AI[ROW]+NROW+CHUNKSIZE),AA+BS2*(AI[ROW]+NROW),BS2*len*sizeof(PetscScalar));  \
    292:         /* free up old matrix storage */ \
    293:         MatSeqXAIJFreeAIJ(A->singlemalloc,&A->a,&A->j,&A->i);\
    294:         AA = A->a = new_a; AI = A->i = new_i; AJ = A->j = new_j;  \
    295:         A->singlemalloc = PETSC_TRUE; \
    296:  \
    297:         RP        = AJ + AI[ROW]; AP = AA + BS2*AI[ROW]; \
    298:         RMAX      = AIMAX[ROW] = AIMAX[ROW] + CHUNKSIZE; \
    299:         A->maxnz += CHUNKSIZE; \
    300:         A->reallocs++; \
    301:       } \

    303: /*
    304:     Object for partitioning graphs
    305: */

    307: typedef struct _MatPartitioningOps *MatPartitioningOps;
    308: struct _MatPartitioningOps {
    309:   PetscErrorCode (*apply)(MatPartitioning,IS*);
    310:   PetscErrorCode (*setfromoptions)(MatPartitioning);
    311:   PetscErrorCode (*destroy)(MatPartitioning);
    312:   PetscErrorCode (*view)(MatPartitioning,PetscViewer);
    313: };

    315: struct _p_MatPartitioning {
    316:   PETSCHEADER(struct _MatPartitioningOps);
    317:   Mat         adj;
    318:   PetscInt    *vertex_weights;
    319:   PetscReal   *part_weights;
    320:   PetscInt    n;                                 /* number of partitions */
    321:   void        *data;
    322:   PetscInt    setupcalled;
    323: };

    325: /*
    326:     MatFDColoring is used to compute Jacobian matrices efficiently
    327:   via coloring. The data structure is explained below in an example.

    329:    Color =   0    1     0    2   |   2      3       0
    330:    ---------------------------------------------------
    331:             00   01              |          05
    332:             10   11              |   14     15               Processor  0
    333:                        22    23  |          25
    334:                        32    33  |
    335:    ===================================================
    336:                                  |   44     45     46
    337:             50                   |          55               Processor 1
    338:                                  |   64            66
    339:    ---------------------------------------------------

    341:     ncolors = 4;

    343:     ncolumns      = {2,1,1,0}
    344:     columns       = {{0,2},{1},{3},{}}
    345:     nrows         = {4,2,3,3}
    346:     rows          = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
    347:     columnsforrow = {{0,0,2,2},{1,1},{4,3,3},{5,5,5}}
    348:     vscaleforrow  = {{,,,},{,},{,,},{,,}}
    349:     vwscale       = {dx(0),dx(1),dx(2),dx(3)}               MPI Vec
    350:     vscale        = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)}   Seq Vec

    352:     ncolumns      = {1,0,1,1}
    353:     columns       = {{6},{},{4},{5}}
    354:     nrows         = {3,0,2,2}
    355:     rows          = {{0,1,2},{},{1,2},{1,2}}
    356:     columnsforrow = {{6,0,6},{},{4,4},{5,5}}
    357:     vscaleforrow =  {{,,},{},{,},{,}}
    358:     vwscale       = {dx(4),dx(5),dx(6)}              MPI Vec
    359:     vscale        = {dx(0),dx(4),dx(5),dx(6)}        Seq Vec

    361:     See the routine MatFDColoringApply() for how this data is used
    362:     to compute the Jacobian.

    364: */

    366: struct  _p_MatFDColoring{
    367:   PETSCHEADER(int);
    368:   PetscInt   M,N,m;            /* total rows, columns; local rows */
    369:   PetscInt   rstart;           /* first row owned by local processor */
    370:   PetscInt   ncolors;          /* number of colors */
    371:   PetscInt   *ncolumns;        /* number of local columns for a color */
    372:   PetscInt   **columns;        /* lists the local columns of each color (using global column numbering) */
    373:   PetscInt   *nrows;           /* number of local rows for each color */
    374:   PetscInt   **rows;           /* lists the local rows for each color (using the local row numbering) */
    375:   PetscInt   **columnsforrow;  /* lists the corresponding columns for those rows (using the global column) */
    376:   PetscReal  error_rel;        /* square root of relative error in computing function */
    377:   PetscReal  umin;             /* minimum allowable u'dx value */
    378:   PetscInt   freq;             /* frequency at which new Jacobian is computed */
    379:   Vec        w1,w2,w3;         /* work vectors used in computing Jacobian */
    380:   PetscErrorCode (*f)(void);       /* function that defines Jacobian */
    381:   void       *fctx;            /* optional user-defined context for use by the function f */
    382:   PetscInt   **vscaleforrow;   /* location in vscale for each columnsforrow[] entry */
    383:   Vec        vscale;           /* holds FD scaling, i.e. 1/dx for each perturbed column */
    384:   PetscTruth usersetsrecompute;/* user determines when Jacobian is recomputed, via MatFDColoringSetRecompute() */
    385:   PetscTruth recompute;        /* used with usersetrecompute to determine if Jacobian should be recomputed */
    386:   Vec        F;                /* current value of user provided function; can set with MatFDColoringSetF() */
    387:   PetscInt   currentcolor;     /* color for which function evaluation is being done now */
    388: };

    390: /*
    391:    Null space context for preconditioner/operators
    392: */
    393: struct _p_MatNullSpace {
    394:   PETSCHEADER(int);
    395:   PetscTruth     has_cnst;
    396:   PetscInt       n;
    397:   Vec*           vecs;
    398:   Vec            vec;                   /* for out of place removals */
    399:   PetscErrorCode (*remove)(Vec,void*);  /* for user provided removal function */
    400:   void*          rmctx;                 /* context for remove() function */
    401: };

    403: /*
    404:    Checking zero pivot for LU, ILU preconditioners.
    405: */
    406: typedef struct {
    407:   PetscInt       nshift,nshift_max;
    408:   PetscReal      shift_amount,shift_lo,shift_hi,shift_top;
    409:   PetscTruth     lushift;
    410:   PetscReal      rs;  /* active row sum of abs(offdiagonals) */
    411:   PetscScalar    pv;  /* pivot of the active row */
    412: } LUShift_Ctx;

    416: /*@C
    417:    MatLUCheckShift_inline - shift the diagonals when zero pivot is detected on LU factor

    419:    Collective on Mat

    421:    Input Parameters:
    422: +  info - information about the matrix factorization
    423: .  sctx - pointer to the struct LUShift_Ctx
    424: -  newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot

    426:    Level: developer
    427: @*/
    428: #define MatLUCheckShift_inline(info,sctx,newshift) 0;\
    429: {\
    430:   PetscInt _newshift;\
    431:   PetscReal _zero = info->zeropivot*rs;\
    432:   if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
    433:     /* force |diag| > zeropivot*rs */\
    434:     if (!sctx.nshift){\
    435:       sctx.shift_amount = info->shiftnz;\
    436:     } else {\
    437:       sctx.shift_amount *= 2.0;\
    438:     }\
    439:     sctx.lushift = PETSC_TRUE;\
    440:     (sctx.nshift)++;\
    441:     _newshift = 1;\
    442:   } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
    443:     /* force matfactor to be diagonally dominant */\
    444:     if (sctx.nshift > sctx.nshift_max) {\
    445:       SETERRQ(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner");\
    446:     } else if (sctx.nshift == sctx.nshift_max) {\
    447:       info->shift_fraction = sctx.shift_hi;\
    448:       sctx.lushift        = PETSC_TRUE;\
    449:     } else {\
    450:       sctx.shift_lo = info->shift_fraction;\
    451:       info->shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;\
    452:       sctx.lushift  = PETSC_TRUE;\
    453:     }\
    454:     sctx.shift_amount = info->shift_fraction * sctx.shift_top;\
    455:     sctx.nshift++;\
    456:     _newshift = 1;\
    457:   } else if (PetscAbsScalar(sctx.pv) <= _zero){\
    458:     _newshift = -1;\
    459:   } else {\
    460:     _newshift = 0;\
    461:   }\
    462:   newshift = _newshift;\
    463: }

    465: /*
    466:    Checking zero pivot for Cholesky, ICC preconditioners.
    467: */
    468: typedef struct {
    469:   PetscInt       nshift;
    470:   PetscReal      shift_amount;
    471:   PetscTruth     chshift;
    472:   PetscReal      rs;  /* active row sum of abs(offdiagonals) */
    473:   PetscScalar    pv;  /* pivot of the active row */
    474: } ChShift_Ctx;

    478: /*@C
    479:    MatCholeskyCheckShift_inline -  shift the diagonals when zero pivot is detected on Cholesky factor

    481:    Collective on Mat

    483:    Input Parameters:
    484: +  info - information about the matrix factorization
    485: .  sctx - pointer to the struct CholeskyShift_Ctx
    486: -  newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot

    488:    Level: developer
    489:    Note: Unlike in the ILU case there is no exit condition on nshift:
    490:        we increase the shift until it converges. There is no guarantee that
    491:        this algorithm converges faster or slower, or is better or worse
    492:        than the ILU algorithm.
    493: @*/
    494: #define MatCholeskyCheckShift_inline(info,sctx,newshift) 0;\
    495: {\
    496:   PetscInt _newshift;\
    497:   PetscReal _zero = info->zeropivot*rs;\
    498:   if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
    499:     /* force |diag| > zeropivot*sctx.rs */\
    500:     if (!sctx.nshift){\
    501:       sctx.shift_amount = info->shiftnz;\
    502:     } else {\
    503:       sctx.shift_amount *= 2.0;\
    504:     }\
    505:     sctx.chshift = PETSC_TRUE;\
    506:     sctx.nshift++;\
    507:     _newshift = 1;\
    508:   } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
    509:     /* calculate a shift that would make this row diagonally dominant */\
    510:     sctx.shift_amount = PetscMax(sctx.rs+PetscAbs(PetscRealPart(sctx.pv)),1.1*sctx.shift_amount);\
    511:     sctx.chshift      = PETSC_TRUE;\
    512:     sctx.nshift++;\
    513:     _newshift = 1;\
    514:   } else if (PetscAbsScalar(sctx.pv) <= _zero){\
    515:     _newshift = -1;\
    516:   } else {\
    517:     _newshift = 0; \
    518:   }\
    519:   newshift = _newshift;\
    520: }

    522: #endif