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 PetscScalar[],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,const MatFactorInfo*);
33: PetscErrorCode (*choleskyfactor)(Mat,IS,const MatFactorInfo*);
34: PetscErrorCode (*sor)(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
35: PetscErrorCode (*transpose)(Mat,MatReuse,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 (*setoption)(Mat,MatOption,PetscTruth);
46: PetscErrorCode (*zeroentries)(Mat);
47: /*24*/
48: PetscErrorCode (*zerorows)(Mat,PetscInt,const PetscInt[],PetscScalar);
49: PetscErrorCode (*lufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
50: PetscErrorCode (*lufactornumeric)(Mat,Mat,const MatFactorInfo*);
51: PetscErrorCode (*choleskyfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
52: PetscErrorCode (*choleskyfactornumeric)(Mat,Mat,const MatFactorInfo*);
53: /*29*/
54: PetscErrorCode (*setuppreallocation)(Mat);
55: PetscErrorCode (*ilufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
56: PetscErrorCode (*iccfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
57: PetscErrorCode (*getarray)(Mat,PetscScalar**);
58: PetscErrorCode (*restorearray)(Mat,PetscScalar**);
59: /*34*/
60: PetscErrorCode (*duplicate)(Mat,MatDuplicateOption,Mat*);
61: PetscErrorCode (*forwardsolve)(Mat,Vec,Vec);
62: PetscErrorCode (*backwardsolve)(Mat,Vec,Vec);
63: PetscErrorCode (*ilufactor)(Mat,IS,IS,const MatFactorInfo*);
64: PetscErrorCode (*iccfactor)(Mat,IS,const MatFactorInfo*);
65: /*39*/
66: PetscErrorCode (*axpy)(Mat,PetscScalar,Mat,MatStructure);
67: PetscErrorCode (*getsubmatrices)(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *[]);
68: PetscErrorCode (*increaseoverlap)(Mat,PetscInt,IS[],PetscInt);
69: PetscErrorCode (*getvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscScalar []);
70: PetscErrorCode (*copy)(Mat,Mat,MatStructure);
71: /*44*/
72: PetscErrorCode (*getrowmax)(Mat,Vec,PetscInt[]);
73: PetscErrorCode (*scale)(Mat,PetscScalar);
74: PetscErrorCode (*shift)(Mat,PetscScalar);
75: PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
76: PetscErrorCode (*dummy)(void);
77: /*49*/
78: PetscErrorCode (*setblocksize)(Mat,PetscInt);
79: PetscErrorCode (*getrowij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
80: PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt *,PetscInt *[],PetscInt *[],PetscTruth *);
81: PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
82: PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscTruth,PetscTruth,PetscInt*,PetscInt *[],PetscInt *[],PetscTruth *);
83: /*54*/
84: PetscErrorCode (*fdcoloringcreate)(Mat,ISColoring,MatFDColoring);
85: PetscErrorCode (*coloringpatch)(Mat,PetscInt,PetscInt,ISColoringValue[],ISColoring*);
86: PetscErrorCode (*setunfactored)(Mat);
87: PetscErrorCode (*permute)(Mat,IS,IS,Mat*);
88: PetscErrorCode (*setvaluesblocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
89: /*59*/
90: PetscErrorCode (*getsubmatrix)(Mat,IS,IS,MatReuse,Mat*);
91: PetscErrorCode (*destroy)(Mat);
92: PetscErrorCode (*view)(Mat,PetscViewer);
93: PetscErrorCode (*convertfrom)(Mat, const MatType,MatReuse,Mat*);
94: PetscErrorCode (*usescaledform)(Mat,PetscTruth);
95: /*64*/
96: PetscErrorCode (*scalesystem)(Mat,Vec,Vec);
97: PetscErrorCode (*unscalesystem)(Mat,Vec,Vec);
98: PetscErrorCode (*setlocaltoglobalmapping)(Mat,ISLocalToGlobalMapping);
99: PetscErrorCode (*setvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
100: PetscErrorCode (*zerorowslocal)(Mat,PetscInt,const PetscInt[],PetscScalar);
101: /*69*/
102: PetscErrorCode (*getrowmaxabs)(Mat,Vec,PetscInt[]);
103: PetscErrorCode (*getrowminabs)(Mat,Vec,PetscInt[]);
104: PetscErrorCode (*convert)(Mat, const MatType,MatReuse,Mat*);
105: PetscErrorCode (*setcoloring)(Mat,ISColoring);
106: PetscErrorCode (*setvaluesadic)(Mat,void*);
107: /*74*/
108: PetscErrorCode (*setvaluesadifor)(Mat,PetscInt,void*);
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: /*79*/
114: PetscErrorCode (*permutesparsify)(Mat, PetscInt, double, double, IS, IS, Mat *);
115: PetscErrorCode (*mults)(Mat, Vecs, Vecs);
116: PetscErrorCode (*solves)(Mat, Vecs, Vecs);
117: PetscErrorCode (*getinertia)(Mat,PetscInt*,PetscInt*,PetscInt*);
118: PetscErrorCode (*load)(PetscViewer, const MatType,Mat*);
119: /*84*/
120: PetscErrorCode (*issymmetric)(Mat,PetscReal,PetscTruth*);
121: PetscErrorCode (*ishermitian)(Mat,PetscReal,PetscTruth*);
122: PetscErrorCode (*isstructurallysymmetric)(Mat,PetscTruth*);
123: PetscErrorCode (*dummy4)(void);
124: PetscErrorCode (*getvecs)(Mat,Vec*,Vec*);
125: /*89*/
126: PetscErrorCode (*matmult)(Mat,Mat,MatReuse,PetscReal,Mat*);
127: PetscErrorCode (*matmultsymbolic)(Mat,Mat,PetscReal,Mat*);
128: PetscErrorCode (*matmultnumeric)(Mat,Mat,Mat);
129: PetscErrorCode (*ptap)(Mat,Mat,MatReuse,PetscReal,Mat*);
130: PetscErrorCode (*ptapsymbolic)(Mat,Mat,PetscReal,Mat*); /* double dispatch wrapper routine */
131: /*94*/
132: PetscErrorCode (*ptapnumeric)(Mat,Mat,Mat); /* double dispatch wrapper routine */
133: PetscErrorCode (*matmulttranspose)(Mat,Mat,MatReuse,PetscReal,Mat*);
134: PetscErrorCode (*matmulttransposesymbolic)(Mat,Mat,PetscReal,Mat*);
135: PetscErrorCode (*matmulttransposenumeric)(Mat,Mat,Mat);
136: PetscErrorCode (*ptapsymbolic_seqaij)(Mat,Mat,PetscReal,Mat*); /* actual implememtation, A=seqaij */
137: /*99*/
138: PetscErrorCode (*ptapnumeric_seqaij)(Mat,Mat,Mat); /* actual implememtation, A=seqaij */
139: PetscErrorCode (*ptapsymbolic_mpiaij)(Mat,Mat,PetscReal,Mat*); /* actual implememtation, A=mpiaij */
140: PetscErrorCode (*ptapnumeric_mpiaij)(Mat,Mat,Mat); /* actual implememtation, A=mpiaij */
141: PetscErrorCode (*conjugate)(Mat); /* complex conjugate */
142: PetscErrorCode (*setsizes)(Mat,PetscInt,PetscInt,PetscInt,PetscInt);
143: /*104*/
144: PetscErrorCode (*setvaluesrow)(Mat,PetscInt,const PetscScalar[]);
145: PetscErrorCode (*realpart)(Mat);
146: PetscErrorCode (*imaginarypart)(Mat);
147: PetscErrorCode (*getrowuppertriangular)(Mat);
148: PetscErrorCode (*restorerowuppertriangular)(Mat);
149: /*109*/
150: PetscErrorCode (*matsolve)(Mat,Mat,Mat);
151: PetscErrorCode (*getredundantmatrix)(Mat,PetscInt,MPI_Comm,PetscInt,MatReuse,Mat*);
152: PetscErrorCode (*getrowmin)(Mat,Vec,PetscInt[]);
153: PetscErrorCode (*getcolumnvector)(Mat,Vec,PetscInt);
154: PetscErrorCode (*missingdiagonal)(Mat,PetscTruth*,PetscInt*);
155: /*114*/
156: PetscErrorCode (*getseqnonzerostructure)(Mat,Mat *);
157: PetscErrorCode (*create)(Mat);
158: PetscErrorCode (*getghosts)(Mat,PetscInt*,const PetscInt *[]);
159: PetscErrorCode (*dummy2)(void);
160: PetscErrorCode (*dummy3)(void);
161: /*119*/
162: PetscErrorCode (*multdiagonalblock)(Mat,Vec,Vec);
163: PetscErrorCode (*hermitiantranspose)(Mat,MatReuse,Mat*);
164: PetscErrorCode (*multhermitiantranspose)(Mat,Vec,Vec);
165: PetscErrorCode (*multhermitiantransposeadd)(Mat,Vec,Vec,Vec);
166: };
167: /*
168: If you add MatOps entries above also add them to the MATOP enum
169: in include/petscmat.h and include/finclude/petscmat.h
170: */
172: /*
173: Utility private matrix routines
174: */
175: EXTERN PetscErrorCode MatConvert_Basic(Mat, const MatType,MatReuse,Mat*);
176: EXTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
177: EXTERN PetscErrorCode MatView_Private(Mat);
179: EXTERN PetscErrorCode MatHeaderCopy(Mat,Mat);
180: EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
181: EXTERN PetscErrorCode MatAXPYGetxtoy_Private(PetscInt,PetscInt*,PetscInt*,PetscInt*, PetscInt*,PetscInt*,PetscInt*, PetscInt**);
182: EXTERN PetscErrorCode MatPtAP_Basic(Mat,Mat,MatReuse,PetscReal,Mat*);
183: EXTERN PetscErrorCode MatDiagonalSet_Default(Mat,Vec,InsertMode);
185: /*
186: The stash is used to temporarily store inserted matrix values that
187: belong to another processor. During the assembly phase the stashed
188: values are moved to the correct processor and
189: */
191: typedef struct _MatStashSpace *PetscMatStashSpace;
193: struct _MatStashSpace {
194: PetscMatStashSpace next;
195: PetscScalar *space_head,*val;
196: PetscInt *idx,*idy;
197: PetscInt total_space_size;
198: PetscInt local_used;
199: PetscInt local_remaining;
200: };
202: EXTERN PetscErrorCode PetscMatStashSpaceGet(PetscInt,PetscInt,PetscMatStashSpace *);
203: EXTERN PetscErrorCode PetscMatStashSpaceContiguous(PetscInt,PetscMatStashSpace *,PetscScalar *,PetscInt *,PetscInt *);
204: EXTERN PetscErrorCode PetscMatStashSpaceDestroy(PetscMatStashSpace);
206: typedef struct {
207: PetscInt nmax; /* maximum stash size */
208: PetscInt umax; /* user specified max-size */
209: PetscInt oldnmax; /* the nmax value used previously */
210: PetscInt n; /* stash size */
211: PetscInt bs; /* block size of the stash */
212: PetscInt reallocs; /* preserve the no of mallocs invoked */
213: PetscMatStashSpace space_head,space; /* linked list to hold stashed global row/column numbers and matrix values */
214: /* The following variables are used for communication */
215: MPI_Comm comm;
216: PetscMPIInt size,rank;
217: PetscMPIInt tag1,tag2;
218: MPI_Request *send_waits; /* array of send requests */
219: MPI_Request *recv_waits; /* array of receive requests */
220: MPI_Status *send_status; /* array of send status */
221: PetscInt nsends,nrecvs; /* numbers of sends and receives */
222: PetscScalar *svalues; /* sending data */
223: PetscInt *sindices;
224: PetscScalar **rvalues; /* receiving data (values) */
225: PetscInt **rindices; /* receiving data (indices) */
226: PetscInt nprocessed; /* number of messages already processed */
227: PetscMPIInt *flg_v; /* indicates what messages have arrived so far and from whom */
228: } MatStash;
230: EXTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
231: EXTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
232: EXTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
233: EXTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
234: EXTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
235: EXTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscTruth);
236: EXTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscTruth);
237: EXTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
238: EXTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
239: EXTERN PetscErrorCode MatStashScatterBegin_Private(Mat,MatStash*,PetscInt*);
240: EXTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,PetscScalar**,PetscInt*);
242: typedef struct {
243: PetscInt dim;
244: PetscInt dims[4];
245: PetscInt starts[4];
246: PetscTruth noc; /* this is a single component problem, hence user will not set MatStencil.c */
247: } MatStencilInfo;
249: /* Info about using compressed row format */
250: typedef struct {
251: PetscTruth use;
252: PetscInt nrows; /* number of non-zero rows */
253: PetscInt *i; /* compressed row pointer */
254: PetscInt *rindex; /* compressed row index */
255: PetscTruth checked; /* if compressed row format have been checked for */
256: } Mat_CompressedRow;
257: EXTERN PetscErrorCode Mat_CheckCompressedRow(Mat,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);
259: struct _p_Mat {
260: PETSCHEADER(struct _MatOps);
261: PetscLayout rmap,cmap;
262: void *data; /* implementation-specific data */
263: MatFactorType factor; /* MAT_FACTOR_LU, or MAT_FACTOR_CHOLESKY */
264: PetscTruth assembled; /* is the matrix assembled? */
265: PetscTruth was_assembled; /* new values inserted into assembled mat */
266: PetscInt num_ass; /* number of times matrix has been assembled */
267: PetscTruth same_nonzero; /* matrix has same nonzero pattern as previous */
268: MatInfo info; /* matrix information */
269: ISLocalToGlobalMapping mapping; /* mapping used in MatSetValuesLocal() */
270: ISLocalToGlobalMapping bmapping; /* mapping used in MatSetValuesBlockedLocal() */
271: InsertMode insertmode; /* have values been inserted in matrix or added? */
272: MatStash stash,bstash; /* used for assembling off-proc mat emements */
273: MatNullSpace nullsp;
274: PetscTruth preallocated;
275: MatStencilInfo stencil; /* information for structured grid */
276: PetscTruth symmetric,hermitian,structurally_symmetric;
277: PetscTruth symmetric_set,hermitian_set,structurally_symmetric_set; /* if true, then corresponding flag is correct*/
278: PetscTruth symmetric_eternal;
279: void *spptr; /* pointer for special library like SuperLU */
280: MatSolverPackage solvertype;
281: };
283: #define MatPreallocated(A) ((!(A)->preallocated) ? MatSetUpPreallocation(A) : 0)
286: /*
287: Object for partitioning graphs
288: */
290: typedef struct _MatPartitioningOps *MatPartitioningOps;
291: struct _MatPartitioningOps {
292: PetscErrorCode (*apply)(MatPartitioning,IS*);
293: PetscErrorCode (*setfromoptions)(MatPartitioning);
294: PetscErrorCode (*destroy)(MatPartitioning);
295: PetscErrorCode (*view)(MatPartitioning,PetscViewer);
296: };
298: struct _p_MatPartitioning {
299: PETSCHEADER(struct _MatPartitioningOps);
300: Mat adj;
301: PetscInt *vertex_weights;
302: PetscReal *part_weights;
303: PetscInt n; /* number of partitions */
304: void *data;
305: PetscInt setupcalled;
306: };
308: /*
309: MatFDColoring is used to compute Jacobian matrices efficiently
310: via coloring. The data structure is explained below in an example.
312: Color = 0 1 0 2 | 2 3 0
313: ---------------------------------------------------
314: 00 01 | 05
315: 10 11 | 14 15 Processor 0
316: 22 23 | 25
317: 32 33 |
318: ===================================================
319: | 44 45 46
320: 50 | 55 Processor 1
321: | 64 66
322: ---------------------------------------------------
324: ncolors = 4;
326: ncolumns = {2,1,1,0}
327: columns = {{0,2},{1},{3},{}}
328: nrows = {4,2,3,3}
329: rows = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
330: columnsforrow = {{0,0,2,2},{1,1},{4,3,3},{5,5,5}}
331: vscaleforrow = {{,,,},{,},{,,},{,,}}
332: vwscale = {dx(0),dx(1),dx(2),dx(3)} MPI Vec
333: vscale = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)} Seq Vec
335: ncolumns = {1,0,1,1}
336: columns = {{6},{},{4},{5}}
337: nrows = {3,0,2,2}
338: rows = {{0,1,2},{},{1,2},{1,2}}
339: columnsforrow = {{6,0,6},{},{4,4},{5,5}}
340: vscaleforrow = {{,,},{},{,},{,}}
341: vwscale = {dx(4),dx(5),dx(6)} MPI Vec
342: vscale = {dx(0),dx(4),dx(5),dx(6)} Seq Vec
344: See the routine MatFDColoringApply() for how this data is used
345: to compute the Jacobian.
347: */
349: struct _p_MatFDColoring{
350: PETSCHEADER(int);
351: PetscInt M,N,m; /* total rows, columns; local rows */
352: PetscInt rstart; /* first row owned by local processor */
353: PetscInt ncolors; /* number of colors */
354: PetscInt *ncolumns; /* number of local columns for a color */
355: PetscInt **columns; /* lists the local columns of each color (using global column numbering) */
356: PetscInt *nrows; /* number of local rows for each color */
357: PetscInt **rows; /* lists the local rows for each color (using the local row numbering) */
358: PetscInt **columnsforrow; /* lists the corresponding columns for those rows (using the global column) */
359: PetscReal error_rel; /* square root of relative error in computing function */
360: PetscReal umin; /* minimum allowable u'dx value */
361: Vec w1,w2,w3; /* work vectors used in computing Jacobian */
362: PetscErrorCode (*f)(void); /* function that defines Jacobian */
363: void *fctx; /* optional user-defined context for use by the function f */
364: PetscInt **vscaleforrow; /* location in vscale for each columnsforrow[] entry */
365: Vec vscale; /* holds FD scaling, i.e. 1/dx for each perturbed column */
366: Vec F; /* current value of user provided function; can set with MatFDColoringSetF() */
367: PetscInt currentcolor; /* color for which function evaluation is being done now */
368: const char *htype; /* "wp" or "ds" */
369: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
371: void *ftn_func_pointer,*ftn_func_cntx; /* serve the same purpose as *fortran_func_pointers in PETSc objects */
372: };
374: /*
375: Null space context for preconditioner/operators
376: */
377: struct _p_MatNullSpace {
378: PETSCHEADER(int);
379: PetscTruth has_cnst;
380: PetscInt n;
381: Vec* vecs;
382: PetscScalar* alpha; /* for projections */
383: Vec vec; /* for out of place removals */
384: PetscErrorCode (*remove)(MatNullSpace,Vec,void*); /* for user provided removal function */
385: void* rmctx; /* context for remove() function */
386: };
388: /*
389: Checking zero pivot for LU, ILU preconditioners.
390: */
391: typedef struct {
392: PetscInt nshift,nshift_max;
393: PetscReal shift_amount,shift_lo,shift_hi,shift_top,shift_fraction;
394: PetscTruth useshift;
395: PetscReal rs; /* active row sum of abs(offdiagonals) */
396: PetscScalar pv; /* pivot of the active row */
397: } FactorShiftCtx;
399: EXTERN PetscErrorCode MatFactorDumpMatrix(Mat);
403: /*@C
404: MatLUCheckShift_inline - shift the diagonals when zero pivot is detected on LU factor
406: Collective on Mat
408: Input Parameters:
409: + info - information about the matrix factorization
410: . sctx - pointer to the struct FactorShiftCtx
411: - row - active row index
413: Output Parameter:
414: + newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot
416: Level: developer
417: @*/
418: #define MatLUCheckShift_inline(info,sctx,row,newshift) 0;\
419: {\
420: PetscInt _newshift;\
421: PetscReal _rs = sctx.rs;\
422: PetscReal _zero = info->zeropivot*_rs;\
423: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO && PetscAbsScalar(sctx.pv) <= _zero){\
424: /* force |diag| > zeropivot*rs */\
425: if (!sctx.nshift){\
426: sctx.shift_amount = info->shiftamount;\
427: } else {\
428: sctx.shift_amount *= 2.0;\
429: }\
430: sctx.useshift = PETSC_TRUE;\
431: (sctx.nshift)++;\
432: _newshift = 1;\
433: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && PetscRealPart(sctx.pv) <= _zero){\
434: /* force matfactor to be diagonally dominant */\
435: if (sctx.nshift > sctx.nshift_max) {\
436: MatFactorDumpMatrix(A);\
437: SETERRQ1(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner after %d tries",sctx.nshift);\
438: } else if (sctx.nshift == sctx.nshift_max) {\
439: sctx.shift_fraction = sctx.shift_hi;\
440: sctx.useshift = PETSC_TRUE;\
441: } else {\
442: sctx.shift_lo = sctx.shift_fraction;\
443: sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;\
444: sctx.useshift = PETSC_TRUE;\
445: }\
446: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;\
447: sctx.nshift++;\
448: _newshift = 1;\
449: } else if (PetscAbsScalar(sctx.pv) <= _zero){\
450: MatFactorDumpMatrix(A);\
451: SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
452: } else {\
453: _newshift = 0;\
454: }\
455: newshift = _newshift;\
456: }
458: #define MatPivotCheck_nz(info,sctx,row) 0;\
459: {\
460: PetscReal _rs = sctx.rs;\
461: PetscReal _zero = info->zeropivot*_rs;\
462: if (PetscAbsScalar(sctx.pv) <= _zero){\
463: /* force |diag| > zeropivot*rs */\
464: if (!sctx.nshift){\
465: sctx.shift_amount = info->shiftamount;\
466: } else {\
467: sctx.shift_amount *= 2.0;\
468: }\
469: sctx.useshift = PETSC_TRUE;\
470: (sctx.nshift)++;\
471: break; \
472: } \
473: }
475: #define MatPivotCheck_pd(info,sctx,row) 0;\
476: {\
477: PetscReal _rs = sctx.rs;\
478: PetscReal _zero = info->zeropivot*_rs;\
479: if (PetscRealPart(sctx.pv) <= _zero){\
480: /* force matfactor to be diagonally dominant */\
481: if (sctx.nshift > sctx.nshift_max) {\
482: MatFactorDumpMatrix(A);\
483: SETERRQ1(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner after %d tries",sctx.nshift);\
484: } else if (sctx.nshift == sctx.nshift_max) {\
485: sctx.shift_fraction = sctx.shift_hi;\
486: sctx.useshift = PETSC_TRUE;\
487: } else {\
488: sctx.shift_lo = sctx.shift_fraction;\
489: sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;\
490: sctx.useshift = PETSC_TRUE;\
491: }\
492: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;\
493: sctx.nshift++;\
494: break; \
495: }\
496: }
498: #define MatPivotCheck_inblocks(info,sctx,row) 0;\
499: {\
500: PetscReal _zero = info->zeropivot;\
501: if (PetscAbsScalar(sctx.pv) <= _zero){\
502: sctx.pv += info->shiftamount;\
503: sctx.shift_amount = 0.0;\
504: sctx.nshift++;\
505: }\
506: }
508: #define MatPivotCheck_none(info,sctx,row) 0;\
509: {\
510: PetscReal _zero = info->zeropivot;\
511: if (PetscAbsScalar(sctx.pv) <= _zero){\
512: MatFactorDumpMatrix(A);\
513: SETERRQ3(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %G tolerance %G",row,PetscAbsScalar(sctx.pv),_zero); \
514: } \
515: }
517: #define MatPivotCheck(info,sctx,row) 0;\
518: {\
519: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO){\
520: MatPivotCheck_nz(info,sctx,row);\
521: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE){\
522: MatPivotCheck_pd(info,sctx,row);\
523: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS){\
524: MatPivotCheck_inblocks(info,sctx,srow);\
525: } else {\
526: MatPivotCheck_none(info,sctx,row);\
527: }\
528: }
530: /*
531: Checking zero pivot for Cholesky, ICC preconditioners.
532: */
533: typedef struct {
534: PetscInt nshift;
535: PetscReal shift_amount;
536: PetscTruth chshift;
537: PetscReal rs; /* active row sum of abs(offdiagonals) */
538: PetscScalar pv; /* pivot of the active row */
539: } ChShift_Ctx;
543: /*@C
544: MatCholeskyCheckShift_inline - shift the diagonals when zero pivot is detected on Cholesky factor
546: Collective on Mat
548: Input Parameters:
549: + info - information about the matrix factorization
550: . sctx - pointer to the struct CholeskyShift_Ctx
551: . row - pivot row
552: - newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot
554: Level: developer
555: Note: Unlike in the ILU case there is no exit condition on nshift:
556: we increase the shift until it converges. There is no guarantee that
557: this algorithm converges faster or slower, or is better or worse
558: than the ILU algorithm.
559: @*/
560: #define MatCholeskyCheckShift_inline(info,sctx,row,newshift) 0; \
561: {\
562: PetscInt _newshift;\
563: PetscReal _rs = sctx.rs;\
564: PetscReal _zero = info->zeropivot*_rs;\
565: if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO && PetscAbsScalar(sctx.pv) <= _zero){\
566: /* force |diag| > zeropivot*sctx.rs */\
567: if (!sctx.nshift){\
568: sctx.shift_amount = info->shiftamount;\
569: } else {\
570: sctx.shift_amount *= 2.0;\
571: }\
572: sctx.chshift = PETSC_TRUE;\
573: sctx.nshift++;\
574: _newshift = 1;\
575: } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && PetscRealPart(sctx.pv) <= _zero){\
576: /* calculate a shift that would make this row diagonally dominant */\
577: sctx.shift_amount = PetscMax(_rs+PetscAbs(PetscRealPart(sctx.pv)),1.1*sctx.shift_amount);\
578: sctx.chshift = PETSC_TRUE;\
579: sctx.nshift++;\
580: _newshift = 1;\
581: } else if (PetscAbsScalar(sctx.pv) <= _zero){\
582: SETERRQ4(PETSC_ERR_MAT_CH_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
583: } else {\
584: _newshift = 0; \
585: }\
586: newshift = _newshift;\
587: }
589: /*
590: Create and initialize a linked list
591: Input Parameters:
592: idx_start - starting index of the list
593: lnk_max - max value of lnk indicating the end of the list
594: nlnk - max length of the list
595: Output Parameters:
596: lnk - list initialized
597: bt - PetscBT (bitarray) with all bits set to false
598: */
599: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
600: (PetscMalloc(nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,0))
602: /*
603: Add an index set into a sorted linked list
604: Input Parameters:
605: nidx - number of input indices
606: indices - interger array
607: idx_start - starting index of the list
608: lnk - linked list(an integer array) that is created
609: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
610: output Parameters:
611: nlnk - number of newly added indices
612: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
613: bt - updated PetscBT (bitarray)
614: */
615: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
616: {\
617: PetscInt _k,_entry,_location,_lnkdata;\
618: nlnk = 0;\
619: _lnkdata = idx_start;\
620: for (_k=0; _k<nidx; _k++){\
621: _entry = indices[_k];\
622: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
623: /* search for insertion location */\
624: /* start from the beginning if _entry < previous _entry */\
625: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
626: do {\
627: _location = _lnkdata;\
628: _lnkdata = lnk[_location];\
629: } while (_entry > _lnkdata);\
630: /* insertion location is found, add entry into lnk */\
631: lnk[_location] = _entry;\
632: lnk[_entry] = _lnkdata;\
633: nlnk++;\
634: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
635: }\
636: }\
637: }
639: /*
640: Add a permuted index set into a sorted linked list
641: Input Parameters:
642: nidx - number of input indices
643: indices - interger array
644: perm - permutation of indices
645: idx_start - starting index of the list
646: lnk - linked list(an integer array) that is created
647: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
648: output Parameters:
649: nlnk - number of newly added indices
650: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
651: bt - updated PetscBT (bitarray)
652: */
653: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
654: {\
655: PetscInt _k,_entry,_location,_lnkdata;\
656: nlnk = 0;\
657: _lnkdata = idx_start;\
658: for (_k=0; _k<nidx; _k++){\
659: _entry = perm[indices[_k]];\
660: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
661: /* search for insertion location */\
662: /* start from the beginning if _entry < previous _entry */\
663: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
664: do {\
665: _location = _lnkdata;\
666: _lnkdata = lnk[_location];\
667: } while (_entry > _lnkdata);\
668: /* insertion location is found, add entry into lnk */\
669: lnk[_location] = _entry;\
670: lnk[_entry] = _lnkdata;\
671: nlnk++;\
672: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
673: }\
674: }\
675: }
677: /*
678: Add a SORTED index set into a sorted linked list
679: Input Parameters:
680: nidx - number of input indices
681: indices - sorted interger array
682: idx_start - starting index of the list
683: lnk - linked list(an integer array) that is created
684: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
685: output Parameters:
686: nlnk - number of newly added indices
687: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
688: bt - updated PetscBT (bitarray)
689: */
690: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
691: {\
692: PetscInt _k,_entry,_location,_lnkdata;\
693: nlnk = 0;\
694: _lnkdata = idx_start;\
695: for (_k=0; _k<nidx; _k++){\
696: _entry = indices[_k];\
697: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
698: /* search for insertion location */\
699: do {\
700: _location = _lnkdata;\
701: _lnkdata = lnk[_location];\
702: } while (_entry > _lnkdata);\
703: /* insertion location is found, add entry into lnk */\
704: lnk[_location] = _entry;\
705: lnk[_entry] = _lnkdata;\
706: nlnk++;\
707: _lnkdata = _entry; /* next search starts from here */\
708: }\
709: }\
710: }
712: /*
713: Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
714: Same as PetscLLAddSorted() with an additional operation:
715: count the number of input indices that are no larger than 'diag'
716: Input Parameters:
717: indices - sorted interger array
718: idx_start - starting index of the list, index of pivot row
719: lnk - linked list(an integer array) that is created
720: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
721: diag - index of the active row in LUFactorSymbolic
722: nzbd - number of input indices with indices <= idx_start
723: im - im[idx_start] is initialized as num of nonzero entries in row=idx_start
724: output Parameters:
725: nlnk - number of newly added indices
726: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
727: bt - updated PetscBT (bitarray)
728: im - im[idx_start]: unchanged if diag is not an entry
729: : num of entries with indices <= diag if diag is an entry
730: */
731: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
732: {\
733: PetscInt _k,_entry,_location,_lnkdata,_nidx;\
734: nlnk = 0;\
735: _lnkdata = idx_start;\
736: _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
737: for (_k=0; _k<_nidx; _k++){\
738: _entry = indices[_k];\
739: nzbd++;\
740: if ( _entry== diag) im[idx_start] = nzbd;\
741: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
742: /* search for insertion location */\
743: do {\
744: _location = _lnkdata;\
745: _lnkdata = lnk[_location];\
746: } while (_entry > _lnkdata);\
747: /* insertion location is found, add entry into lnk */\
748: lnk[_location] = _entry;\
749: lnk[_entry] = _lnkdata;\
750: nlnk++;\
751: _lnkdata = _entry; /* next search starts from here */\
752: }\
753: }\
754: }
756: /*
757: Copy data on the list into an array, then initialize the list
758: Input Parameters:
759: idx_start - starting index of the list
760: lnk_max - max value of lnk indicating the end of the list
761: nlnk - number of data on the list to be copied
762: lnk - linked list
763: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
764: output Parameters:
765: indices - array that contains the copied data
766: lnk - linked list that is cleaned and initialize
767: bt - PetscBT (bitarray) with all bits set to false
768: */
769: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
770: {\
771: PetscInt _j,_idx=idx_start;\
772: for (_j=0; _j<nlnk; _j++){\
773: _idx = lnk[_idx];\
774: *(indices+_j) = _idx;\
775: PetscBTClear(bt,_idx);\
776: }\
777: lnk[idx_start] = lnk_max;\
778: }
779: /*
780: Free memories used by the list
781: */
782: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
784: /* Routines below are used for incomplete matrix factorization */
785: /*
786: Create and initialize a linked list and its levels
787: Input Parameters:
788: idx_start - starting index of the list
789: lnk_max - max value of lnk indicating the end of the list
790: nlnk - max length of the list
791: Output Parameters:
792: lnk - list initialized
793: lnk_lvl - array of size nlnk for storing levels of lnk
794: bt - PetscBT (bitarray) with all bits set to false
795: */
796: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
797: (PetscMalloc(2*nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))
799: /*
800: Initialize a sorted linked list used for ILU and ICC
801: Input Parameters:
802: nidx - number of input idx
803: idx - interger array used for storing column indices
804: idx_start - starting index of the list
805: perm - indices of an IS
806: lnk - linked list(an integer array) that is created
807: lnklvl - levels of lnk
808: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
809: output Parameters:
810: nlnk - number of newly added idx
811: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
812: lnklvl - levels of lnk
813: bt - updated PetscBT (bitarray)
814: */
815: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
816: {\
817: PetscInt _k,_entry,_location,_lnkdata;\
818: nlnk = 0;\
819: _lnkdata = idx_start;\
820: for (_k=0; _k<nidx; _k++){\
821: _entry = perm[idx[_k]];\
822: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
823: /* search for insertion location */\
824: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
825: do {\
826: _location = _lnkdata;\
827: _lnkdata = lnk[_location];\
828: } while (_entry > _lnkdata);\
829: /* insertion location is found, add entry into lnk */\
830: lnk[_location] = _entry;\
831: lnk[_entry] = _lnkdata;\
832: lnklvl[_entry] = 0;\
833: nlnk++;\
834: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
835: }\
836: }\
837: }
839: /*
840: Add a SORTED index set into a sorted linked list for ILU
841: Input Parameters:
842: nidx - number of input indices
843: idx - sorted interger array used for storing column indices
844: level - level of fill, e.g., ICC(level)
845: idxlvl - level of idx
846: idx_start - starting index of the list
847: lnk - linked list(an integer array) that is created
848: lnklvl - levels of lnk
849: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
850: prow - the row number of idx
851: output Parameters:
852: nlnk - number of newly added idx
853: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
854: lnklvl - levels of lnk
855: bt - updated PetscBT (bitarray)
857: Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
858: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
859: */
860: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
861: {\
862: PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
863: nlnk = 0;\
864: _lnkdata = idx_start;\
865: for (_k=0; _k<nidx; _k++){\
866: _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
867: if (_incrlev > level) continue;\
868: _entry = idx[_k];\
869: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
870: /* search for insertion location */\
871: do {\
872: _location = _lnkdata;\
873: _lnkdata = lnk[_location];\
874: } while (_entry > _lnkdata);\
875: /* insertion location is found, add entry into lnk */\
876: lnk[_location] = _entry;\
877: lnk[_entry] = _lnkdata;\
878: lnklvl[_entry] = _incrlev;\
879: nlnk++;\
880: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
881: } else { /* existing entry: update lnklvl */\
882: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
883: }\
884: }\
885: }
887: /*
888: Add a index set into a sorted linked list
889: Input Parameters:
890: nidx - number of input idx
891: idx - interger array used for storing column indices
892: level - level of fill, e.g., ICC(level)
893: idxlvl - level of idx
894: idx_start - starting index of the list
895: lnk - linked list(an integer array) that is created
896: lnklvl - levels of lnk
897: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
898: output Parameters:
899: nlnk - number of newly added idx
900: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
901: lnklvl - levels of lnk
902: bt - updated PetscBT (bitarray)
903: */
904: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
905: {\
906: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
907: nlnk = 0;\
908: _lnkdata = idx_start;\
909: for (_k=0; _k<nidx; _k++){\
910: _incrlev = idxlvl[_k] + 1;\
911: if (_incrlev > level) continue;\
912: _entry = idx[_k];\
913: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
914: /* search for insertion location */\
915: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
916: do {\
917: _location = _lnkdata;\
918: _lnkdata = lnk[_location];\
919: } while (_entry > _lnkdata);\
920: /* insertion location is found, add entry into lnk */\
921: lnk[_location] = _entry;\
922: lnk[_entry] = _lnkdata;\
923: lnklvl[_entry] = _incrlev;\
924: nlnk++;\
925: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
926: } else { /* existing entry: update lnklvl */\
927: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
928: }\
929: }\
930: }
932: /*
933: Add a SORTED index set into a sorted linked list
934: Input Parameters:
935: nidx - number of input indices
936: idx - sorted interger array used for storing column indices
937: level - level of fill, e.g., ICC(level)
938: idxlvl - level of idx
939: idx_start - starting index of the list
940: lnk - linked list(an integer array) that is created
941: lnklvl - levels of lnk
942: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
943: output Parameters:
944: nlnk - number of newly added idx
945: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
946: lnklvl - levels of lnk
947: bt - updated PetscBT (bitarray)
948: */
949: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
950: {\
951: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
952: nlnk = 0;\
953: _lnkdata = idx_start;\
954: for (_k=0; _k<nidx; _k++){\
955: _incrlev = idxlvl[_k] + 1;\
956: if (_incrlev > level) continue;\
957: _entry = idx[_k];\
958: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
959: /* search for insertion location */\
960: do {\
961: _location = _lnkdata;\
962: _lnkdata = lnk[_location];\
963: } while (_entry > _lnkdata);\
964: /* insertion location is found, add entry into lnk */\
965: lnk[_location] = _entry;\
966: lnk[_entry] = _lnkdata;\
967: lnklvl[_entry] = _incrlev;\
968: nlnk++;\
969: _lnkdata = _entry; /* next search starts from here */\
970: } else { /* existing entry: update lnklvl */\
971: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
972: }\
973: }\
974: }
976: /*
977: Add a SORTED index set into a sorted linked list for ICC
978: Input Parameters:
979: nidx - number of input indices
980: idx - sorted interger array used for storing column indices
981: level - level of fill, e.g., ICC(level)
982: idxlvl - level of idx
983: idx_start - starting index of the list
984: lnk - linked list(an integer array) that is created
985: lnklvl - levels of lnk
986: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
987: idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
988: output Parameters:
989: nlnk - number of newly added indices
990: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
991: lnklvl - levels of lnk
992: bt - updated PetscBT (bitarray)
993: Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
994: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
995: */
996: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
997: {\
998: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
999: nlnk = 0;\
1000: _lnkdata = idx_start;\
1001: for (_k=0; _k<nidx; _k++){\
1002: _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
1003: if (_incrlev > level) continue;\
1004: _entry = idx[_k];\
1005: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1006: /* search for insertion location */\
1007: do {\
1008: _location = _lnkdata;\
1009: _lnkdata = lnk[_location];\
1010: } while (_entry > _lnkdata);\
1011: /* insertion location is found, add entry into lnk */\
1012: lnk[_location] = _entry;\
1013: lnk[_entry] = _lnkdata;\
1014: lnklvl[_entry] = _incrlev;\
1015: nlnk++;\
1016: _lnkdata = _entry; /* next search starts from here */\
1017: } else { /* existing entry: update lnklvl */\
1018: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1019: }\
1020: }\
1021: }
1023: /*
1024: Copy data on the list into an array, then initialize the list
1025: Input Parameters:
1026: idx_start - starting index of the list
1027: lnk_max - max value of lnk indicating the end of the list
1028: nlnk - number of data on the list to be copied
1029: lnk - linked list
1030: lnklvl - level of lnk
1031: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1032: output Parameters:
1033: indices - array that contains the copied data
1034: lnk - linked list that is cleaned and initialize
1035: lnklvl - level of lnk that is reinitialized
1036: bt - PetscBT (bitarray) with all bits set to false
1037: */
1038: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
1039: {\
1040: PetscInt _j,_idx=idx_start;\
1041: for (_j=0; _j<nlnk; _j++){\
1042: _idx = lnk[_idx];\
1043: *(indices+_j) = _idx;\
1044: *(indiceslvl+_j) = lnklvl[_idx];\
1045: lnklvl[_idx] = -1;\
1046: PetscBTClear(bt,_idx);\
1047: }\
1048: lnk[idx_start] = lnk_max;\
1049: }
1050: /*
1051: Free memories used by the list
1052: */
1053: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
1072: #endif