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