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 (*dummy)(void);
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 (*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 (*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: 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: };
153: /*
154: If you add MatOps entries above also add them to the MATOP enum
155: in include/petscmat.h and include/finclude/petscmat.h
156: */
158: /*
159: Utility private matrix routines
160: */
161: EXTERN PetscErrorCode MatConvert_Basic(Mat, MatType,MatReuse,Mat*);
162: EXTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
163: EXTERN PetscErrorCode MatView_Private(Mat);
165: EXTERN PetscErrorCode MatHeaderCopy(Mat,Mat);
166: EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
167: EXTERN PetscErrorCode MatAXPYGetxtoy_Private(PetscInt,PetscInt*,PetscInt*,PetscInt*, PetscInt*,PetscInt*,PetscInt*, PetscInt**);
168: EXTERN PetscErrorCode MatPtAP_Basic(Mat,Mat,MatReuse,PetscReal,Mat*);
169: EXTERN PetscErrorCode MatDiagonalSet_Default(Mat,Vec,InsertMode);
171: /*
172: The stash is used to temporarily store inserted matrix values that
173: belong to another processor. During the assembly phase the stashed
174: values are moved to the correct processor and
175: */
177: typedef struct _MatStashSpace *PetscMatStashSpace;
179: struct _MatStashSpace {
180: PetscMatStashSpace next;
181: MatScalar *space_head,*val;
182: PetscInt *idx,*idy;
183: PetscInt total_space_size;
184: PetscInt local_used;
185: PetscInt local_remaining;
186: };
188: EXTERN PetscErrorCode PetscMatStashSpaceGet(PetscInt,PetscInt,PetscMatStashSpace *);
189: EXTERN PetscErrorCode PetscMatStashSpaceContiguous(PetscInt,PetscMatStashSpace *,PetscScalar *,PetscInt *,PetscInt *);
190: EXTERN PetscErrorCode PetscMatStashSpaceDestroy(PetscMatStashSpace);
192: typedef struct {
193: PetscInt nmax; /* maximum stash size */
194: PetscInt umax; /* user specified max-size */
195: PetscInt oldnmax; /* the nmax value used previously */
196: PetscInt n; /* stash size */
197: PetscInt bs; /* block size of the stash */
198: PetscInt reallocs; /* preserve the no of mallocs invoked */
199: PetscMatStashSpace space_head,space; /* linked list to hold stashed global row/column numbers and matrix values */
200: /* The following variables are used for communication */
201: MPI_Comm comm;
202: PetscMPIInt size,rank;
203: PetscMPIInt tag1,tag2;
204: MPI_Request *send_waits; /* array of send requests */
205: MPI_Request *recv_waits; /* array of receive requests */
206: MPI_Status *send_status; /* array of send status */
207: PetscInt nsends,nrecvs; /* numbers of sends and receives */
208: MatScalar *svalues; /* sending data */
209: MatScalar **rvalues; /* receiving data (values) */
210: PetscInt **rindices; /* receiving data (indices) */
211: PetscMPIInt *nprocs; /* tmp data used both during scatterbegin and end */
212: PetscInt nprocessed; /* number of messages already processed */
213: } MatStash;
215: EXTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
216: EXTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
217: EXTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
218: EXTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
219: EXTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
220: EXTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[]);
221: EXTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt);
222: EXTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
223: EXTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const MatScalar[],PetscInt,PetscInt,PetscInt);
224: EXTERN PetscErrorCode MatStashScatterBegin_Private(MatStash*,PetscInt*);
225: EXTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,MatScalar**,PetscInt*);
227: #define FACTOR_LU 1
228: #define FACTOR_CHOLESKY 2
230: typedef struct {
231: PetscInt dim;
232: PetscInt dims[4];
233: PetscInt starts[4];
234: PetscTruth noc; /* this is a single component problem, hence user will not set MatStencil.c */
235: } MatStencilInfo;
237: /* Info about using compressed row format */
238: typedef struct {
239: PetscTruth use;
240: PetscInt nrows; /* number of non-zero rows */
241: PetscInt *i; /* compressed row pointer */
242: PetscInt *rindex; /* compressed row index */
243: PetscTruth checked; /* if compressed row format have been checked for */
244: } Mat_CompressedRow;
245: EXTERN PetscErrorCode Mat_CheckCompressedRow(Mat,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);
247: struct _p_Mat {
248: PETSCHEADER(struct _MatOps);
249: PetscMap rmap,cmap;
250: void *data; /* implementation-specific data */
251: PetscInt factor; /* 0, FACTOR_LU, or FACTOR_CHOLESKY */
252: PetscTruth assembled; /* is the matrix assembled? */
253: PetscTruth was_assembled; /* new values inserted into assembled mat */
254: PetscInt num_ass; /* number of times matrix has been assembled */
255: PetscTruth same_nonzero; /* matrix has same nonzero pattern as previous */
256: MatInfo info; /* matrix information */
257: ISLocalToGlobalMapping mapping; /* mapping used in MatSetValuesLocal() */
258: ISLocalToGlobalMapping bmapping; /* mapping used in MatSetValuesBlockedLocal() */
259: InsertMode insertmode; /* have values been inserted in matrix or added? */
260: MatStash stash,bstash; /* used for assembling off-proc mat emements */
261: MatNullSpace nullsp;
262: PetscTruth preallocated;
263: MatStencilInfo stencil; /* information for structured grid */
264: PetscTruth symmetric,hermitian,structurally_symmetric;
265: PetscTruth symmetric_set,hermitian_set,structurally_symmetric_set; /* if true, then corresponding flag is correct*/
266: PetscTruth symmetric_eternal;
267: void *spptr; /* pointer for special library like SuperLU */
268: };
270: #define MatPreallocated(A) ((!(A)->preallocated) ? MatSetUpPreallocation(A) : 0)
273: /*
274: Object for partitioning graphs
275: */
277: typedef struct _MatPartitioningOps *MatPartitioningOps;
278: struct _MatPartitioningOps {
279: PetscErrorCode (*apply)(MatPartitioning,IS*);
280: PetscErrorCode (*setfromoptions)(MatPartitioning);
281: PetscErrorCode (*destroy)(MatPartitioning);
282: PetscErrorCode (*view)(MatPartitioning,PetscViewer);
283: };
285: struct _p_MatPartitioning {
286: PETSCHEADER(struct _MatPartitioningOps);
287: Mat adj;
288: PetscInt *vertex_weights;
289: PetscReal *part_weights;
290: PetscInt n; /* number of partitions */
291: void *data;
292: PetscInt setupcalled;
293: };
295: /*
296: MatFDColoring is used to compute Jacobian matrices efficiently
297: via coloring. The data structure is explained below in an example.
299: Color = 0 1 0 2 | 2 3 0
300: ---------------------------------------------------
301: 00 01 | 05
302: 10 11 | 14 15 Processor 0
303: 22 23 | 25
304: 32 33 |
305: ===================================================
306: | 44 45 46
307: 50 | 55 Processor 1
308: | 64 66
309: ---------------------------------------------------
311: ncolors = 4;
313: ncolumns = {2,1,1,0}
314: columns = {{0,2},{1},{3},{}}
315: nrows = {4,2,3,3}
316: rows = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
317: columnsforrow = {{0,0,2,2},{1,1},{4,3,3},{5,5,5}}
318: vscaleforrow = {{,,,},{,},{,,},{,,}}
319: vwscale = {dx(0),dx(1),dx(2),dx(3)} MPI Vec
320: vscale = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)} Seq Vec
322: ncolumns = {1,0,1,1}
323: columns = {{6},{},{4},{5}}
324: nrows = {3,0,2,2}
325: rows = {{0,1,2},{},{1,2},{1,2}}
326: columnsforrow = {{6,0,6},{},{4,4},{5,5}}
327: vscaleforrow = {{,,},{},{,},{,}}
328: vwscale = {dx(4),dx(5),dx(6)} MPI Vec
329: vscale = {dx(0),dx(4),dx(5),dx(6)} Seq Vec
331: See the routine MatFDColoringApply() for how this data is used
332: to compute the Jacobian.
334: */
336: struct _p_MatFDColoring{
337: PETSCHEADER(int);
338: PetscInt M,N,m; /* total rows, columns; local rows */
339: PetscInt rstart; /* first row owned by local processor */
340: PetscInt ncolors; /* number of colors */
341: PetscInt *ncolumns; /* number of local columns for a color */
342: PetscInt **columns; /* lists the local columns of each color (using global column numbering) */
343: PetscInt *nrows; /* number of local rows for each color */
344: PetscInt **rows; /* lists the local rows for each color (using the local row numbering) */
345: PetscInt **columnsforrow; /* lists the corresponding columns for those rows (using the global column) */
346: PetscReal error_rel; /* square root of relative error in computing function */
347: PetscReal umin; /* minimum allowable u'dx value */
348: PetscInt freq; /* frequency at which new Jacobian is computed */
349: Vec w1,w2,w3; /* work vectors used in computing Jacobian */
350: PetscErrorCode (*f)(void); /* function that defines Jacobian */
351: void *fctx; /* optional user-defined context for use by the function f */
352: PetscInt **vscaleforrow; /* location in vscale for each columnsforrow[] entry */
353: Vec vscale; /* holds FD scaling, i.e. 1/dx for each perturbed column */
354: PetscTruth usersetsrecompute;/* user determines when Jacobian is recomputed, via MatFDColoringSetRecompute() */
355: PetscTruth recompute; /* used with usersetrecompute to determine if Jacobian should be recomputed */
356: Vec F; /* current value of user provided function; can set with MatFDColoringSetF() */
357: PetscInt currentcolor; /* color for which function evaluation is being done now */
358: const char *htype; /* "wp" or "ds" */
359: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
360: };
362: /*
363: Null space context for preconditioner/operators
364: */
365: struct _p_MatNullSpace {
366: PETSCHEADER(int);
367: PetscTruth has_cnst;
368: PetscInt n;
369: Vec* vecs;
370: Vec vec; /* for out of place removals */
371: PetscErrorCode (*remove)(Vec,void*); /* for user provided removal function */
372: void* rmctx; /* context for remove() function */
373: };
375: /*
376: Checking zero pivot for LU, ILU preconditioners.
377: */
378: typedef struct {
379: PetscInt nshift,nshift_max;
380: PetscReal shift_amount,shift_lo,shift_hi,shift_top;
381: PetscTruth lushift;
382: PetscReal rs; /* active row sum of abs(offdiagonals) */
383: PetscScalar pv; /* pivot of the active row */
384: } LUShift_Ctx;
386: EXTERN PetscErrorCode MatFactorDumpMatrix(Mat);
390: /*@C
391: MatLUCheckShift_inline - shift the diagonals when zero pivot is detected on LU factor
393: Collective on Mat
395: Input Parameters:
396: + info - information about the matrix factorization
397: . sctx - pointer to the struct LUShift_Ctx
398: - row - active row index
400: Output Parameter:
401: + newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot
403: Level: developer
404: @*/
405: #define MatLUCheckShift_inline(info,sctx,row,newshift) 0;\
406: {\
407: PetscInt _newshift;\
408: PetscReal _rs = sctx.rs;\
409: PetscReal _zero = info->zeropivot*_rs;\
410: if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
411: /* force |diag| > zeropivot*rs */\
412: if (!sctx.nshift){\
413: sctx.shift_amount = info->shiftnz;\
414: } else {\
415: sctx.shift_amount *= 2.0;\
416: }\
417: sctx.lushift = PETSC_TRUE;\
418: (sctx.nshift)++;\
419: _newshift = 1;\
420: } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
421: /* force matfactor to be diagonally dominant */\
422: if (sctx.nshift > sctx.nshift_max) {\
423: MatFactorDumpMatrix(A);\
424: SETERRQ1(PETSC_ERR_CONV_FAILED,"Unable to determine shift to enforce positive definite preconditioner after %d tries",sctx.nshift);\
425: } else if (sctx.nshift == sctx.nshift_max) {\
426: info->shift_fraction = sctx.shift_hi;\
427: sctx.lushift = PETSC_TRUE;\
428: } else {\
429: sctx.shift_lo = info->shift_fraction;\
430: info->shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;\
431: sctx.lushift = PETSC_TRUE;\
432: }\
433: sctx.shift_amount = info->shift_fraction * sctx.shift_top;\
434: sctx.nshift++;\
435: _newshift = 1;\
436: } else if (PetscAbsScalar(sctx.pv) <= _zero){\
437: MatFactorDumpMatrix(A);\
438: SETERRQ4(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
439: } else {\
440: _newshift = 0;\
441: }\
442: newshift = _newshift;\
443: }
445: /*
446: Checking zero pivot for Cholesky, ICC preconditioners.
447: */
448: typedef struct {
449: PetscInt nshift;
450: PetscReal shift_amount;
451: PetscTruth chshift;
452: PetscReal rs; /* active row sum of abs(offdiagonals) */
453: PetscScalar pv; /* pivot of the active row */
454: } ChShift_Ctx;
458: /*@C
459: MatCholeskyCheckShift_inline - shift the diagonals when zero pivot is detected on Cholesky factor
461: Collective on Mat
463: Input Parameters:
464: + info - information about the matrix factorization
465: . sctx - pointer to the struct CholeskyShift_Ctx
466: . row - pivot row
467: - newshift - 0: shift is unchanged; 1: shft is updated; -1: zeropivot
469: Level: developer
470: Note: Unlike in the ILU case there is no exit condition on nshift:
471: we increase the shift until it converges. There is no guarantee that
472: this algorithm converges faster or slower, or is better or worse
473: than the ILU algorithm.
474: @*/
475: #define MatCholeskyCheckShift_inline(info,sctx,row,newshift) 0; \
476: {\
477: PetscInt _newshift;\
478: PetscReal _rs = sctx.rs;\
479: PetscReal _zero = info->zeropivot*_rs;\
480: if (info->shiftnz && PetscAbsScalar(sctx.pv) <= _zero){\
481: /* force |diag| > zeropivot*sctx.rs */\
482: if (!sctx.nshift){\
483: sctx.shift_amount = info->shiftnz;\
484: } else {\
485: sctx.shift_amount *= 2.0;\
486: }\
487: sctx.chshift = PETSC_TRUE;\
488: sctx.nshift++;\
489: _newshift = 1;\
490: } else if (info->shiftpd && PetscRealPart(sctx.pv) <= _zero){\
491: /* calculate a shift that would make this row diagonally dominant */\
492: sctx.shift_amount = PetscMax(_rs+PetscAbs(PetscRealPart(sctx.pv)),1.1*sctx.shift_amount);\
493: sctx.chshift = PETSC_TRUE;\
494: sctx.nshift++;\
495: _newshift = 1;\
496: } else if (PetscAbsScalar(sctx.pv) <= _zero){\
497: SETERRQ4(PETSC_ERR_MAT_CH_ZRPVT,"Zero pivot row %D value %G tolerance %G * rowsum %G",row,PetscAbsScalar(sctx.pv),_zero,_rs); \
498: } else {\
499: _newshift = 0; \
500: }\
501: newshift = _newshift;\
502: }
504: /*
505: Create and initialize a linked list
506: Input Parameters:
507: idx_start - starting index of the list
508: lnk_max - max value of lnk indicating the end of the list
509: nlnk - max length of the list
510: Output Parameters:
511: lnk - list initialized
512: bt - PetscBT (bitarray) with all bits set to false
513: */
514: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
515: (PetscMalloc(nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,0))
517: /*
518: Add an index set into a sorted linked list
519: Input Parameters:
520: nidx - number of input indices
521: indices - interger array
522: idx_start - starting index of the list
523: lnk - linked list(an integer array) that is created
524: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
525: output Parameters:
526: nlnk - number of newly added indices
527: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
528: bt - updated PetscBT (bitarray)
529: */
530: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
531: {\
532: PetscInt _k,_entry,_location,_lnkdata;\
533: nlnk = 0;\
534: _lnkdata = idx_start;\
535: for (_k=0; _k<nidx; _k++){\
536: _entry = indices[_k];\
537: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
538: /* search for insertion location */\
539: /* start from the beginning if _entry < previous _entry */\
540: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
541: do {\
542: _location = _lnkdata;\
543: _lnkdata = lnk[_location];\
544: } while (_entry > _lnkdata);\
545: /* insertion location is found, add entry into lnk */\
546: lnk[_location] = _entry;\
547: lnk[_entry] = _lnkdata;\
548: nlnk++;\
549: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
550: }\
551: }\
552: }
554: /*
555: Add a permuted index set into a sorted linked list
556: Input Parameters:
557: nidx - number of input indices
558: indices - interger array
559: perm - permutation of indices
560: idx_start - starting index of the list
561: lnk - linked list(an integer array) that is created
562: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
563: output Parameters:
564: nlnk - number of newly added indices
565: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
566: bt - updated PetscBT (bitarray)
567: */
568: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
569: {\
570: PetscInt _k,_entry,_location,_lnkdata;\
571: nlnk = 0;\
572: _lnkdata = idx_start;\
573: for (_k=0; _k<nidx; _k++){\
574: _entry = perm[indices[_k]];\
575: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
576: /* search for insertion location */\
577: /* start from the beginning if _entry < previous _entry */\
578: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
579: do {\
580: _location = _lnkdata;\
581: _lnkdata = lnk[_location];\
582: } while (_entry > _lnkdata);\
583: /* insertion location is found, add entry into lnk */\
584: lnk[_location] = _entry;\
585: lnk[_entry] = _lnkdata;\
586: nlnk++;\
587: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
588: }\
589: }\
590: }
592: /*
593: Add a SORTED index set into a sorted linked list
594: Input Parameters:
595: nidx - number of input indices
596: indices - sorted interger array
597: idx_start - starting index of the list
598: lnk - linked list(an integer array) that is created
599: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
600: output Parameters:
601: nlnk - number of newly added indices
602: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
603: bt - updated PetscBT (bitarray)
604: */
605: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
606: {\
607: PetscInt _k,_entry,_location,_lnkdata;\
608: nlnk = 0;\
609: _lnkdata = idx_start;\
610: for (_k=0; _k<nidx; _k++){\
611: _entry = indices[_k];\
612: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
613: /* search for insertion location */\
614: do {\
615: _location = _lnkdata;\
616: _lnkdata = lnk[_location];\
617: } while (_entry > _lnkdata);\
618: /* insertion location is found, add entry into lnk */\
619: lnk[_location] = _entry;\
620: lnk[_entry] = _lnkdata;\
621: nlnk++;\
622: _lnkdata = _entry; /* next search starts from here */\
623: }\
624: }\
625: }
627: /*
628: Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
629: Same as PetscLLAddSorted() with an additional operation:
630: count the number of input indices that are no larger than 'diag'
631: Input Parameters:
632: indices - sorted interger array
633: idx_start - starting index of the list
634: lnk - linked list(an integer array) that is created
635: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
636: diag - index of the active row in LUFactorSymbolic
637: nzbd - number of input indices with indices <= idx_start
638: output Parameters:
639: nlnk - number of newly added indices
640: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
641: bt - updated PetscBT (bitarray)
642: im - im[idx_start] = num of entries with indices <= diag
643: */
644: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
645: {\
646: PetscInt _k,_entry,_location,_lnkdata,_nidx;\
647: nlnk = 0;\
648: _lnkdata = idx_start;\
649: _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
650: for (_k=0; _k<_nidx; _k++){\
651: _entry = indices[_k];\
652: nzbd++;\
653: if ( _entry== diag) im[idx_start] = nzbd;\
654: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
655: /* search for insertion location */\
656: do {\
657: _location = _lnkdata;\
658: _lnkdata = lnk[_location];\
659: } while (_entry > _lnkdata);\
660: /* insertion location is found, add entry into lnk */\
661: lnk[_location] = _entry;\
662: lnk[_entry] = _lnkdata;\
663: nlnk++;\
664: _lnkdata = _entry; /* next search starts from here */\
665: }\
666: }\
667: }
669: /*
670: Copy data on the list into an array, then initialize the list
671: Input Parameters:
672: idx_start - starting index of the list
673: lnk_max - max value of lnk indicating the end of the list
674: nlnk - number of data on the list to be copied
675: lnk - linked list
676: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
677: output Parameters:
678: indices - array that contains the copied data
679: lnk - linked list that is cleaned and initialize
680: bt - PetscBT (bitarray) with all bits set to false
681: */
682: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
683: {\
684: PetscInt _j,_idx=idx_start;\
685: for (_j=0; _j<nlnk; _j++){\
686: _idx = lnk[_idx];\
687: *(indices+_j) = _idx;\
688: PetscBTClear(bt,_idx);\
689: }\
690: lnk[idx_start] = lnk_max;\
691: }
692: /*
693: Free memories used by the list
694: */
695: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
697: /* Routines below are used for incomplete matrix factorization */
698: /*
699: Create and initialize a linked list and its levels
700: Input Parameters:
701: idx_start - starting index of the list
702: lnk_max - max value of lnk indicating the end of the list
703: nlnk - max length of the list
704: Output Parameters:
705: lnk - list initialized
706: lnk_lvl - array of size nlnk for storing levels of lnk
707: bt - PetscBT (bitarray) with all bits set to false
708: */
709: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
710: (PetscMalloc(2*nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))
712: /*
713: Initialize a sorted linked list used for ILU and ICC
714: Input Parameters:
715: nidx - number of input idx
716: idx - interger array used for storing column indices
717: idx_start - starting index of the list
718: perm - indices of an IS
719: lnk - linked list(an integer array) that is created
720: lnklvl - levels of lnk
721: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
722: output Parameters:
723: nlnk - number of newly added idx
724: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
725: lnklvl - levels of lnk
726: bt - updated PetscBT (bitarray)
727: */
728: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
729: {\
730: PetscInt _k,_entry,_location,_lnkdata;\
731: nlnk = 0;\
732: _lnkdata = idx_start;\
733: for (_k=0; _k<nidx; _k++){\
734: _entry = perm[idx[_k]];\
735: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
736: /* search for insertion location */\
737: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
738: do {\
739: _location = _lnkdata;\
740: _lnkdata = lnk[_location];\
741: } while (_entry > _lnkdata);\
742: /* insertion location is found, add entry into lnk */\
743: lnk[_location] = _entry;\
744: lnk[_entry] = _lnkdata;\
745: lnklvl[_entry] = 0;\
746: nlnk++;\
747: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
748: }\
749: }\
750: }
752: /*
753: Add a SORTED index set into a sorted linked list for ILU
754: Input Parameters:
755: nidx - number of input indices
756: idx - sorted interger array used for storing column indices
757: level - level of fill, e.g., ICC(level)
758: idxlvl - level of idx
759: idx_start - starting index of the list
760: lnk - linked list(an integer array) that is created
761: lnklvl - levels of lnk
762: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
763: prow - the row number of idx
764: output Parameters:
765: nlnk - number of newly added idx
766: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
767: lnklvl - levels of lnk
768: bt - updated PetscBT (bitarray)
770: Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
771: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
772: */
773: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
774: {\
775: PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
776: nlnk = 0;\
777: _lnkdata = idx_start;\
778: for (_k=0; _k<nidx; _k++){\
779: _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
780: if (_incrlev > level) continue;\
781: _entry = idx[_k];\
782: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
783: /* search for insertion location */\
784: do {\
785: _location = _lnkdata;\
786: _lnkdata = lnk[_location];\
787: } while (_entry > _lnkdata);\
788: /* insertion location is found, add entry into lnk */\
789: lnk[_location] = _entry;\
790: lnk[_entry] = _lnkdata;\
791: lnklvl[_entry] = _incrlev;\
792: nlnk++;\
793: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
794: } else { /* existing entry: update lnklvl */\
795: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
796: }\
797: }\
798: }
800: /*
801: Add a index set into a sorted linked list
802: Input Parameters:
803: nidx - number of input idx
804: idx - interger array used for storing column indices
805: level - level of fill, e.g., ICC(level)
806: idxlvl - level of idx
807: idx_start - starting index of the list
808: lnk - linked list(an integer array) that is created
809: lnklvl - levels of lnk
810: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
811: output Parameters:
812: nlnk - number of newly added idx
813: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
814: lnklvl - levels of lnk
815: bt - updated PetscBT (bitarray)
816: */
817: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
818: {\
819: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
820: nlnk = 0;\
821: _lnkdata = idx_start;\
822: for (_k=0; _k<nidx; _k++){\
823: _incrlev = idxlvl[_k] + 1;\
824: if (_incrlev > level) continue;\
825: _entry = idx[_k];\
826: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
827: /* search for insertion location */\
828: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
829: do {\
830: _location = _lnkdata;\
831: _lnkdata = lnk[_location];\
832: } while (_entry > _lnkdata);\
833: /* insertion location is found, add entry into lnk */\
834: lnk[_location] = _entry;\
835: lnk[_entry] = _lnkdata;\
836: lnklvl[_entry] = _incrlev;\
837: nlnk++;\
838: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
839: } else { /* existing entry: update lnklvl */\
840: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
841: }\
842: }\
843: }
845: /*
846: Add a SORTED index set into a sorted linked list
847: Input Parameters:
848: nidx - number of input indices
849: idx - sorted interger array used for storing column indices
850: level - level of fill, e.g., ICC(level)
851: idxlvl - level of idx
852: idx_start - starting index of the list
853: lnk - linked list(an integer array) that is created
854: lnklvl - levels of lnk
855: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
856: output Parameters:
857: nlnk - number of newly added idx
858: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
859: lnklvl - levels of lnk
860: bt - updated PetscBT (bitarray)
861: */
862: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
863: {\
864: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
865: nlnk = 0;\
866: _lnkdata = idx_start;\
867: for (_k=0; _k<nidx; _k++){\
868: _incrlev = idxlvl[_k] + 1;\
869: if (_incrlev > level) continue;\
870: _entry = idx[_k];\
871: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
872: /* search for insertion location */\
873: do {\
874: _location = _lnkdata;\
875: _lnkdata = lnk[_location];\
876: } while (_entry > _lnkdata);\
877: /* insertion location is found, add entry into lnk */\
878: lnk[_location] = _entry;\
879: lnk[_entry] = _lnkdata;\
880: lnklvl[_entry] = _incrlev;\
881: nlnk++;\
882: _lnkdata = _entry; /* next search starts from here */\
883: } else { /* existing entry: update lnklvl */\
884: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
885: }\
886: }\
887: }
889: /*
890: Add a SORTED index set into a sorted linked list for ICC
891: Input Parameters:
892: nidx - number of input indices
893: idx - sorted interger array used for storing column indices
894: level - level of fill, e.g., ICC(level)
895: idxlvl - level of idx
896: idx_start - starting index of the list
897: lnk - linked list(an integer array) that is created
898: lnklvl - levels of lnk
899: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
900: idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
901: output Parameters:
902: nlnk - number of newly added indices
903: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
904: lnklvl - levels of lnk
905: bt - updated PetscBT (bitarray)
906: Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
907: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
908: */
909: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
910: {\
911: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
912: nlnk = 0;\
913: _lnkdata = idx_start;\
914: for (_k=0; _k<nidx; _k++){\
915: _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
916: if (_incrlev > level) continue;\
917: _entry = idx[_k];\
918: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
919: /* search for insertion location */\
920: do {\
921: _location = _lnkdata;\
922: _lnkdata = lnk[_location];\
923: } while (_entry > _lnkdata);\
924: /* insertion location is found, add entry into lnk */\
925: lnk[_location] = _entry;\
926: lnk[_entry] = _lnkdata;\
927: lnklvl[_entry] = _incrlev;\
928: nlnk++;\
929: _lnkdata = _entry; /* next search starts from here */\
930: } else { /* existing entry: update lnklvl */\
931: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
932: }\
933: }\
934: }
936: /*
937: Copy data on the list into an array, then initialize the list
938: Input Parameters:
939: idx_start - starting index of the list
940: lnk_max - max value of lnk indicating the end of the list
941: nlnk - number of data on the list to be copied
942: lnk - linked list
943: lnklvl - level of lnk
944: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
945: output Parameters:
946: indices - array that contains the copied data
947: lnk - linked list that is cleaned and initialize
948: lnklvl - level of lnk that is reinitialized
949: bt - PetscBT (bitarray) with all bits set to false
950: */
951: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
952: {\
953: PetscInt _j,_idx=idx_start;\
954: for (_j=0; _j<nlnk; _j++){\
955: _idx = lnk[_idx];\
956: *(indices+_j) = _idx;\
957: *(indiceslvl+_j) = lnklvl[_idx];\
958: lnklvl[_idx] = -1;\
959: PetscBTClear(bt,_idx);\
960: }\
961: lnk[idx_start] = lnk_max;\
962: }
963: /*
964: Free memories used by the list
965: */
966: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
980: #endif