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,PetscBool *);
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,PetscBool );
46: PetscErrorCode (*zeroentries)(Mat);
47: /*24*/
48: PetscErrorCode (*zerorows)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
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 (*zerorowscolumns)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
77: /*49*/
78: PetscErrorCode (*setblocksize)(Mat,PetscInt);
79: PetscErrorCode (*getrowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,PetscInt *[],PetscInt *[],PetscBool *);
80: PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt *,PetscInt *[],PetscInt *[],PetscBool *);
81: PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,PetscInt *[],PetscInt *[],PetscBool *);
82: PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,PetscInt *[],PetscInt *[],PetscBool *);
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,PetscBool );
95: /*64*/
96: PetscErrorCode (*scalesystem)(Mat,Vec,Vec);
97: PetscErrorCode (*unscalesystem)(Mat,Vec,Vec);
98: PetscErrorCode (*setlocaltoglobalmapping)(Mat,ISLocalToGlobalMapping,ISLocalToGlobalMapping);
99: PetscErrorCode (*setvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
100: PetscErrorCode (*zerorowslocal)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
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 (*findzerodiagonals)(Mat,IS*);
115: PetscErrorCode (*mults)(Mat, Vecs, Vecs);
116: PetscErrorCode (*solves)(Mat, Vecs, Vecs);
117: PetscErrorCode (*getinertia)(Mat,PetscInt*,PetscInt*,PetscInt*);
118: PetscErrorCode (*load)(Mat, PetscViewer);
119: /*84*/
120: PetscErrorCode (*issymmetric)(Mat,PetscReal,PetscBool *);
121: PetscErrorCode (*ishermitian)(Mat,PetscReal,PetscBool *);
122: PetscErrorCode (*isstructurallysymmetric)(Mat,PetscBool *);
123: PetscErrorCode (*setvaluesblockedlocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
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,PetscBool *,PetscInt*);
155: /*114*/
156: PetscErrorCode (*getseqnonzerostructure)(Mat,Mat *);
157: PetscErrorCode (*create)(Mat);
158: PetscErrorCode (*getghosts)(Mat,PetscInt*,const PetscInt *[]);
159: PetscErrorCode (*getlocalsubmatrix)(Mat,IS,IS,Mat*);
160: PetscErrorCode (*restorelocalsubmatrix)(Mat,IS,IS,Mat*);
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: PetscErrorCode (*getmultiprocblock)(Mat,MPI_Comm,Mat*);
167: /*124*/
168: PetscErrorCode (*findnonzerorows)(Mat,IS*);
169: PetscErrorCode (*getcolumnnorms)(Mat,NormType,PetscReal*);
170: PetscErrorCode (*invertblockdiagonal)(Mat,PetscScalar**);
171: PetscErrorCode (*dummy4)(Mat,Vec,Vec,Vec);
172: PetscErrorCode (*getsubmatricesparallel)(Mat,PetscInt,const IS[], const IS[], MatReuse, Mat**);
173: /*129*/
174: PetscErrorCode (*setvaluesbatch)(Mat,PetscInt,PetscInt,PetscInt*,const PetscScalar*);
175: };
176: /*
177: If you add MatOps entries above also add them to the MATOP enum
178: in include/petscmat.h and include/finclude/petscmat.h
179: */
181: typedef struct _p_MatBaseName* MatBaseName;
182: struct _p_MatBaseName {
183: char *bname,*sname,*mname;
184: MatBaseName next;
185: };
189: /*
190: Utility private matrix routines
191: */
202: /*
203: The stash is used to temporarily store inserted matrix values that
204: belong to another processor. During the assembly phase the stashed
205: values are moved to the correct processor and
206: */
208: typedef struct _MatStashSpace *PetscMatStashSpace;
210: struct _MatStashSpace {
211: PetscMatStashSpace next;
212: PetscScalar *space_head,*val;
213: PetscInt *idx,*idy;
214: PetscInt total_space_size;
215: PetscInt local_used;
216: PetscInt local_remaining;
217: };
223: typedef struct {
224: PetscInt nmax; /* maximum stash size */
225: PetscInt umax; /* user specified max-size */
226: PetscInt oldnmax; /* the nmax value used previously */
227: PetscInt n; /* stash size */
228: PetscInt bs; /* block size of the stash */
229: PetscInt reallocs; /* preserve the no of mallocs invoked */
230: PetscMatStashSpace space_head,space; /* linked list to hold stashed global row/column numbers and matrix values */
231: /* The following variables are used for communication */
232: MPI_Comm comm;
233: PetscMPIInt size,rank;
234: PetscMPIInt tag1,tag2;
235: MPI_Request *send_waits; /* array of send requests */
236: MPI_Request *recv_waits; /* array of receive requests */
237: MPI_Status *send_status; /* array of send status */
238: PetscInt nsends,nrecvs; /* numbers of sends and receives */
239: PetscScalar *svalues; /* sending data */
240: PetscInt *sindices;
241: PetscScalar **rvalues; /* receiving data (values) */
242: PetscInt **rindices; /* receiving data (indices) */
243: PetscInt nprocessed; /* number of messages already processed */
244: PetscMPIInt *flg_v; /* indicates what messages have arrived so far and from whom */
245: PetscBool reproduce;
246: PetscInt reproduce_count;
247: } MatStash;
261: typedef struct {
262: PetscInt dim;
263: PetscInt dims[4];
264: PetscInt starts[4];
265: PetscBool noc; /* this is a single component problem, hence user will not set MatStencil.c */
266: } MatStencilInfo;
268: /* Info about using compressed row format */
269: typedef struct {
270: PetscBool check; /* indicates that at MatAssembly() it should check if compressed rows will be efficient */
271: PetscBool use; /* indicates compressed rows have been checked and will be used */
272: PetscInt nrows; /* number of non-zero rows */
273: PetscInt *i; /* compressed row pointer */
274: PetscInt *rindex; /* compressed row index */
275: } Mat_CompressedRow;
278: struct _p_Mat {
279: PETSCHEADER(struct _MatOps);
280: PetscLayout rmap,cmap;
281: void *data; /* implementation-specific data */
282: MatFactorType factortype; /* MAT_FACTOR_LU, ILU, CHOLESKY or ICC */
283: PetscBool assembled; /* is the matrix assembled? */
284: PetscBool was_assembled; /* new values inserted into assembled mat */
285: PetscInt num_ass; /* number of times matrix has been assembled */
286: PetscBool same_nonzero; /* matrix has same nonzero pattern as previous */
287: MatInfo info; /* matrix information */
288: InsertMode insertmode; /* have values been inserted in matrix or added? */
289: MatStash stash,bstash; /* used for assembling off-proc mat emements */
290: MatNullSpace nullsp;
291: PetscBool preallocated;
292: MatStencilInfo stencil; /* information for structured grid */
293: PetscBool symmetric,hermitian,structurally_symmetric,spd;
294: PetscBool symmetric_set,hermitian_set,structurally_symmetric_set,spd_set; /* if true, then corresponding flag is correct*/
295: PetscBool symmetric_eternal;
296: PetscBool nooffprocentries,nooffproczerorows;
297: #if defined(PETSC_HAVE_CUSP)
298: PetscCUSPFlag valid_GPU_matrix; /* flag pointing to the matrix on the gpu*/
299: #endif
300: void *spptr; /* pointer for special library like SuperLU */
301: MatSolverPackage solvertype;
302: };
304: #define MatPreallocated(A) ((!(A)->preallocated) ? MatSetUpPreallocation(A) : 0)
307: /*
308: Object for partitioning graphs
309: */
311: typedef struct _MatPartitioningOps *MatPartitioningOps;
312: struct _MatPartitioningOps {
313: PetscErrorCode (*apply)(MatPartitioning,IS*);
314: PetscErrorCode (*setfromoptions)(MatPartitioning);
315: PetscErrorCode (*destroy)(MatPartitioning);
316: PetscErrorCode (*view)(MatPartitioning,PetscViewer);
317: };
319: struct _p_MatPartitioning {
320: PETSCHEADER(struct _MatPartitioningOps);
321: Mat adj;
322: PetscInt *vertex_weights;
323: PetscReal *part_weights;
324: PetscInt n; /* number of partitions */
325: void *data;
326: PetscInt setupcalled;
327: };
329: /*
330: MatFDColoring is used to compute Jacobian matrices efficiently
331: via coloring. The data structure is explained below in an example.
333: Color = 0 1 0 2 | 2 3 0
334: ---------------------------------------------------
335: 00 01 | 05
336: 10 11 | 14 15 Processor 0
337: 22 23 | 25
338: 32 33 |
339: ===================================================
340: | 44 45 46
341: 50 | 55 Processor 1
342: | 64 66
343: ---------------------------------------------------
345: ncolors = 4;
347: ncolumns = {2,1,1,0}
348: columns = {{0,2},{1},{3},{}}
349: nrows = {4,2,3,3}
350: rows = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
351: columnsforrow = {{0,0,2,2},{1,1},{4,3,3},{5,5,5}}
352: vscaleforrow = {{,,,},{,},{,,},{,,}}
353: vwscale = {dx(0),dx(1),dx(2),dx(3)} MPI Vec
354: vscale = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)} Seq Vec
356: ncolumns = {1,0,1,1}
357: columns = {{6},{},{4},{5}}
358: nrows = {3,0,2,2}
359: rows = {{0,1,2},{},{1,2},{1,2}}
360: columnsforrow = {{6,0,6},{},{4,4},{5,5}}
361: vscaleforrow = {{,,},{},{,},{,}}
362: vwscale = {dx(4),dx(5),dx(6)} MPI Vec
363: vscale = {dx(0),dx(4),dx(5),dx(6)} Seq Vec
365: See the routine MatFDColoringApply() for how this data is used
366: to compute the Jacobian.
368: */
370: struct _p_MatFDColoring{
371: PETSCHEADER(int);
372: PetscInt M,N,m; /* total rows, columns; local rows */
373: PetscInt rstart; /* first row owned by local processor */
374: PetscInt ncolors; /* number of colors */
375: PetscInt *ncolumns; /* number of local columns for a color */
376: PetscInt **columns; /* lists the local columns of each color (using global column numbering) */
377: PetscInt *nrows; /* number of local rows for each color */
378: PetscInt **rows; /* lists the local rows for each color (using the local row numbering) */
379: PetscInt **columnsforrow; /* lists the corresponding columns for those rows (using the global column) */
380: PetscReal error_rel; /* square root of relative error in computing function */
381: PetscReal umin; /* minimum allowable u'dx value */
382: Vec w1,w2,w3; /* work vectors used in computing Jacobian */
383: PetscErrorCode (*f)(void); /* function that defines Jacobian */
384: void *fctx; /* optional user-defined context for use by the function f */
385: PetscInt **vscaleforrow; /* location in vscale for each columnsforrow[] entry */
386: Vec vscale; /* holds FD scaling, i.e. 1/dx for each perturbed column */
387: Vec F; /* current value of user provided function; can set with MatFDColoringSetF() */
388: PetscInt currentcolor; /* color for which function evaluation is being done now */
389: const char *htype; /* "wp" or "ds" */
390: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
392: void *ftn_func_pointer,*ftn_func_cntx; /* serve the same purpose as *fortran_func_pointers in PETSc objects */
393: };
395: /*
396: Null space context for preconditioner/operators
397: */
398: struct _p_MatNullSpace {
399: PETSCHEADER(int);
400: PetscBool has_cnst;
401: PetscInt n;
402: Vec* vecs;
403: PetscScalar* alpha; /* for projections */
404: Vec vec; /* for out of place removals */
405: PetscErrorCode (*remove)(MatNullSpace,Vec,void*); /* for user provided removal function */
406: void* rmctx; /* context for remove() function */
407: };
409: /*
410: Checking zero pivot for LU, ILU preconditioners.
411: */
412: typedef struct {
413: PetscInt nshift,nshift_max;
414: PetscReal shift_amount,shift_lo,shift_hi,shift_top,shift_fraction;
415: PetscBool newshift;
416: PetscReal rs; /* active row sum of abs(offdiagonals) */
417: PetscScalar pv; /* pivot of the active row */
418: } FactorShiftCtx;
424: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_nz(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
425: {
426: PetscReal _rs = sctx->rs;
427: PetscReal _zero = info->zeropivot*_rs;
430: if (PetscAbsScalar(sctx->pv) <= _zero){
431: /* force |diag| > zeropivot*rs */
432: if (!sctx->nshift) sctx->shift_amount = info->shiftamount;
433: else sctx->shift_amount *= 2.0;
434: sctx->newshift = PETSC_TRUE;
435: (sctx->nshift)++;
436: } else {
437: sctx->newshift = PETSC_FALSE;
438: }
439: return(0);
440: }
444: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_pd(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
445: {
446: PetscReal _rs = sctx->rs;
447: PetscReal _zero = info->zeropivot*_rs;
450: if (PetscRealPart(sctx->pv) <= _zero){
451: /* force matfactor to be diagonally dominant */
452: if (sctx->nshift == sctx->nshift_max) {
453: sctx->shift_fraction = sctx->shift_hi;
454: } else {
455: sctx->shift_lo = sctx->shift_fraction;
456: sctx->shift_fraction = (sctx->shift_hi+sctx->shift_lo)/2.;
457: }
458: sctx->shift_amount = sctx->shift_fraction * sctx->shift_top;
459: sctx->nshift++;
460: sctx->newshift = PETSC_TRUE;
461: } else {
462: sctx->newshift = PETSC_FALSE;
463: }
464: return(0);
465: }
469: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_inblocks(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
470: {
471: PetscReal _zero = info->zeropivot;
474: if (PetscAbsScalar(sctx->pv) <= _zero){
475: sctx->pv += info->shiftamount;
476: sctx->shift_amount = 0.0;
477: sctx->nshift++;
478: }
479: sctx->newshift = PETSC_FALSE;
480: return(0);
481: }
485: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_none(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
486: {
487: PetscReal _zero = info->zeropivot;
490: sctx->newshift = PETSC_FALSE;
491: if (PetscAbsScalar(sctx->pv) <= _zero) {
493: PetscBool flg = PETSC_FALSE;
494:
495: PetscOptionsGetBool(PETSC_NULL,"-mat_dump",&flg,PETSC_NULL);
496: if (flg) {
497: MatView(mat,PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
498: }
499: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %G tolerance %G",row,PetscAbsScalar(sctx->pv),_zero);
500: }
501: return(0);
502: }
506: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
507: {
511: if (info->shifttype == (PetscReal) MAT_SHIFT_NONZERO){
512: MatPivotCheck_nz(mat,info,sctx,row);
513: } else if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE){
514: MatPivotCheck_pd(mat,info,sctx,row);
515: } else if (info->shifttype == (PetscReal) MAT_SHIFT_INBLOCKS){
516: MatPivotCheck_inblocks(mat,info,sctx,row);
517: } else {
518: MatPivotCheck_none(mat,info,sctx,row);
519: }
520: return(0);
521: }
523: /*
524: Create and initialize a linked list
525: Input Parameters:
526: idx_start - starting index of the list
527: lnk_max - max value of lnk indicating the end of the list
528: nlnk - max length of the list
529: Output Parameters:
530: lnk - list initialized
531: bt - PetscBT (bitarray) with all bits set to false
532: */
533: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
534: (PetscMalloc(nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,0))
536: /*
537: Add an index set into a sorted linked list
538: Input Parameters:
539: nidx - number of input indices
540: indices - interger array
541: idx_start - starting index of the list
542: lnk - linked list(an integer array) that is created
543: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
544: output Parameters:
545: nlnk - number of newly added indices
546: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
547: bt - updated PetscBT (bitarray)
548: */
549: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
550: {\
551: PetscInt _k,_entry,_location,_lnkdata;\
552: nlnk = 0;\
553: _lnkdata = idx_start;\
554: for (_k=0; _k<nidx; _k++){\
555: _entry = indices[_k];\
556: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
557: /* search for insertion location */\
558: /* start from the beginning if _entry < previous _entry */\
559: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
560: do {\
561: _location = _lnkdata;\
562: _lnkdata = lnk[_location];\
563: } while (_entry > _lnkdata);\
564: /* insertion location is found, add entry into lnk */\
565: lnk[_location] = _entry;\
566: lnk[_entry] = _lnkdata;\
567: nlnk++;\
568: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
569: }\
570: }\
571: }
573: /*
574: Add a permuted index set into a sorted linked list
575: Input Parameters:
576: nidx - number of input indices
577: indices - interger array
578: perm - permutation of indices
579: idx_start - starting index of the list
580: lnk - linked list(an integer array) that is created
581: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
582: output Parameters:
583: nlnk - number of newly added indices
584: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
585: bt - updated PetscBT (bitarray)
586: */
587: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
588: {\
589: PetscInt _k,_entry,_location,_lnkdata;\
590: nlnk = 0;\
591: _lnkdata = idx_start;\
592: for (_k=0; _k<nidx; _k++){\
593: _entry = perm[indices[_k]];\
594: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
595: /* search for insertion location */\
596: /* start from the beginning if _entry < previous _entry */\
597: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
598: do {\
599: _location = _lnkdata;\
600: _lnkdata = lnk[_location];\
601: } while (_entry > _lnkdata);\
602: /* insertion location is found, add entry into lnk */\
603: lnk[_location] = _entry;\
604: lnk[_entry] = _lnkdata;\
605: nlnk++;\
606: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
607: }\
608: }\
609: }
611: /*
612: Add a SORTED index set into a sorted linked list
613: Input Parameters:
614: nidx - number of input indices
615: indices - sorted interger array
616: idx_start - starting index of the list
617: lnk - linked list(an integer array) that is created
618: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
619: output Parameters:
620: nlnk - number of newly added indices
621: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
622: bt - updated PetscBT (bitarray)
623: */
624: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
625: {\
626: PetscInt _k,_entry,_location,_lnkdata;\
627: nlnk = 0;\
628: _lnkdata = idx_start;\
629: for (_k=0; _k<nidx; _k++){\
630: _entry = indices[_k];\
631: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
632: /* search for insertion location */\
633: do {\
634: _location = _lnkdata;\
635: _lnkdata = lnk[_location];\
636: } while (_entry > _lnkdata);\
637: /* insertion location is found, add entry into lnk */\
638: lnk[_location] = _entry;\
639: lnk[_entry] = _lnkdata;\
640: nlnk++;\
641: _lnkdata = _entry; /* next search starts from here */\
642: }\
643: }\
644: }
646: /*
647: Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
648: Same as PetscLLAddSorted() with an additional operation:
649: count the number of input indices that are no larger than 'diag'
650: Input Parameters:
651: indices - sorted interger array
652: idx_start - starting index of the list, index of pivot row
653: lnk - linked list(an integer array) that is created
654: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
655: diag - index of the active row in LUFactorSymbolic
656: nzbd - number of input indices with indices <= idx_start
657: im - im[idx_start] is initialized as num of nonzero entries in row=idx_start
658: output Parameters:
659: nlnk - number of newly added indices
660: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
661: bt - updated PetscBT (bitarray)
662: im - im[idx_start]: unchanged if diag is not an entry
663: : num of entries with indices <= diag if diag is an entry
664: */
665: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
666: {\
667: PetscInt _k,_entry,_location,_lnkdata,_nidx;\
668: nlnk = 0;\
669: _lnkdata = idx_start;\
670: _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
671: for (_k=0; _k<_nidx; _k++){\
672: _entry = indices[_k];\
673: nzbd++;\
674: if ( _entry== diag) im[idx_start] = nzbd;\
675: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
676: /* search for insertion location */\
677: do {\
678: _location = _lnkdata;\
679: _lnkdata = lnk[_location];\
680: } while (_entry > _lnkdata);\
681: /* insertion location is found, add entry into lnk */\
682: lnk[_location] = _entry;\
683: lnk[_entry] = _lnkdata;\
684: nlnk++;\
685: _lnkdata = _entry; /* next search starts from here */\
686: }\
687: }\
688: }
690: /*
691: Copy data on the list into an array, then initialize the list
692: Input Parameters:
693: idx_start - starting index of the list
694: lnk_max - max value of lnk indicating the end of the list
695: nlnk - number of data on the list to be copied
696: lnk - linked list
697: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
698: output Parameters:
699: indices - array that contains the copied data
700: lnk - linked list that is cleaned and initialize
701: bt - PetscBT (bitarray) with all bits set to false
702: */
703: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
704: {\
705: PetscInt _j,_idx=idx_start;\
706: for (_j=0; _j<nlnk; _j++){\
707: _idx = lnk[_idx];\
708: *(indices+_j) = _idx;\
709: PetscBTClear(bt,_idx);\
710: }\
711: lnk[idx_start] = lnk_max;\
712: }
713: /*
714: Free memories used by the list
715: */
716: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
718: /* Routines below are used for incomplete matrix factorization */
719: /*
720: Create and initialize a linked list and its levels
721: Input Parameters:
722: idx_start - starting index of the list
723: lnk_max - max value of lnk indicating the end of the list
724: nlnk - max length of the list
725: Output Parameters:
726: lnk - list initialized
727: lnk_lvl - array of size nlnk for storing levels of lnk
728: bt - PetscBT (bitarray) with all bits set to false
729: */
730: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
731: (PetscMalloc(2*nlnk*sizeof(PetscInt),&lnk) || PetscBTCreate(nlnk,bt) || PetscBTMemzero(nlnk,bt) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))
733: /*
734: Initialize a sorted linked list used for ILU and ICC
735: Input Parameters:
736: nidx - number of input idx
737: idx - interger array used for storing column indices
738: idx_start - starting index of the list
739: perm - indices of an IS
740: lnk - linked list(an integer array) that is created
741: lnklvl - levels of lnk
742: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
743: output Parameters:
744: nlnk - number of newly added idx
745: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
746: lnklvl - levels of lnk
747: bt - updated PetscBT (bitarray)
748: */
749: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
750: {\
751: PetscInt _k,_entry,_location,_lnkdata;\
752: nlnk = 0;\
753: _lnkdata = idx_start;\
754: for (_k=0; _k<nidx; _k++){\
755: _entry = perm[idx[_k]];\
756: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
757: /* search for insertion location */\
758: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
759: do {\
760: _location = _lnkdata;\
761: _lnkdata = lnk[_location];\
762: } while (_entry > _lnkdata);\
763: /* insertion location is found, add entry into lnk */\
764: lnk[_location] = _entry;\
765: lnk[_entry] = _lnkdata;\
766: lnklvl[_entry] = 0;\
767: nlnk++;\
768: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
769: }\
770: }\
771: }
773: /*
774: Add a SORTED index set into a sorted linked list for ILU
775: Input Parameters:
776: nidx - number of input indices
777: idx - sorted interger array used for storing column indices
778: level - level of fill, e.g., ICC(level)
779: idxlvl - level of idx
780: idx_start - starting index of the list
781: lnk - linked list(an integer array) that is created
782: lnklvl - levels of lnk
783: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
784: prow - the row number of idx
785: output Parameters:
786: nlnk - number of newly added idx
787: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
788: lnklvl - levels of lnk
789: bt - updated PetscBT (bitarray)
791: Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
792: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
793: */
794: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
795: {\
796: PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
797: nlnk = 0;\
798: _lnkdata = idx_start;\
799: for (_k=0; _k<nidx; _k++){\
800: _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
801: if (_incrlev > level) continue;\
802: _entry = idx[_k];\
803: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
804: /* search for insertion location */\
805: do {\
806: _location = _lnkdata;\
807: _lnkdata = lnk[_location];\
808: } while (_entry > _lnkdata);\
809: /* insertion location is found, add entry into lnk */\
810: lnk[_location] = _entry;\
811: lnk[_entry] = _lnkdata;\
812: lnklvl[_entry] = _incrlev;\
813: nlnk++;\
814: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
815: } else { /* existing entry: update lnklvl */\
816: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
817: }\
818: }\
819: }
821: /*
822: Add a index set into a sorted linked list
823: Input Parameters:
824: nidx - number of input idx
825: idx - interger array used for storing column indices
826: level - level of fill, e.g., ICC(level)
827: idxlvl - level of idx
828: idx_start - starting index of the list
829: lnk - linked list(an integer array) that is created
830: lnklvl - levels of lnk
831: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
832: output Parameters:
833: nlnk - number of newly added idx
834: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
835: lnklvl - levels of lnk
836: bt - updated PetscBT (bitarray)
837: */
838: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
839: {\
840: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
841: nlnk = 0;\
842: _lnkdata = idx_start;\
843: for (_k=0; _k<nidx; _k++){\
844: _incrlev = idxlvl[_k] + 1;\
845: if (_incrlev > level) continue;\
846: _entry = idx[_k];\
847: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
848: /* search for insertion location */\
849: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
850: do {\
851: _location = _lnkdata;\
852: _lnkdata = lnk[_location];\
853: } while (_entry > _lnkdata);\
854: /* insertion location is found, add entry into lnk */\
855: lnk[_location] = _entry;\
856: lnk[_entry] = _lnkdata;\
857: lnklvl[_entry] = _incrlev;\
858: nlnk++;\
859: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
860: } else { /* existing entry: update lnklvl */\
861: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
862: }\
863: }\
864: }
866: /*
867: Add a SORTED index set into a sorted linked list
868: Input Parameters:
869: nidx - number of input indices
870: idx - sorted interger array used for storing column indices
871: level - level of fill, e.g., ICC(level)
872: idxlvl - level of idx
873: idx_start - starting index of the list
874: lnk - linked list(an integer array) that is created
875: lnklvl - levels of lnk
876: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
877: output Parameters:
878: nlnk - number of newly added idx
879: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
880: lnklvl - levels of lnk
881: bt - updated PetscBT (bitarray)
882: */
883: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
884: {\
885: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
886: nlnk = 0;\
887: _lnkdata = idx_start;\
888: for (_k=0; _k<nidx; _k++){\
889: _incrlev = idxlvl[_k] + 1;\
890: if (_incrlev > level) continue;\
891: _entry = idx[_k];\
892: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
893: /* search for insertion location */\
894: do {\
895: _location = _lnkdata;\
896: _lnkdata = lnk[_location];\
897: } while (_entry > _lnkdata);\
898: /* insertion location is found, add entry into lnk */\
899: lnk[_location] = _entry;\
900: lnk[_entry] = _lnkdata;\
901: lnklvl[_entry] = _incrlev;\
902: nlnk++;\
903: _lnkdata = _entry; /* next search starts from here */\
904: } else { /* existing entry: update lnklvl */\
905: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
906: }\
907: }\
908: }
910: /*
911: Add a SORTED index set into a sorted linked list for ICC
912: Input Parameters:
913: nidx - number of input indices
914: idx - sorted interger array used for storing column indices
915: level - level of fill, e.g., ICC(level)
916: idxlvl - level of idx
917: idx_start - starting index of the list
918: lnk - linked list(an integer array) that is created
919: lnklvl - levels of lnk
920: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
921: idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
922: output Parameters:
923: nlnk - number of newly added indices
924: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
925: lnklvl - levels of lnk
926: bt - updated PetscBT (bitarray)
927: Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
928: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
929: */
930: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
931: {\
932: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
933: nlnk = 0;\
934: _lnkdata = idx_start;\
935: for (_k=0; _k<nidx; _k++){\
936: _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
937: if (_incrlev > level) continue;\
938: _entry = idx[_k];\
939: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
940: /* search for insertion location */\
941: do {\
942: _location = _lnkdata;\
943: _lnkdata = lnk[_location];\
944: } while (_entry > _lnkdata);\
945: /* insertion location is found, add entry into lnk */\
946: lnk[_location] = _entry;\
947: lnk[_entry] = _lnkdata;\
948: lnklvl[_entry] = _incrlev;\
949: nlnk++;\
950: _lnkdata = _entry; /* next search starts from here */\
951: } else { /* existing entry: update lnklvl */\
952: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
953: }\
954: }\
955: }
957: /*
958: Copy data on the list into an array, then initialize the list
959: Input Parameters:
960: idx_start - starting index of the list
961: lnk_max - max value of lnk indicating the end of the list
962: nlnk - number of data on the list to be copied
963: lnk - linked list
964: lnklvl - level of lnk
965: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
966: output Parameters:
967: indices - array that contains the copied data
968: lnk - linked list that is cleaned and initialize
969: lnklvl - level of lnk that is reinitialized
970: bt - PetscBT (bitarray) with all bits set to false
971: */
972: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
973: {\
974: PetscInt _j,_idx=idx_start;\
975: for (_j=0; _j<nlnk; _j++){\
976: _idx = lnk[_idx];\
977: *(indices+_j) = _idx;\
978: *(indiceslvl+_j) = lnklvl[_idx];\
979: lnklvl[_idx] = -1;\
980: PetscBTClear(bt,_idx);\
981: }\
982: lnk[idx_start] = lnk_max;\
983: }
984: /*
985: Free memories used by the list
986: */
987: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(bt))
1008: #endif