Actual source code: matimpl.h

      2: #ifndef __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 (*getrowmax)(Mat,Vec,PetscInt[]);
     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,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
     81:   PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt *,PetscInt *[],PetscInt *[],PetscTruth *);
     82:   PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
     83:   PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscTruth,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 (*convertfrom)(Mat, MatType,MatReuse,Mat*);
     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 (*getrowmaxabs)(Mat,Vec,PetscInt[]);
    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:   PetscErrorCode (*setsizes)(Mat,PetscInt,PetscInt,PetscInt,PetscInt);
    144:   /*105*/
    145:   PetscErrorCode (*setvaluesrow)(Mat,PetscInt,const MatScalar[]);
    146:   PetscErrorCode (*realpart)(Mat);
    147:   PetscErrorCode (*imaginarypart)(Mat);
    148:   PetscErrorCode (*getrowuppertriangular)(Mat);
    149:   PetscErrorCode (*restorerowuppertriangular)(Mat);
    150:   /*110*/
    151:   PetscErrorCode (*matsolve)(Mat,Mat,Mat);
    152:   PetscErrorCode (*getredundantmatrix)(Mat,PetscInt,MPI_Comm,PetscInt,MatReuse,Mat*);
    153:   PetscErrorCode (*getrowmin)(Mat,Vec,PetscInt[]);
    154:   PetscErrorCode (*getcolumnvector)(Mat,Vec,PetscInt);
    155: };
    156: /*
    157:     If you add MatOps entries above also add them to the MATOP enum
    158:     in include/petscmat.h and include/finclude/petscmat.h
    159: */

    161: /*
    162:    Utility private matrix routines
    163: */
    164: EXTERN PetscErrorCode MatConvert_Basic(Mat, MatType,MatReuse,Mat*);
    165: EXTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
    166: EXTERN PetscErrorCode MatView_Private(Mat);

    168: EXTERN PetscErrorCode MatHeaderCopy(Mat,Mat);
    169: EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
    170: EXTERN PetscErrorCode MatAXPYGetxtoy_Private(PetscInt,PetscInt*,PetscInt*,PetscInt*, PetscInt*,PetscInt*,PetscInt*, PetscInt**);
    171: EXTERN PetscErrorCode MatPtAP_Basic(Mat,Mat,MatReuse,PetscReal,Mat*);
    172: EXTERN PetscErrorCode MatDiagonalSet_Default(Mat,Vec,InsertMode);

    174: /*
    175:   The stash is used to temporarily store inserted matrix values that
    176:   belong to another processor. During the assembly phase the stashed
    177:   values are moved to the correct processor and
    178: */

    180: typedef struct _MatStashSpace *PetscMatStashSpace;

    182: struct _MatStashSpace {
    183:   PetscMatStashSpace next;
    184:   MatScalar          *space_head,*val;
    185:   PetscInt           *idx,*idy;
    186:   PetscInt           total_space_size;
    187:   PetscInt           local_used;
    188:   PetscInt           local_remaining;
    189: };

    191: EXTERN PetscErrorCode PetscMatStashSpaceGet(PetscInt,PetscInt,PetscMatStashSpace *);
    192: EXTERN PetscErrorCode PetscMatStashSpaceContiguous(PetscInt,PetscMatStashSpace *,PetscScalar *,PetscInt *,PetscInt *);
    193: EXTERN PetscErrorCode PetscMatStashSpaceDestroy(PetscMatStashSpace);

    195: typedef struct {
    196:   PetscInt      nmax;                   /* maximum stash size */
    197:   PetscInt      umax;                   /* user specified max-size */
    198:   PetscInt      oldnmax;                /* the nmax value used previously */
    199:   PetscInt      n;                      /* stash size */
    200:   PetscInt      bs;                     /* block size of the stash */
    201:   PetscInt      reallocs;               /* preserve the no of mallocs invoked */
    202:   PetscMatStashSpace space_head,space;  /* linked list to hold stashed global row/column numbers and matrix values */
    203:   /* The following variables are used for communication */
    204:   MPI_Comm      comm;
    205:   PetscMPIInt   size,rank;
    206:   PetscMPIInt   tag1,tag2;
    207:   MPI_Request   *send_waits;            /* array of send requests */
    208:   MPI_Request   *recv_waits;            /* array of receive requests */
    209:   MPI_Status    *send_status;           /* array of send status */
    210:   PetscInt      nsends,nrecvs;          /* numbers of sends and receives */
    211:   MatScalar     *svalues;               /* sending data */
    212:   MatScalar     **rvalues;              /* receiving data (values) */
    213:   PetscInt      **rindices;             /* receiving data (indices) */
    214:   PetscMPIInt   *nprocs;                /* tmp data used both during scatterbegin and end */
    215:   PetscInt      nprocessed;             /* number of messages already processed */
    216: } MatStash;

    218: EXTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
    219: EXTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
    220: EXTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
    221: EXTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
    222: EXTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
    223: EXTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[]);
    224: EXTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt);
    225: EXTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
    226: EXTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
    227: EXTERN PetscErrorCode MatStashScatterBegin_Private(MatStash*,PetscInt*);
    228: EXTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,MatScalar**,PetscInt*);

    230: #define FACTOR_LU       1
    231: #define FACTOR_CHOLESKY 2

    233: typedef struct {
    234:   PetscInt   dim;
    235:   PetscInt   dims[4];
    236:   PetscInt   starts[4];
    237:   PetscTruth noc;        /* this is a single component problem, hence user will not set MatStencil.c */
    238: } MatStencilInfo;

    240: /* Info about using compressed row format */
    241: typedef struct {
    242:   PetscTruth use;
    243:   PetscInt   nrows;                         /* number of non-zero rows */
    244:   PetscInt   *i;                            /* compressed row pointer  */
    245:   PetscInt   *rindex;                       /* compressed row index               */
    246:   PetscTruth checked;                       /* if compressed row format have been checked for */
    247: } Mat_CompressedRow;
    248: EXTERN PetscErrorCode Mat_CheckCompressedRow(Mat,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);

    250: struct _p_Mat {
    251:   PETSCHEADER(struct _MatOps);
    252:   PetscMap               rmap,cmap;
    253:   void                   *data;            /* implementation-specific data */
    254:   PetscInt               factor;           /* 0, FACTOR_LU, or FACTOR_CHOLESKY */
    255:   PetscTruth             assembled;        /* is the matrix assembled? */
    256:   PetscTruth             was_assembled;    /* new values inserted into assembled mat */
    257:   PetscInt               num_ass;          /* number of times matrix has been assembled */
    258:   PetscTruth             same_nonzero;     /* matrix has same nonzero pattern as previous */
    259:   MatInfo                info;             /* matrix information */
    260:   ISLocalToGlobalMapping mapping;          /* mapping used in MatSetValuesLocal() */
    261:   ISLocalToGlobalMapping bmapping;         /* mapping used in MatSetValuesBlockedLocal() */
    262:   InsertMode             insertmode;       /* have values been inserted in matrix or added? */
    263:   MatStash               stash,bstash;     /* used for assembling off-proc mat emements */
    264:   MatNullSpace           nullsp;
    265:   PetscTruth             preallocated;
    266:   MatStencilInfo         stencil;          /* information for structured grid */
    267:   PetscTruth             symmetric,hermitian,structurally_symmetric;
    268:   PetscTruth             symmetric_set,hermitian_set,structurally_symmetric_set; /* if true, then corresponding flag is correct*/
    269:   PetscTruth             symmetric_eternal;
    270:   void                   *spptr;          /* pointer for special library like SuperLU */
    271: };

    273: #define MatPreallocated(A)  ((!(A)->preallocated) ? MatSetUpPreallocation(A) : 0)

    276: /*
    277:     Object for partitioning graphs
    278: */

    280: typedef struct _MatPartitioningOps *MatPartitioningOps;
    281: struct _MatPartitioningOps {
    282:   PetscErrorCode (*apply)(MatPartitioning,IS*);
    283:   PetscErrorCode (*setfromoptions)(MatPartitioning);
    284:   PetscErrorCode (*destroy)(MatPartitioning);
    285:   PetscErrorCode (*view)(MatPartitioning,PetscViewer);
    286: };

    288: struct _p_MatPartitioning {
    289:   PETSCHEADER(struct _MatPartitioningOps);
    290:   Mat         adj;
    291:   PetscInt    *vertex_weights;
    292:   PetscReal   *part_weights;
    293:   PetscInt    n;                                 /* number of partitions */
    294:   void        *data;
    295:   PetscInt    setupcalled;
    296: };

    298: /*
    299:     MatFDColoring is used to compute Jacobian matrices efficiently
    300:   via coloring. The data structure is explained below in an example.

    302:    Color =   0    1     0    2   |   2      3       0
    303:    ---------------------------------------------------
    304:             00   01              |          05
    305:             10   11              |   14     15               Processor  0
    306:                        22    23  |          25
    307:                        32    33  |
    308:    ===================================================
    309:                                  |   44     45     46
    310:             50                   |          55               Processor 1
    311:                                  |   64            66
    312:    ---------------------------------------------------

    314:     ncolors = 4;

    316:     ncolumns      = {2,1,1,0}
    317:     columns       = {{0,2},{1},{3},{}}
    318:     nrows         = {4,2,3,3}
    319:     rows          = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
    320:     columnsforrow = {{0,0,2,2},{1,1},{4,3,3},{5,5,5}}
    321:     vscaleforrow  = {{,,,},{,},{,,},{,,}}
    322:     vwscale       = {dx(0),dx(1),dx(2),dx(3)}               MPI Vec
    323:     vscale        = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)}   Seq Vec

    325:     ncolumns      = {1,0,1,1}
    326:     columns       = {{6},{},{4},{5}}
    327:     nrows         = {3,0,2,2}
    328:     rows          = {{0,1,2},{},{1,2},{1,2}}
    329:     columnsforrow = {{6,0,6},{},{4,4},{5,5}}
    330:     vscaleforrow =  {{,,},{},{,},{,}}
    331:     vwscale       = {dx(4),dx(5),dx(6)}              MPI Vec
    332:     vscale        = {dx(0),dx(4),dx(5),dx(6)}        Seq Vec

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

    337: */

    339: struct  _p_MatFDColoring{
    340:   PETSCHEADER(int);
    341:   PetscInt       M,N,m;            /* total rows, columns; local rows */
    342:   PetscInt       rstart;           /* first row owned by local processor */
    343:   PetscInt       ncolors;          /* number of colors */
    344:   PetscInt       *ncolumns;        /* number of local columns for a color */
    345:   PetscInt       **columns;        /* lists the local columns of each color (using global column numbering) */
    346:   PetscInt       *nrows;           /* number of local rows for each color */
    347:   PetscInt       **rows;           /* lists the local rows for each color (using the local row numbering) */
    348:   PetscInt       **columnsforrow;  /* lists the corresponding columns for those rows (using the global column) */
    349:   PetscReal      error_rel;        /* square root of relative error in computing function */
    350:   PetscReal      umin;             /* minimum allowable u'dx value */
    351:   PetscInt       freq;             /* frequency at which new Jacobian is computed */
    352:   Vec            w1,w2,w3;         /* work vectors used in computing Jacobian */
    353:   PetscErrorCode (*f)(void);       /* function that defines Jacobian */
    354:   void           *fctx;            /* optional user-defined context for use by the function f */
    355:   PetscInt       **vscaleforrow;   /* location in vscale for each columnsforrow[] entry */
    356:   Vec            vscale;           /* holds FD scaling, i.e. 1/dx for each perturbed column */
    357:   PetscTruth     usersetsrecompute;/* user determines when Jacobian is recomputed, via MatFDColoringSetRecompute() */
    358:   PetscTruth     recompute;        /* used with usersetrecompute to determine if Jacobian should be recomputed */
    359:   Vec            F;                /* current value of user provided function; can set with MatFDColoringSetF() */
    360:   PetscInt       currentcolor;     /* color for which function evaluation is being done now */
    361:   const char     *htype;            /* "wp" or "ds" */
    362:   ISColoringType ctype;            /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
    363: };

    365: /*
    366:    Null space context for preconditioner/operators
    367: */
    368: struct _p_MatNullSpace {
    369:   PETSCHEADER(int);
    370:   PetscTruth     has_cnst;
    371:   PetscInt       n;
    372:   Vec*           vecs;
    373:   PetscScalar*   alpha;                 /* for projections */
    374:   Vec            vec;                   /* for out of place removals */
    375:   PetscErrorCode (*remove)(Vec,void*);  /* for user provided removal function */
    376:   void*          rmctx;                 /* context for remove() function */
    377: };

    379: /*
    380:    Checking zero pivot for LU, ILU preconditioners.
    381: */
    382: typedef struct {
    383:   PetscInt       nshift,nshift_max;
    384:   PetscReal      shift_amount,shift_lo,shift_hi,shift_top;
    385:   PetscTruth     lushift;
    386:   PetscReal      rs;  /* active row sum of abs(offdiagonals) */
    387:   PetscScalar    pv;  /* pivot of the active row */
    388: } LUShift_Ctx;

    390: EXTERN PetscErrorCode MatFactorDumpMatrix(Mat);

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

    397:    Collective on Mat

    399:    Input Parameters:
    400: +  info - information about the matrix factorization
    401: .  sctx - pointer to the struct LUShift_Ctx
    402: -  row  - active row index

    404:    Output  Parameter:
    405: +  newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot

    407:    Level: developer
    408: @*/
    409: #define MatLUCheckShift_inline(info,sctx,row,newshift) 0;\
    410: {\
    411:   PetscInt  _newshift;\
    412:   PetscReal _rs   = sctx.rs;\
    413:   PetscReal _zero = info->zeropivot*_rs;\
    414:   if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
    415:     /* force |diag| > zeropivot*rs */\
    416:     if (!sctx.nshift){\
    417:       sctx.shift_amount = info->shiftnz;\
    418:     } else {\
    419:       sctx.shift_amount *= 2.0;\
    420:     }\
    421:     sctx.lushift = PETSC_TRUE;\
    422:     (sctx.nshift)++;\
    423:     _newshift = 1;\
    424:   } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
    425:     /* force matfactor to be diagonally dominant */\
    426:     if (sctx.nshift > sctx.nshift_max) {\
    427:       MatFactorDumpMatrix(A);\
    428:       SETERRQ1(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner after %d tries",sctx.nshift);\
    429:     } else if (sctx.nshift == sctx.nshift_max) {\
    430:       info->shift_fraction = sctx.shift_hi;\
    431:       sctx.lushift        = PETSC_TRUE;\
    432:     } else {\
    433:       sctx.shift_lo = info->shift_fraction;\
    434:       info->shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;\
    435:       sctx.lushift  = PETSC_TRUE;\
    436:     }\
    437:     sctx.shift_amount = info->shift_fraction * sctx.shift_top;\
    438:     sctx.nshift++;\
    439:     _newshift = 1;\
    440:   } else if (PetscAbsScalar(sctx.pv) <= _zero){\
    441:     MatFactorDumpMatrix(A);\
    442:     SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
    443:   } else {\
    444:     _newshift = 0;\
    445:   }\
    446:   newshift = _newshift;\
    447: }

    449: /*
    450:    Checking zero pivot for Cholesky, ICC preconditioners.
    451: */
    452: typedef struct {
    453:   PetscInt       nshift;
    454:   PetscReal      shift_amount;
    455:   PetscTruth     chshift;
    456:   PetscReal      rs;  /* active row sum of abs(offdiagonals) */
    457:   PetscScalar    pv;  /* pivot of the active row */
    458: } ChShift_Ctx;

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

    465:    Collective on Mat

    467:    Input Parameters:
    468: +  info - information about the matrix factorization
    469: .  sctx - pointer to the struct CholeskyShift_Ctx
    470: .  row  - pivot row
    471: -  newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot

    473:    Level: developer
    474:    Note: Unlike in the ILU case there is no exit condition on nshift:
    475:        we increase the shift until it converges. There is no guarantee that
    476:        this algorithm converges faster or slower, or is better or worse
    477:        than the ILU algorithm.
    478: @*/
    479: #define MatCholeskyCheckShift_inline(info,sctx,row,newshift) 0;        \
    480: {\
    481:   PetscInt  _newshift;\
    482:   PetscReal _rs   = sctx.rs;\
    483:   PetscReal _zero = info->zeropivot*_rs;\
    484:   if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
    485:     /* force |diag| > zeropivot*sctx.rs */\
    486:     if (!sctx.nshift){\
    487:       sctx.shift_amount = info->shiftnz;\
    488:     } else {\
    489:       sctx.shift_amount *= 2.0;\
    490:     }\
    491:     sctx.chshift = PETSC_TRUE;\
    492:     sctx.nshift++;\
    493:     _newshift = 1;\
    494:   } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
    495:     /* calculate a shift that would make this row diagonally dominant */\
    496:     sctx.shift_amount = PetscMax(_rs+PetscAbs(PetscRealPart(sctx.pv)),1.1*sctx.shift_amount);\
    497:     sctx.chshift      = PETSC_TRUE;\
    498:     sctx.nshift++;\
    499:     _newshift = 1;\
    500:   } else if (PetscAbsScalar(sctx.pv) <= _zero){\
    501:     SETERRQ4(PETSC_ERR_MAT_CH_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
    502:   } else {\
    503:     _newshift = 0; \
    504:   }\
    505:   newshift = _newshift;\
    506: }

    508: /*
    509:   Create and initialize a linked list
    510:   Input Parameters:
    511:     idx_start - starting index of the list
    512:     lnk_max   - max value of lnk indicating the end of the list
    513:     nlnk      - max length of the list
    514:   Output Parameters:
    515:     lnk       - list initialized
    516:     bt        - PetscBT (bitarray) with all bits set to false
    517: */
    518: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
    519:   (PetscMalloc(nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,0))

    521: /*
    522:   Add an index set into a sorted linked list
    523:   Input Parameters:
    524:     nidx      - number of input indices
    525:     indices   - interger array
    526:     idx_start - starting index of the list
    527:     lnk       - linked list(an integer array) that is created
    528:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    529:   output Parameters:
    530:     nlnk      - number of newly added indices
    531:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from indices
    532:     bt        - updated PetscBT (bitarray)
    533: */
    534: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
    535: {\
    536:   PetscInt _k,_entry,_location,_lnkdata;\
    537:   nlnk     = 0;\
    538:   _lnkdata = idx_start;\
    539:   for (_k=0; _k<nidx; _k++){\
    540:     _entry = indices[_k];\
    541:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    542:       /* search for insertion location */\
    543:       /* start from the beginning if _entry < previous _entry */\
    544:       if (_k && _entry < _lnkdata) _lnkdata  = idx_start;\
    545:       do {\
    546:         _location = _lnkdata;\
    547:         _lnkdata  = lnk[_location];\
    548:       } while (_entry > _lnkdata);\
    549:       /* insertion location is found, add entry into lnk */\
    550:       lnk[_location] = _entry;\
    551:       lnk[_entry]    = _lnkdata;\
    552:       nlnk++;\
    553:       _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
    554:     }\
    555:   }\
    556: }

    558: /*
    559:   Add a permuted index set into a sorted linked list
    560:   Input Parameters:
    561:     nidx      - number of input indices
    562:     indices   - interger array
    563:     perm      - permutation of indices
    564:     idx_start - starting index of the list
    565:     lnk       - linked list(an integer array) that is created
    566:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    567:   output Parameters:
    568:     nlnk      - number of newly added indices
    569:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from indices
    570:     bt        - updated PetscBT (bitarray)
    571: */
    572: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
    573: {\
    574:   PetscInt _k,_entry,_location,_lnkdata;\
    575:   nlnk     = 0;\
    576:   _lnkdata = idx_start;\
    577:   for (_k=0; _k<nidx; _k++){\
    578:     _entry = perm[indices[_k]];\
    579:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    580:       /* search for insertion location */\
    581:       /* start from the beginning if _entry < previous _entry */\
    582:       if (_k && _entry < _lnkdata) _lnkdata  = idx_start;\
    583:       do {\
    584:         _location = _lnkdata;\
    585:         _lnkdata  = lnk[_location];\
    586:       } while (_entry > _lnkdata);\
    587:       /* insertion location is found, add entry into lnk */\
    588:       lnk[_location] = _entry;\
    589:       lnk[_entry]    = _lnkdata;\
    590:       nlnk++;\
    591:       _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
    592:     }\
    593:   }\
    594: }

    596: /*
    597:   Add a SORTED index set into a sorted linked list
    598:   Input Parameters:
    599:     nidx      - number of input indices
    600:     indices   - sorted interger array
    601:     idx_start - starting index of the list
    602:     lnk       - linked list(an integer array) that is created
    603:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    604:   output Parameters:
    605:     nlnk      - number of newly added indices
    606:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from indices
    607:     bt        - updated PetscBT (bitarray)
    608: */
    609: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
    610: {\
    611:   PetscInt _k,_entry,_location,_lnkdata;\
    612:   nlnk      = 0;\
    613:   _lnkdata  = idx_start;\
    614:   for (_k=0; _k<nidx; _k++){\
    615:     _entry = indices[_k];\
    616:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    617:       /* search for insertion location */\
    618:       do {\
    619:         _location = _lnkdata;\
    620:         _lnkdata  = lnk[_location];\
    621:       } while (_entry > _lnkdata);\
    622:       /* insertion location is found, add entry into lnk */\
    623:       lnk[_location] = _entry;\
    624:       lnk[_entry]    = _lnkdata;\
    625:       nlnk++;\
    626:       _lnkdata = _entry; /* next search starts from here */\
    627:     }\
    628:   }\
    629: }

    631: /*
    632:   Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
    633:   Same as PetscLLAddSorted() with an additional operation:
    634:        count the number of input indices that are no larger than 'diag'
    635:   Input Parameters:
    636:     indices   - sorted interger array
    637:     idx_start - starting index of the list
    638:     lnk       - linked list(an integer array) that is created
    639:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    640:     diag      - index of the active row in LUFactorSymbolic
    641:     nzbd      - number of input indices with indices <= idx_start
    642:   output Parameters:
    643:     nlnk      - number of newly added indices
    644:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from indices
    645:     bt        - updated PetscBT (bitarray)
    646:     im        - im[idx_start] =  num of entries with indices <= diag
    647: */
    648: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
    649: {\
    650:   PetscInt _k,_entry,_location,_lnkdata,_nidx;\
    651:   nlnk     = 0;\
    652:   _lnkdata = idx_start;\
    653:   _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
    654:   for (_k=0; _k<_nidx; _k++){\
    655:     _entry = indices[_k];\
    656:     nzbd++;\
    657:     if ( _entry== diag) im[idx_start] = nzbd;\
    658:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    659:       /* search for insertion location */\
    660:       do {\
    661:         _location = _lnkdata;\
    662:         _lnkdata  = lnk[_location];\
    663:       } while (_entry > _lnkdata);\
    664:       /* insertion location is found, add entry into lnk */\
    665:       lnk[_location] = _entry;\
    666:       lnk[_entry]    = _lnkdata;\
    667:       nlnk++;\
    668:       _lnkdata = _entry; /* next search starts from here */\
    669:     }\
    670:   }\
    671: }

    673: /*
    674:   Copy data on the list into an array, then initialize the list
    675:   Input Parameters:
    676:     idx_start - starting index of the list
    677:     lnk_max   - max value of lnk indicating the end of the list
    678:     nlnk      - number of data on the list to be copied
    679:     lnk       - linked list
    680:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    681:   output Parameters:
    682:     indices   - array that contains the copied data
    683:     lnk       - linked list that is cleaned and initialize
    684:     bt        - PetscBT (bitarray) with all bits set to false
    685: */
    686: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
    687: {\
    688:   PetscInt _j,_idx=idx_start;\
    689:   for (_j=0; _j<nlnk; _j++){\
    690:     _idx = lnk[_idx];\
    691:     *(indices+_j) = _idx;\
    692:     PetscBTClear(bt,_idx);\
    693:   }\
    694:   lnk[idx_start] = lnk_max;\
    695: }
    696: /*
    697:   Free memories used by the list
    698: */
    699: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))

    701: /* Routines below are used for incomplete matrix factorization */
    702: /*
    703:   Create and initialize a linked list and its levels
    704:   Input Parameters:
    705:     idx_start - starting index of the list
    706:     lnk_max   - max value of lnk indicating the end of the list
    707:     nlnk      - max length of the list
    708:   Output Parameters:
    709:     lnk       - list initialized
    710:     lnk_lvl   - array of size nlnk for storing levels of lnk
    711:     bt        - PetscBT (bitarray) with all bits set to false
    712: */
    713: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
    714:   (PetscMalloc(2*nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))

    716: /*
    717:   Initialize a sorted linked list used for ILU and ICC
    718:   Input Parameters:
    719:     nidx      - number of input idx
    720:     idx       - interger array used for storing column indices
    721:     idx_start - starting index of the list
    722:     perm      - indices of an IS
    723:     lnk       - linked list(an integer array) that is created
    724:     lnklvl    - levels of lnk
    725:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    726:   output Parameters:
    727:     nlnk     - number of newly added idx
    728:     lnk      - the sorted(increasing order) linked list containing new and non-redundate entries from idx
    729:     lnklvl   - levels of lnk
    730:     bt       - updated PetscBT (bitarray)
    731: */
    732: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
    733: {\
    734:   PetscInt _k,_entry,_location,_lnkdata;\
    735:   nlnk     = 0;\
    736:   _lnkdata = idx_start;\
    737:   for (_k=0; _k<nidx; _k++){\
    738:     _entry = perm[idx[_k]];\
    739:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    740:       /* search for insertion location */\
    741:       if (_k && _entry < _lnkdata) _lnkdata  = idx_start;\
    742:       do {\
    743:         _location = _lnkdata;\
    744:         _lnkdata  = lnk[_location];\
    745:       } while (_entry > _lnkdata);\
    746:       /* insertion location is found, add entry into lnk */\
    747:       lnk[_location]  = _entry;\
    748:       lnk[_entry]     = _lnkdata;\
    749:       lnklvl[_entry] = 0;\
    750:       nlnk++;\
    751:       _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
    752:     }\
    753:   }\
    754: }

    756: /*
    757:   Add a SORTED index set into a sorted linked list for ILU
    758:   Input Parameters:
    759:     nidx      - number of input indices
    760:     idx       - sorted interger array used for storing column indices
    761:     level     - level of fill, e.g., ICC(level)
    762:     idxlvl    - level of idx
    763:     idx_start - starting index of the list
    764:     lnk       - linked list(an integer array) that is created
    765:     lnklvl    - levels of lnk
    766:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    767:     prow      - the row number of idx
    768:   output Parameters:
    769:     nlnk     - number of newly added idx
    770:     lnk      - the sorted(increasing order) linked list containing new and non-redundate entries from idx
    771:     lnklvl   - levels of lnk
    772:     bt       - updated PetscBT (bitarray)

    774:   Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
    775:         where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
    776: */
    777: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
    778: {\
    779:   PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
    780:   nlnk     = 0;\
    781:   _lnkdata = idx_start;\
    782:   for (_k=0; _k<nidx; _k++){\
    783:     _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
    784:     if (_incrlev > level) continue;\
    785:     _entry = idx[_k];\
    786:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    787:       /* search for insertion location */\
    788:       do {\
    789:         _location = _lnkdata;\
    790:         _lnkdata  = lnk[_location];\
    791:       } while (_entry > _lnkdata);\
    792:       /* insertion location is found, add entry into lnk */\
    793:       lnk[_location]  = _entry;\
    794:       lnk[_entry]     = _lnkdata;\
    795:       lnklvl[_entry] = _incrlev;\
    796:       nlnk++;\
    797:       _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
    798:     } else { /* existing entry: update lnklvl */\
    799:       if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
    800:     }\
    801:   }\
    802: }

    804: /*
    805:   Add a index set into a sorted linked list
    806:   Input Parameters:
    807:     nidx      - number of input idx
    808:     idx   - interger array used for storing column indices
    809:     level     - level of fill, e.g., ICC(level)
    810:     idxlvl - level of idx
    811:     idx_start - starting index of the list
    812:     lnk       - linked list(an integer array) that is created
    813:     lnklvl   - levels of lnk
    814:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    815:   output Parameters:
    816:     nlnk      - number of newly added idx
    817:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from idx
    818:     lnklvl   - levels of lnk
    819:     bt        - updated PetscBT (bitarray)
    820: */
    821: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
    822: {\
    823:   PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
    824:   nlnk     = 0;\
    825:   _lnkdata = idx_start;\
    826:   for (_k=0; _k<nidx; _k++){\
    827:     _incrlev = idxlvl[_k] + 1;\
    828:     if (_incrlev > level) continue;\
    829:     _entry = idx[_k];\
    830:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    831:       /* search for insertion location */\
    832:       if (_k && _entry < _lnkdata) _lnkdata  = idx_start;\
    833:       do {\
    834:         _location = _lnkdata;\
    835:         _lnkdata  = lnk[_location];\
    836:       } while (_entry > _lnkdata);\
    837:       /* insertion location is found, add entry into lnk */\
    838:       lnk[_location]  = _entry;\
    839:       lnk[_entry]     = _lnkdata;\
    840:       lnklvl[_entry] = _incrlev;\
    841:       nlnk++;\
    842:       _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
    843:     } else { /* existing entry: update lnklvl */\
    844:       if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
    845:     }\
    846:   }\
    847: }

    849: /*
    850:   Add a SORTED index set into a sorted linked list
    851:   Input Parameters:
    852:     nidx      - number of input indices
    853:     idx   - sorted interger array used for storing column indices
    854:     level     - level of fill, e.g., ICC(level)
    855:     idxlvl - level of idx
    856:     idx_start - starting index of the list
    857:     lnk       - linked list(an integer array) that is created
    858:     lnklvl    - levels of lnk
    859:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    860:   output Parameters:
    861:     nlnk      - number of newly added idx
    862:     lnk       - the sorted(increasing order) linked list containing new and non-redundate entries from idx
    863:     lnklvl    - levels of lnk
    864:     bt        - updated PetscBT (bitarray)
    865: */
    866: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
    867: {\
    868:   PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
    869:   nlnk = 0;\
    870:   _lnkdata = idx_start;\
    871:   for (_k=0; _k<nidx; _k++){\
    872:     _incrlev = idxlvl[_k] + 1;\
    873:     if (_incrlev > level) continue;\
    874:     _entry = idx[_k];\
    875:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    876:       /* search for insertion location */\
    877:       do {\
    878:         _location = _lnkdata;\
    879:         _lnkdata  = lnk[_location];\
    880:       } while (_entry > _lnkdata);\
    881:       /* insertion location is found, add entry into lnk */\
    882:       lnk[_location] = _entry;\
    883:       lnk[_entry]    = _lnkdata;\
    884:       lnklvl[_entry] = _incrlev;\
    885:       nlnk++;\
    886:       _lnkdata = _entry; /* next search starts from here */\
    887:     } else { /* existing entry: update lnklvl */\
    888:       if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
    889:     }\
    890:   }\
    891: }

    893: /*
    894:   Add a SORTED index set into a sorted linked list for ICC
    895:   Input Parameters:
    896:     nidx      - number of input indices
    897:     idx       - sorted interger array used for storing column indices
    898:     level     - level of fill, e.g., ICC(level)
    899:     idxlvl    - level of idx
    900:     idx_start - starting index of the list
    901:     lnk       - linked list(an integer array) that is created
    902:     lnklvl    - levels of lnk
    903:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    904:     idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
    905:   output Parameters:
    906:     nlnk   - number of newly added indices
    907:     lnk    - the sorted(increasing order) linked list containing new and non-redundate entries from idx
    908:     lnklvl - levels of lnk
    909:     bt     - updated PetscBT (bitarray)
    910:   Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
    911:         where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
    912: */
    913: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
    914: {\
    915:   PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
    916:   nlnk = 0;\
    917:   _lnkdata = idx_start;\
    918:   for (_k=0; _k<nidx; _k++){\
    919:     _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
    920:     if (_incrlev > level) continue;\
    921:     _entry = idx[_k];\
    922:     if (!PetscBTLookupSet(bt,_entry)){  /* new entry */\
    923:       /* search for insertion location */\
    924:       do {\
    925:         _location = _lnkdata;\
    926:         _lnkdata  = lnk[_location];\
    927:       } while (_entry > _lnkdata);\
    928:       /* insertion location is found, add entry into lnk */\
    929:       lnk[_location] = _entry;\
    930:       lnk[_entry]    = _lnkdata;\
    931:       lnklvl[_entry] = _incrlev;\
    932:       nlnk++;\
    933:       _lnkdata = _entry; /* next search starts from here */\
    934:     } else { /* existing entry: update lnklvl */\
    935:       if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
    936:     }\
    937:   }\
    938: }

    940: /*
    941:   Copy data on the list into an array, then initialize the list
    942:   Input Parameters:
    943:     idx_start - starting index of the list
    944:     lnk_max   - max value of lnk indicating the end of the list
    945:     nlnk      - number of data on the list to be copied
    946:     lnk       - linked list
    947:     lnklvl    - level of lnk
    948:     bt        - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
    949:   output Parameters:
    950:     indices - array that contains the copied data
    951:     lnk     - linked list that is cleaned and initialize
    952:     lnklvl  - level of lnk that is reinitialized
    953:     bt      - PetscBT (bitarray) with all bits set to false
    954: */
    955: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
    956: {\
    957:   PetscInt _j,_idx=idx_start;\
    958:   for (_j=0; _j<nlnk; _j++){\
    959:     _idx = lnk[_idx];\
    960:     *(indices+_j) = _idx;\
    961:     *(indiceslvl+_j) = lnklvl[_idx];\
    962:     lnklvl[_idx] = -1;\
    963:     PetscBTClear(bt,_idx);\
    964:   }\
    965:   lnk[idx_start] = lnk_max;\
    966: }
    967: /*
    968:   Free memories used by the list
    969: */
    970: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))

    986: #endif