Actual source code: matimpl.h
petsc-3.6.4 2016-04-12
2: #ifndef __MATIMPL_H
5: #include <petscmat.h>
6: #include <petsc/private/petscimpl.h>
8: PETSC_EXTERN PetscBool MatRegisterAllCalled;
9: PETSC_EXTERN PetscBool MatOrderingRegisterAllCalled;
10: PETSC_EXTERN PetscBool MatColoringRegisterAllCalled;
11: PETSC_EXTERN PetscBool MatPartitioningRegisterAllCalled;
12: PETSC_EXTERN PetscBool MatCoarsenRegisterAllCalled;
13: PETSC_EXTERN PetscErrorCode MatRegisterAll(void);
14: PETSC_EXTERN PetscErrorCode MatOrderingRegisterAll(void);
15: PETSC_EXTERN PetscErrorCode MatColoringRegisterAll(void);
16: PETSC_EXTERN PetscErrorCode MatPartitioningRegisterAll(void);
17: PETSC_EXTERN PetscErrorCode MatCoarsenRegisterAll(void);
19: /*
20: This file defines the parts of the matrix data structure that are
21: shared by all matrix types.
22: */
24: /*
25: If you add entries here also add them to the MATOP enum
26: in include/petscmat.h and include/petsc/finclude/petscmat.h
27: */
28: typedef struct _MatOps *MatOps;
29: struct _MatOps {
30: /* 0*/
31: PetscErrorCode (*setvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
32: PetscErrorCode (*getrow)(Mat,PetscInt,PetscInt *,PetscInt*[],PetscScalar*[]);
33: PetscErrorCode (*restorerow)(Mat,PetscInt,PetscInt *,PetscInt *[],PetscScalar *[]);
34: PetscErrorCode (*mult)(Mat,Vec,Vec);
35: PetscErrorCode (*multadd)(Mat,Vec,Vec,Vec);
36: /* 5*/
37: PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
38: PetscErrorCode (*multtransposeadd)(Mat,Vec,Vec,Vec);
39: PetscErrorCode (*solve)(Mat,Vec,Vec);
40: PetscErrorCode (*solveadd)(Mat,Vec,Vec,Vec);
41: PetscErrorCode (*solvetranspose)(Mat,Vec,Vec);
42: /*10*/
43: PetscErrorCode (*solvetransposeadd)(Mat,Vec,Vec,Vec);
44: PetscErrorCode (*lufactor)(Mat,IS,IS,const MatFactorInfo*);
45: PetscErrorCode (*choleskyfactor)(Mat,IS,const MatFactorInfo*);
46: PetscErrorCode (*sor)(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
47: PetscErrorCode (*transpose)(Mat,MatReuse,Mat *);
48: /*15*/
49: PetscErrorCode (*getinfo)(Mat,MatInfoType,MatInfo*);
50: PetscErrorCode (*equal)(Mat,Mat,PetscBool *);
51: PetscErrorCode (*getdiagonal)(Mat,Vec);
52: PetscErrorCode (*diagonalscale)(Mat,Vec,Vec);
53: PetscErrorCode (*norm)(Mat,NormType,PetscReal*);
54: /*20*/
55: PetscErrorCode (*assemblybegin)(Mat,MatAssemblyType);
56: PetscErrorCode (*assemblyend)(Mat,MatAssemblyType);
57: PetscErrorCode (*setoption)(Mat,MatOption,PetscBool );
58: PetscErrorCode (*zeroentries)(Mat);
59: /*24*/
60: PetscErrorCode (*zerorows)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
61: PetscErrorCode (*lufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
62: PetscErrorCode (*lufactornumeric)(Mat,Mat,const MatFactorInfo*);
63: PetscErrorCode (*choleskyfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
64: PetscErrorCode (*choleskyfactornumeric)(Mat,Mat,const MatFactorInfo*);
65: /*29*/
66: PetscErrorCode (*setup)(Mat);
67: PetscErrorCode (*ilufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
68: PetscErrorCode (*iccfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
69: PetscErrorCode (*placeholder_32)(Mat);
70: PetscErrorCode (*placeholder_33)(Mat);
71: /*34*/
72: PetscErrorCode (*duplicate)(Mat,MatDuplicateOption,Mat*);
73: PetscErrorCode (*forwardsolve)(Mat,Vec,Vec);
74: PetscErrorCode (*backwardsolve)(Mat,Vec,Vec);
75: PetscErrorCode (*ilufactor)(Mat,IS,IS,const MatFactorInfo*);
76: PetscErrorCode (*iccfactor)(Mat,IS,const MatFactorInfo*);
77: /*39*/
78: PetscErrorCode (*axpy)(Mat,PetscScalar,Mat,MatStructure);
79: PetscErrorCode (*getsubmatrices)(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *[]);
80: PetscErrorCode (*increaseoverlap)(Mat,PetscInt,IS[],PetscInt);
81: PetscErrorCode (*getvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscScalar []);
82: PetscErrorCode (*copy)(Mat,Mat,MatStructure);
83: /*44*/
84: PetscErrorCode (*getrowmax)(Mat,Vec,PetscInt[]);
85: PetscErrorCode (*scale)(Mat,PetscScalar);
86: PetscErrorCode (*shift)(Mat,PetscScalar);
87: PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
88: PetscErrorCode (*zerorowscolumns)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
89: /*49*/
90: PetscErrorCode (*setrandom)(Mat,PetscRandom);
91: PetscErrorCode (*getrowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
92: PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt *,const PetscInt *[],const PetscInt *[],PetscBool *);
93: PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
94: PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
95: /*54*/
96: PetscErrorCode (*fdcoloringcreate)(Mat,ISColoring,MatFDColoring);
97: PetscErrorCode (*coloringpatch)(Mat,PetscInt,PetscInt,ISColoringValue[],ISColoring*);
98: PetscErrorCode (*setunfactored)(Mat);
99: PetscErrorCode (*permute)(Mat,IS,IS,Mat*);
100: PetscErrorCode (*setvaluesblocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
101: /*59*/
102: PetscErrorCode (*getsubmatrix)(Mat,IS,IS,MatReuse,Mat*);
103: PetscErrorCode (*destroy)(Mat);
104: PetscErrorCode (*view)(Mat,PetscViewer);
105: PetscErrorCode (*convertfrom)(Mat, MatType,MatReuse,Mat*);
106: PetscErrorCode (*matmatmult)(Mat,Mat,Mat,MatReuse,PetscReal,Mat*);
107: /*64*/
108: PetscErrorCode (*matmatmultsymbolic)(Mat,Mat,Mat,PetscReal,Mat*);
109: PetscErrorCode (*matmatmultnumeric)(Mat,Mat,Mat,Mat);
110: PetscErrorCode (*setlocaltoglobalmapping)(Mat,ISLocalToGlobalMapping,ISLocalToGlobalMapping);
111: PetscErrorCode (*setvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
112: PetscErrorCode (*zerorowslocal)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
113: /*69*/
114: PetscErrorCode (*getrowmaxabs)(Mat,Vec,PetscInt[]);
115: PetscErrorCode (*getrowminabs)(Mat,Vec,PetscInt[]);
116: PetscErrorCode (*convert)(Mat, MatType,MatReuse,Mat*);
117: PetscErrorCode (*setcoloring)(Mat,ISColoring);
118: PetscErrorCode (*placeholder_73)(Mat,void*);
119: /*74*/
120: PetscErrorCode (*setvaluesadifor)(Mat,PetscInt,void*);
121: PetscErrorCode (*fdcoloringapply)(Mat,MatFDColoring,Vec,void*);
122: PetscErrorCode (*setfromoptions)(PetscOptions*,Mat);
123: PetscErrorCode (*multconstrained)(Mat,Vec,Vec);
124: PetscErrorCode (*multtransposeconstrained)(Mat,Vec,Vec);
125: /*79*/
126: PetscErrorCode (*findzerodiagonals)(Mat,IS*);
127: PetscErrorCode (*mults)(Mat, Vecs, Vecs);
128: PetscErrorCode (*solves)(Mat, Vecs, Vecs);
129: PetscErrorCode (*getinertia)(Mat,PetscInt*,PetscInt*,PetscInt*);
130: PetscErrorCode (*load)(Mat, PetscViewer);
131: /*84*/
132: PetscErrorCode (*issymmetric)(Mat,PetscReal,PetscBool *);
133: PetscErrorCode (*ishermitian)(Mat,PetscReal,PetscBool *);
134: PetscErrorCode (*isstructurallysymmetric)(Mat,PetscBool *);
135: PetscErrorCode (*setvaluesblockedlocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
136: PetscErrorCode (*getvecs)(Mat,Vec*,Vec*);
137: /*89*/
138: PetscErrorCode (*matmult)(Mat,Mat,MatReuse,PetscReal,Mat*);
139: PetscErrorCode (*matmultsymbolic)(Mat,Mat,PetscReal,Mat*);
140: PetscErrorCode (*matmultnumeric)(Mat,Mat,Mat);
141: PetscErrorCode (*ptap)(Mat,Mat,MatReuse,PetscReal,Mat*);
142: PetscErrorCode (*ptapsymbolic)(Mat,Mat,PetscReal,Mat*); /* double dispatch wrapper routine */
143: /*94*/
144: PetscErrorCode (*ptapnumeric)(Mat,Mat,Mat); /* double dispatch wrapper routine */
145: PetscErrorCode (*mattransposemult)(Mat,Mat,MatReuse,PetscReal,Mat*);
146: PetscErrorCode (*mattransposemultsymbolic)(Mat,Mat,PetscReal,Mat*);
147: PetscErrorCode (*mattransposemultnumeric)(Mat,Mat,Mat);
148: PetscErrorCode (*placeholder_98)(Mat);
149: /*99*/
150: PetscErrorCode (*placeholder_99)(Mat);
151: PetscErrorCode (*placeholder_100)(Mat);
152: PetscErrorCode (*placeholder_101)(Mat);
153: PetscErrorCode (*conjugate)(Mat); /* complex conjugate */
154: PetscErrorCode (*placeholder_103)(void);
155: /*104*/
156: PetscErrorCode (*setvaluesrow)(Mat,PetscInt,const PetscScalar[]);
157: PetscErrorCode (*realpart)(Mat);
158: PetscErrorCode (*imaginarypart)(Mat);
159: PetscErrorCode (*getrowuppertriangular)(Mat);
160: PetscErrorCode (*restorerowuppertriangular)(Mat);
161: /*109*/
162: PetscErrorCode (*matsolve)(Mat,Mat,Mat);
163: PetscErrorCode (*placeholder_110)(Mat);
164: PetscErrorCode (*getrowmin)(Mat,Vec,PetscInt[]);
165: PetscErrorCode (*getcolumnvector)(Mat,Vec,PetscInt);
166: PetscErrorCode (*missingdiagonal)(Mat,PetscBool *,PetscInt*);
167: /*114*/
168: PetscErrorCode (*getseqnonzerostructure)(Mat,Mat *);
169: PetscErrorCode (*create)(Mat);
170: PetscErrorCode (*getghosts)(Mat,PetscInt*,const PetscInt *[]);
171: PetscErrorCode (*getlocalsubmatrix)(Mat,IS,IS,Mat*);
172: PetscErrorCode (*restorelocalsubmatrix)(Mat,IS,IS,Mat*);
173: /*119*/
174: PetscErrorCode (*multdiagonalblock)(Mat,Vec,Vec);
175: PetscErrorCode (*hermitiantranspose)(Mat,MatReuse,Mat*);
176: PetscErrorCode (*multhermitiantranspose)(Mat,Vec,Vec);
177: PetscErrorCode (*multhermitiantransposeadd)(Mat,Vec,Vec,Vec);
178: PetscErrorCode (*getmultiprocblock)(Mat,MPI_Comm,MatReuse,Mat*);
179: /*124*/
180: PetscErrorCode (*findnonzerorows)(Mat,IS*);
181: PetscErrorCode (*getcolumnnorms)(Mat,NormType,PetscReal*);
182: PetscErrorCode (*invertblockdiagonal)(Mat,const PetscScalar**);
183: PetscErrorCode (*placeholder_127)(Mat,Vec,Vec,Vec);
184: PetscErrorCode (*getsubmatricesmpi)(Mat,PetscInt,const IS[], const IS[], MatReuse, Mat**);
185: /*129*/
186: PetscErrorCode (*setvaluesbatch)(Mat,PetscInt,PetscInt,PetscInt*,const PetscScalar*);
187: PetscErrorCode (*transposematmult)(Mat,Mat,MatReuse,PetscReal,Mat*);
188: PetscErrorCode (*transposematmultsymbolic)(Mat,Mat,PetscReal,Mat*);
189: PetscErrorCode (*transposematmultnumeric)(Mat,Mat,Mat);
190: PetscErrorCode (*transposecoloringcreate)(Mat,ISColoring,MatTransposeColoring);
191: /*134*/
192: PetscErrorCode (*transcoloringapplysptoden)(MatTransposeColoring,Mat,Mat);
193: PetscErrorCode (*transcoloringapplydentosp)(MatTransposeColoring,Mat,Mat);
194: PetscErrorCode (*rart)(Mat,Mat,MatReuse,PetscReal,Mat*);
195: PetscErrorCode (*rartsymbolic)(Mat,Mat,PetscReal,Mat*); /* double dispatch wrapper routine */
196: PetscErrorCode (*rartnumeric)(Mat,Mat,Mat); /* double dispatch wrapper routine */
197: /*139*/
198: PetscErrorCode (*setblocksizes)(Mat,PetscInt,PetscInt);
199: PetscErrorCode (*aypx)(Mat,PetscScalar,Mat,MatStructure);
200: PetscErrorCode (*residual)(Mat,Vec,Vec,Vec);
201: PetscErrorCode (*fdcoloringsetup)(Mat,ISColoring,MatFDColoring);
202: PetscErrorCode (*findoffblockdiagonalentries)(Mat,IS*);
203: /*144*/
204: PetscErrorCode (*creatempimatconcatenateseqmat)(MPI_Comm,Mat,PetscInt,MatReuse,Mat*);
206: };
207: /*
208: If you add MatOps entries above also add them to the MATOP enum
209: in include/petscmat.h and include/petsc/finclude/petscmat.h
210: */
212: #include <petscsys.h>
213: PETSC_EXTERN PetscErrorCode MatRegisterOp(MPI_Comm, const char[], PetscVoidFunction, const char[], PetscInt, ...);
214: PETSC_EXTERN PetscErrorCode MatQueryOp(MPI_Comm, PetscVoidFunction*, const char[], PetscInt, ...);
216: typedef struct _p_MatBaseName* MatBaseName;
217: struct _p_MatBaseName {
218: char *bname,*sname,*mname;
219: MatBaseName next;
220: };
222: PETSC_EXTERN MatBaseName MatBaseNameList;
224: /*
225: Utility private matrix routines
226: */
227: PETSC_INTERN PetscErrorCode MatConvert_Basic(Mat, MatType,MatReuse,Mat*);
228: PETSC_INTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
229: PETSC_INTERN PetscErrorCode MatHeaderMerge(Mat,Mat);
230: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat,Mat);
231: PETSC_INTERN PetscErrorCode MatDiagonalSet_Default(Mat,Vec,InsertMode);
233: #if defined(PETSC_USE_DEBUG)
234: # define MatCheckPreallocated(A,arg) do { \
235: if (PetscUnlikely(!(A)->preallocated)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatXXXSetPreallocation() or MatSetUp() on argument %D \"%s\" before %s()",(arg),#A,PETSC_FUNCTION_NAME); \
236: } while (0)
237: #else
238: # define MatCheckPreallocated(A,arg) do {} while (0)
239: #endif
241: /*
242: The stash is used to temporarily store inserted matrix values that
243: belong to another processor. During the assembly phase the stashed
244: values are moved to the correct processor and
245: */
247: typedef struct _MatStashSpace *PetscMatStashSpace;
249: struct _MatStashSpace {
250: PetscMatStashSpace next;
251: PetscScalar *space_head,*val;
252: PetscInt *idx,*idy;
253: PetscInt total_space_size;
254: PetscInt local_used;
255: PetscInt local_remaining;
256: };
258: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceGet(PetscInt,PetscInt,PetscMatStashSpace *);
259: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceContiguous(PetscInt,PetscMatStashSpace *,PetscScalar *,PetscInt *,PetscInt *);
260: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceDestroy(PetscMatStashSpace*);
262: typedef struct {
263: PetscInt nmax; /* maximum stash size */
264: PetscInt umax; /* user specified max-size */
265: PetscInt oldnmax; /* the nmax value used previously */
266: PetscInt n; /* stash size */
267: PetscInt bs; /* block size of the stash */
268: PetscInt reallocs; /* preserve the no of mallocs invoked */
269: PetscMatStashSpace space_head,space; /* linked list to hold stashed global row/column numbers and matrix values */
270: /* The following variables are used for communication */
271: MPI_Comm comm;
272: PetscMPIInt size,rank;
273: PetscMPIInt tag1,tag2;
274: MPI_Request *send_waits; /* array of send requests */
275: MPI_Request *recv_waits; /* array of receive requests */
276: MPI_Status *send_status; /* array of send status */
277: PetscInt nsends,nrecvs; /* numbers of sends and receives */
278: PetscScalar *svalues; /* sending data */
279: PetscInt *sindices;
280: PetscScalar **rvalues; /* receiving data (values) */
281: PetscInt **rindices; /* receiving data (indices) */
282: PetscInt nprocessed; /* number of messages already processed */
283: PetscMPIInt *flg_v; /* indicates what messages have arrived so far and from whom */
284: PetscBool reproduce;
285: PetscInt reproduce_count;
286: } MatStash;
288: PETSC_INTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
289: PETSC_INTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
290: PETSC_INTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
291: PETSC_INTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
292: PETSC_INTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
293: PETSC_INTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscBool );
294: PETSC_INTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscBool );
295: PETSC_INTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
296: PETSC_INTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
297: PETSC_INTERN PetscErrorCode MatStashScatterBegin_Private(Mat,MatStash*,PetscInt*);
298: PETSC_INTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,PetscScalar**,PetscInt*);
300: typedef struct {
301: PetscInt dim;
302: PetscInt dims[4];
303: PetscInt starts[4];
304: PetscBool noc; /* this is a single component problem, hence user will not set MatStencil.c */
305: } MatStencilInfo;
307: /* Info about using compressed row format */
308: typedef struct {
309: PetscBool use; /* indicates compressed rows have been checked and will be used */
310: PetscInt nrows; /* number of non-zero rows */
311: PetscInt *i; /* compressed row pointer */
312: PetscInt *rindex; /* compressed row index */
313: } Mat_CompressedRow;
314: PETSC_EXTERN PetscErrorCode MatCheckCompressedRow(Mat,PetscInt,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);
316: typedef struct { /* used by MatCreateRedundantMatrix() for reusing matredundant */
317: PetscInt nzlocal,nsends,nrecvs;
318: PetscMPIInt *send_rank,*recv_rank;
319: PetscInt *sbuf_nz,*rbuf_nz,*sbuf_j,**rbuf_j;
320: PetscScalar *sbuf_a,**rbuf_a;
321: MPI_Comm subcomm; /* when user does not provide a subcomm */
322: IS isrow,iscol;
323: Mat *matseq;
324: } Mat_Redundant;
326: struct _p_Mat {
327: PETSCHEADER(struct _MatOps);
328: PetscLayout rmap,cmap;
329: void *data; /* implementation-specific data */
330: MatFactorType factortype; /* MAT_FACTOR_LU, ILU, CHOLESKY or ICC */
331: PetscBool assembled; /* is the matrix assembled? */
332: PetscBool was_assembled; /* new values inserted into assembled mat */
333: PetscInt num_ass; /* number of times matrix has been assembled */
334: PetscObjectState nonzerostate; /* each time new nonzeros locations are introduced into the matrix this is updated */
335: MatInfo info; /* matrix information */
336: InsertMode insertmode; /* have values been inserted in matrix or added? */
337: MatStash stash,bstash; /* used for assembling off-proc mat emements */
338: MatNullSpace nullsp; /* null space (operator is singular) */
339: MatNullSpace transnullsp; /* null space of transpose of operator */
340: MatNullSpace nearnullsp; /* near null space to be used by multigrid methods */
341: PetscBool preallocated;
342: MatStencilInfo stencil; /* information for structured grid */
343: PetscBool symmetric,hermitian,structurally_symmetric,spd;
344: PetscBool symmetric_set,hermitian_set,structurally_symmetric_set,spd_set; /* if true, then corresponding flag is correct*/
345: PetscBool symmetric_eternal;
346: PetscBool nooffprocentries,nooffproczerorows;
347: #if defined(PETSC_HAVE_CUSP)
348: PetscCUSPFlag valid_GPU_matrix; /* flag pointing to the matrix on the gpu*/
349: #endif
350: #if defined(PETSC_HAVE_VIENNACL)
351: PetscViennaCLFlag valid_GPU_matrix; /* flag pointing to the matrix on the gpu*/
352: #endif
353: void *spptr; /* pointer for special library like SuperLU */
354: MatSolverPackage solvertype;
355: PetscBool checksymmetryonassembly,checknullspaceonassembly;
356: PetscReal checksymmetrytol;
357: Mat_Redundant *redundant; /* used by MatCreateRedundantMatrix() */
358: PetscBool erroriffpe; /* Generate an error if FPE detected (for example a zero pivot) instead of returning*/
359: };
361: PETSC_INTERN PetscErrorCode MatAXPY_Basic(Mat,PetscScalar,Mat,MatStructure);
362: PETSC_INTERN PetscErrorCode MatAXPY_BasicWithPreallocation(Mat,Mat,PetscScalar,Mat,MatStructure);
363: /*
364: Object for partitioning graphs
365: */
367: typedef struct _MatPartitioningOps *MatPartitioningOps;
368: struct _MatPartitioningOps {
369: PetscErrorCode (*apply)(MatPartitioning,IS*);
370: PetscErrorCode (*setfromoptions)(PetscOptions*,MatPartitioning);
371: PetscErrorCode (*destroy)(MatPartitioning);
372: PetscErrorCode (*view)(MatPartitioning,PetscViewer);
373: };
375: struct _p_MatPartitioning {
376: PETSCHEADER(struct _MatPartitioningOps);
377: Mat adj;
378: PetscInt *vertex_weights;
379: PetscReal *part_weights;
380: PetscInt n; /* number of partitions */
381: void *data;
382: PetscInt setupcalled;
383: };
385: /*
386: Object for coarsen graphs
387: */
388: typedef struct _MatCoarsenOps *MatCoarsenOps;
389: struct _MatCoarsenOps {
390: PetscErrorCode (*apply)(MatCoarsen);
391: PetscErrorCode (*setfromoptions)(PetscOptions*,MatCoarsen);
392: PetscErrorCode (*destroy)(MatCoarsen);
393: PetscErrorCode (*view)(MatCoarsen,PetscViewer);
394: };
396: struct _p_MatCoarsen {
397: PETSCHEADER(struct _MatCoarsenOps);
398: Mat graph;
399: PetscInt setupcalled;
400: void *subctx;
401: /* */
402: PetscBool strict_aggs;
403: IS perm;
404: PetscCoarsenData *agg_lists;
405: };
407: PETSC_EXTERN PetscErrorCode PetscCDCreate(PetscInt,PetscCoarsenData**);
408: PETSC_EXTERN PetscErrorCode PetscCDDestroy(PetscCoarsenData*);
409: PETSC_EXTERN PetscErrorCode PetscLLNSetID(PetscCDIntNd*,PetscInt);
410: PETSC_EXTERN PetscErrorCode PetscLLNGetID(const PetscCDIntNd*,PetscInt*);
411: PETSC_EXTERN PetscErrorCode PetscCDAppendID(PetscCoarsenData*,PetscInt,PetscInt);
412: PETSC_EXTERN PetscErrorCode PetscCDAppendRemove(PetscCoarsenData*,PetscInt,PetscInt);
413: PETSC_EXTERN PetscErrorCode PetscCDAppendNode(PetscCoarsenData*,PetscInt,PetscCDIntNd*);
414: PETSC_EXTERN PetscErrorCode PetscCDRemoveNextNode(PetscCoarsenData*,PetscInt,PetscCDIntNd*);
415: PETSC_EXTERN PetscErrorCode PetscCDRemoveAllAt(PetscCoarsenData*,PetscInt);
416: PETSC_EXTERN PetscErrorCode PetscCDSizeAt(const PetscCoarsenData*,PetscInt,PetscInt*);
417: PETSC_EXTERN PetscErrorCode PetscCDEmptyAt(const PetscCoarsenData*,PetscInt,PetscBool*);
418: PETSC_EXTERN PetscErrorCode PetscCDSetChuckSize(PetscCoarsenData*,PetscInt);
419: PETSC_EXTERN PetscErrorCode PetscCDPrint(const PetscCoarsenData*,MPI_Comm);
420: PETSC_EXTERN PetscErrorCode PetscCDGetMIS(PetscCoarsenData*,IS*);
421: PETSC_EXTERN PetscErrorCode PetscCDGetMat(const PetscCoarsenData*,Mat*);
422: PETSC_EXTERN PetscErrorCode PetscCDSetMat(PetscCoarsenData*,Mat);
424: typedef PetscCDIntNd *PetscCDPos;
425: PETSC_EXTERN PetscErrorCode PetscCDGetHeadPos(const PetscCoarsenData*,PetscInt,PetscCDPos*);
426: PETSC_EXTERN PetscErrorCode PetscCDGetNextPos(const PetscCoarsenData*,PetscInt,PetscCDPos*);
427: PETSC_EXTERN PetscErrorCode PetscCDGetASMBlocks(const PetscCoarsenData*,const PetscInt,PetscInt*,IS**);
428: /* PetscErrorCode PetscCDSetRemovedIS( PetscCoarsenData *ail, MPI_Comm, const PetscInt, PetscInt[] ); */
429: /* PetscErrorCode PetscCDGetRemovedIS( PetscCoarsenData *ail, IS * ); */
431: /*
432: MatFDColoring is used to compute Jacobian matrices efficiently
433: via coloring. The data structure is explained below in an example.
435: Color = 0 1 0 2 | 2 3 0
436: ---------------------------------------------------
437: 00 01 | 05
438: 10 11 | 14 15 Processor 0
439: 22 23 | 25
440: 32 33 |
441: ===================================================
442: | 44 45 46
443: 50 | 55 Processor 1
444: | 64 66
445: ---------------------------------------------------
447: ncolors = 4;
449: ncolumns = {2,1,1,0}
450: columns = {{0,2},{1},{3},{}}
451: nrows = {4,2,3,3}
452: rows = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
453: vwscale = {dx(0),dx(1),dx(2),dx(3)} MPI Vec
454: vscale = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)} Seq Vec
456: ncolumns = {1,0,1,1}
457: columns = {{6},{},{4},{5}}
458: nrows = {3,0,2,2}
459: rows = {{0,1,2},{},{1,2},{1,2}}
460: vwscale = {dx(4),dx(5),dx(6)} MPI Vec
461: vscale = {dx(0),dx(4),dx(5),dx(6)} Seq Vec
463: See the routine MatFDColoringApply() for how this data is used
464: to compute the Jacobian.
466: */
467: typedef struct {
468: PetscInt row;
469: PetscInt col;
470: PetscScalar *valaddr; /* address of value */
471: } MatEntry;
473: typedef struct {
474: PetscInt row;
475: PetscScalar *valaddr; /* address of value */
476: } MatEntry2;
478: struct _p_MatFDColoring{
479: PETSCHEADER(int);
480: PetscInt M,N,m; /* total rows, columns; local rows */
481: PetscInt rstart; /* first row owned by local processor */
482: PetscInt ncolors; /* number of colors */
483: PetscInt *ncolumns; /* number of local columns for a color */
484: PetscInt **columns; /* lists the local columns of each color (using global column numbering) */
485: PetscInt *nrows; /* number of local rows for each color */
486: MatEntry *matentry; /* holds (row, column, address of value) for Jacobian matrix entry */
487: MatEntry2 *matentry2; /* holds (row, address of value) for Jacobian matrix entry */
488: PetscScalar *dy; /* store a block of F(x+dx)-F(x) when J is in BAIJ format */
489: PetscReal error_rel; /* square root of relative error in computing function */
490: PetscReal umin; /* minimum allowable u'dx value */
491: Vec w1,w2,w3; /* work vectors used in computing Jacobian */
492: PetscBool fset; /* indicates that the initial function value F(X) is set */
493: PetscErrorCode (*f)(void); /* function that defines Jacobian */
494: void *fctx; /* optional user-defined context for use by the function f */
495: Vec vscale; /* holds FD scaling, i.e. 1/dx for each perturbed column */
496: PetscInt currentcolor; /* color for which function evaluation is being done now */
497: const char *htype; /* "wp" or "ds" */
498: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
499: PetscInt brows,bcols; /* number of block rows or columns for speedup inserting the dense matrix into sparse Jacobian */
500: PetscBool setupcalled; /* true if setup has been called */
501: void (*ftn_func_pointer)(void),*ftn_func_cntx; /* serve the same purpose as *fortran_func_pointers in PETSc objects */
502: };
504: typedef struct _MatColoringOps *MatColoringOps;
505: struct _MatColoringOps {
506: PetscErrorCode (*destroy)(MatColoring);
507: PetscErrorCode (*setfromoptions)(PetscOptions*,MatColoring);
508: PetscErrorCode (*view)(MatColoring,PetscViewer);
509: PetscErrorCode (*apply)(MatColoring,ISColoring*);
510: PetscErrorCode (*weights)(MatColoring,PetscReal**,PetscInt**);
511: };
513: struct _p_MatColoring {
514: PETSCHEADER(struct _MatColoringOps);
515: Mat mat;
516: PetscInt dist; /* distance of the coloring */
517: PetscInt maxcolors; /* the maximum number of colors returned, maxcolors=1 for MIS */
518: void *data; /* inner context */
519: PetscBool valid; /* check to see if what is produced is a valid coloring */
520: MatColoringWeightType weight_type; /* type of weight computation to be performed */
521: PetscReal *user_weights; /* custom weights and permutation */
522: PetscInt *user_lperm;
523: };
525: struct _p_MatTransposeColoring{
526: PETSCHEADER(int);
527: PetscInt M,N,m; /* total rows, columns; local rows */
528: PetscInt rstart; /* first row owned by local processor */
529: PetscInt ncolors; /* number of colors */
530: PetscInt *ncolumns; /* number of local columns for a color */
531: PetscInt *nrows; /* number of local rows for each color */
532: PetscInt currentcolor; /* color for which function evaluation is being done now */
533: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_GHOSTED */
535: PetscInt *colorforrow,*colorforcol; /* pointer to rows and columns */
536: PetscInt *rows; /* lists the local rows for each color (using the local row numbering) */
537: PetscInt *den2sp; /* maps (row,color) in the dense matrix to index of sparse matrix array a->a */
538: PetscInt *columns; /* lists the local columns of each color (using global column numbering) */
539: PetscInt brows; /* number of rows for efficient implementation of MatTransColoringApplyDenToSp() */
540: PetscInt *lstart; /* array used for loop over row blocks of Csparse */
541: };
543: /*
544: Null space context for preconditioner/operators
545: */
546: struct _p_MatNullSpace {
547: PETSCHEADER(int);
548: PetscBool has_cnst;
549: PetscInt n;
550: Vec* vecs;
551: PetscScalar* alpha; /* for projections */
552: PetscErrorCode (*remove)(MatNullSpace,Vec,void*); /* for user provided removal function */
553: void* rmctx; /* context for remove() function */
554: };
556: /*
557: Checking zero pivot for LU, ILU preconditioners.
558: */
559: typedef struct {
560: PetscInt nshift,nshift_max;
561: PetscReal shift_amount,shift_lo,shift_hi,shift_top,shift_fraction;
562: PetscBool newshift;
563: PetscReal rs; /* active row sum of abs(offdiagonals) */
564: PetscScalar pv; /* pivot of the active row */
565: } FactorShiftCtx;
567: PETSC_EXTERN PetscErrorCode MatFactorDumpMatrix(Mat);
568: PETSC_EXTERN PetscErrorCode MatShift_Basic(Mat,PetscScalar);
572: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_nz(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
573: {
574: PetscReal _rs = sctx->rs;
575: PetscReal _zero = info->zeropivot*_rs;
578: if (PetscAbsScalar(sctx->pv) <= _zero){
579: /* force |diag| > zeropivot*rs */
580: if (!sctx->nshift) sctx->shift_amount = info->shiftamount;
581: else sctx->shift_amount *= 2.0;
582: sctx->newshift = PETSC_TRUE;
583: (sctx->nshift)++;
584: } else {
585: sctx->newshift = PETSC_FALSE;
586: }
587: return(0);
588: }
592: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_pd(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
593: {
594: PetscReal _rs = sctx->rs;
595: PetscReal _zero = info->zeropivot*_rs;
598: if (PetscRealPart(sctx->pv) <= _zero){
599: /* force matfactor to be diagonally dominant */
600: if (sctx->nshift == sctx->nshift_max) {
601: sctx->shift_fraction = sctx->shift_hi;
602: } else {
603: sctx->shift_lo = sctx->shift_fraction;
604: sctx->shift_fraction = (sctx->shift_hi+sctx->shift_lo)/2.;
605: }
606: sctx->shift_amount = sctx->shift_fraction * sctx->shift_top;
607: sctx->nshift++;
608: sctx->newshift = PETSC_TRUE;
609: } else {
610: sctx->newshift = PETSC_FALSE;
611: }
612: return(0);
613: }
617: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_inblocks(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
618: {
619: PetscReal _zero = info->zeropivot;
622: if (PetscAbsScalar(sctx->pv) <= _zero){
623: sctx->pv += info->shiftamount;
624: sctx->shift_amount = 0.0;
625: sctx->nshift++;
626: }
627: sctx->newshift = PETSC_FALSE;
628: return(0);
629: }
633: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_none(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
634: {
635: PetscReal _zero = info->zeropivot;
639: sctx->newshift = PETSC_FALSE;
640: if (PetscAbsScalar(sctx->pv) <= _zero) {
641: if (!mat->erroriffpe) {
642: PetscInfo3(mat,"Detected zero pivot in factorization in row %D value %g tolerance %g",row,(double)PetscAbsScalar(sctx->pv),(double)_zero);
643: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %g tolerance %g",row,(double)PetscAbsScalar(sctx->pv),(double)_zero);
644: }
645: return(0);
646: }
650: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
651: {
655: if (info->shifttype == (PetscReal) MAT_SHIFT_NONZERO){
656: MatPivotCheck_nz(mat,info,sctx,row);
657: } else if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE){
658: MatPivotCheck_pd(mat,info,sctx,row);
659: } else if (info->shifttype == (PetscReal) MAT_SHIFT_INBLOCKS){
660: MatPivotCheck_inblocks(mat,info,sctx,row);
661: } else {
662: MatPivotCheck_none(mat,info,sctx,row);
663: }
664: return(0);
665: }
667: /*
668: Create and initialize a linked list
669: Input Parameters:
670: idx_start - starting index of the list
671: lnk_max - max value of lnk indicating the end of the list
672: nlnk - max length of the list
673: Output Parameters:
674: lnk - list initialized
675: bt - PetscBT (bitarray) with all bits set to false
676: lnk_empty - flg indicating the list is empty
677: */
678: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
679: (PetscMalloc1(nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk[idx_start] = lnk_max,0))
681: #define PetscLLCreate_new(idx_start,lnk_max,nlnk,lnk,bt,lnk_empty)\
682: (PetscMalloc1(nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk_empty = PETSC_TRUE,0) ||(lnk[idx_start] = lnk_max,0))
684: /*
685: Add an index set into a sorted linked list
686: Input Parameters:
687: nidx - number of input indices
688: indices - interger array
689: idx_start - starting index of the list
690: lnk - linked list(an integer array) that is created
691: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
692: output Parameters:
693: nlnk - number of newly added indices
694: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
695: bt - updated PetscBT (bitarray)
696: */
697: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
698: {\
699: PetscInt _k,_entry,_location,_lnkdata;\
700: nlnk = 0;\
701: _lnkdata = idx_start;\
702: for (_k=0; _k<nidx; _k++){\
703: _entry = indices[_k];\
704: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
705: /* search for insertion location */\
706: /* start from the beginning if _entry < previous _entry */\
707: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
708: do {\
709: _location = _lnkdata;\
710: _lnkdata = lnk[_location];\
711: } while (_entry > _lnkdata);\
712: /* insertion location is found, add entry into lnk */\
713: lnk[_location] = _entry;\
714: lnk[_entry] = _lnkdata;\
715: nlnk++;\
716: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
717: }\
718: }\
719: }
721: /*
722: Add a permuted index set into a sorted linked list
723: Input Parameters:
724: nidx - number of input indices
725: indices - interger array
726: perm - permutation of indices
727: idx_start - starting index of the list
728: lnk - linked list(an integer array) that is created
729: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
730: output Parameters:
731: nlnk - number of newly added indices
732: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
733: bt - updated PetscBT (bitarray)
734: */
735: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
736: {\
737: PetscInt _k,_entry,_location,_lnkdata;\
738: nlnk = 0;\
739: _lnkdata = idx_start;\
740: for (_k=0; _k<nidx; _k++){\
741: _entry = perm[indices[_k]];\
742: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
743: /* search for insertion location */\
744: /* start from the beginning if _entry < previous _entry */\
745: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
746: do {\
747: _location = _lnkdata;\
748: _lnkdata = lnk[_location];\
749: } while (_entry > _lnkdata);\
750: /* insertion location is found, add entry into lnk */\
751: lnk[_location] = _entry;\
752: lnk[_entry] = _lnkdata;\
753: nlnk++;\
754: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
755: }\
756: }\
757: }
759: /*
760: Add a SORTED ascending index set into a sorted linked list - same as PetscLLAdd() bus skip 'if (_k && _entry < _lnkdata) _lnkdata = idx_start;'
761: Input Parameters:
762: nidx - number of input indices
763: indices - sorted interger array
764: idx_start - starting index of the list
765: lnk - linked list(an integer array) that is created
766: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
767: output Parameters:
768: nlnk - number of newly added indices
769: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
770: bt - updated PetscBT (bitarray)
771: */
772: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
773: {\
774: PetscInt _k,_entry,_location,_lnkdata;\
775: nlnk = 0;\
776: _lnkdata = idx_start;\
777: for (_k=0; _k<nidx; _k++){\
778: _entry = indices[_k];\
779: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
780: /* search for insertion location */\
781: do {\
782: _location = _lnkdata;\
783: _lnkdata = lnk[_location];\
784: } while (_entry > _lnkdata);\
785: /* insertion location is found, add entry into lnk */\
786: lnk[_location] = _entry;\
787: lnk[_entry] = _lnkdata;\
788: nlnk++;\
789: _lnkdata = _entry; /* next search starts from here */\
790: }\
791: }\
792: }
794: #define PetscLLAddSorted_new(nidx,indices,idx_start,lnk_empty,nlnk,lnk,bt) 0; \
795: {\
796: PetscInt _k,_entry,_location,_lnkdata;\
797: if (lnk_empty){\
798: _lnkdata = idx_start; \
799: for (_k=0; _k<nidx; _k++){ \
800: _entry = indices[_k]; \
801: PetscBTSet(bt,_entry); /* mark the new entry */ \
802: _location = _lnkdata; \
803: _lnkdata = lnk[_location]; \
804: /* insertion location is found, add entry into lnk */ \
805: lnk[_location] = _entry; \
806: lnk[_entry] = _lnkdata; \
807: _lnkdata = _entry; /* next search starts from here */ \
808: } \
809: /*\
810: lnk[indices[nidx-1]] = lnk[idx_start];\
811: lnk[idx_start] = indices[0];\
812: PetscBTSet(bt,indices[0]); \
813: for (_k=1; _k<nidx; _k++){ \
814: PetscBTSet(bt,indices[_k]); \
815: lnk[indices[_k-1]] = indices[_k]; \
816: } \
817: */\
818: nlnk = nidx;\
819: lnk_empty = PETSC_FALSE;\
820: } else {\
821: nlnk = 0; \
822: _lnkdata = idx_start; \
823: for (_k=0; _k<nidx; _k++){ \
824: _entry = indices[_k]; \
825: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */ \
826: /* search for insertion location */ \
827: do { \
828: _location = _lnkdata; \
829: _lnkdata = lnk[_location]; \
830: } while (_entry > _lnkdata); \
831: /* insertion location is found, add entry into lnk */ \
832: lnk[_location] = _entry; \
833: lnk[_entry] = _lnkdata; \
834: nlnk++; \
835: _lnkdata = _entry; /* next search starts from here */ \
836: } \
837: } \
838: } \
839: }
841: /*
842: Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
843: Same as PetscLLAddSorted() with an additional operation:
844: count the number of input indices that are no larger than 'diag'
845: Input Parameters:
846: indices - sorted interger array
847: idx_start - starting index of the list, index of pivot row
848: lnk - linked list(an integer array) that is created
849: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
850: diag - index of the active row in LUFactorSymbolic
851: nzbd - number of input indices with indices <= idx_start
852: im - im[idx_start] is initialized as num of nonzero entries in row=idx_start
853: output Parameters:
854: nlnk - number of newly added indices
855: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
856: bt - updated PetscBT (bitarray)
857: im - im[idx_start]: unchanged if diag is not an entry
858: : num of entries with indices <= diag if diag is an entry
859: */
860: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
861: {\
862: PetscInt _k,_entry,_location,_lnkdata,_nidx;\
863: nlnk = 0;\
864: _lnkdata = idx_start;\
865: _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
866: for (_k=0; _k<_nidx; _k++){\
867: _entry = indices[_k];\
868: nzbd++;\
869: if ( _entry== diag) im[idx_start] = nzbd;\
870: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
871: /* search for insertion location */\
872: do {\
873: _location = _lnkdata;\
874: _lnkdata = lnk[_location];\
875: } while (_entry > _lnkdata);\
876: /* insertion location is found, add entry into lnk */\
877: lnk[_location] = _entry;\
878: lnk[_entry] = _lnkdata;\
879: nlnk++;\
880: _lnkdata = _entry; /* next search starts from here */\
881: }\
882: }\
883: }
885: /*
886: Copy data on the list into an array, then initialize the list
887: Input Parameters:
888: idx_start - starting index of the list
889: lnk_max - max value of lnk indicating the end of the list
890: nlnk - number of data on the list to be copied
891: lnk - linked list
892: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
893: output Parameters:
894: indices - array that contains the copied data
895: lnk - linked list that is cleaned and initialize
896: bt - PetscBT (bitarray) with all bits set to false
897: */
898: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
899: {\
900: PetscInt _j,_idx=idx_start;\
901: for (_j=0; _j<nlnk; _j++){\
902: _idx = lnk[_idx];\
903: indices[_j] = _idx;\
904: PetscBTClear(bt,_idx);\
905: }\
906: lnk[idx_start] = lnk_max;\
907: }
908: /*
909: Free memories used by the list
910: */
911: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(&(bt)))
913: /* Routines below are used for incomplete matrix factorization */
914: /*
915: Create and initialize a linked list and its levels
916: Input Parameters:
917: idx_start - starting index of the list
918: lnk_max - max value of lnk indicating the end of the list
919: nlnk - max length of the list
920: Output Parameters:
921: lnk - list initialized
922: lnk_lvl - array of size nlnk for storing levels of lnk
923: bt - PetscBT (bitarray) with all bits set to false
924: */
925: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
926: (PetscMalloc1(2*nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))
928: /*
929: Initialize a sorted linked list used for ILU and ICC
930: Input Parameters:
931: nidx - number of input idx
932: idx - interger array used for storing column indices
933: idx_start - starting index of the list
934: perm - indices of an IS
935: lnk - linked list(an integer array) that is created
936: lnklvl - levels of lnk
937: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
938: output Parameters:
939: nlnk - number of newly added idx
940: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
941: lnklvl - levels of lnk
942: bt - updated PetscBT (bitarray)
943: */
944: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
945: {\
946: PetscInt _k,_entry,_location,_lnkdata;\
947: nlnk = 0;\
948: _lnkdata = idx_start;\
949: for (_k=0; _k<nidx; _k++){\
950: _entry = perm[idx[_k]];\
951: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
952: /* search for insertion location */\
953: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
954: do {\
955: _location = _lnkdata;\
956: _lnkdata = lnk[_location];\
957: } while (_entry > _lnkdata);\
958: /* insertion location is found, add entry into lnk */\
959: lnk[_location] = _entry;\
960: lnk[_entry] = _lnkdata;\
961: lnklvl[_entry] = 0;\
962: nlnk++;\
963: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
964: }\
965: }\
966: }
968: /*
969: Add a SORTED index set into a sorted linked list for ILU
970: Input Parameters:
971: nidx - number of input indices
972: idx - sorted interger array used for storing column indices
973: level - level of fill, e.g., ICC(level)
974: idxlvl - level of idx
975: idx_start - starting index of the list
976: lnk - linked list(an integer array) that is created
977: lnklvl - levels of lnk
978: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
979: prow - the row number of idx
980: output Parameters:
981: nlnk - number of newly added idx
982: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
983: lnklvl - levels of lnk
984: bt - updated PetscBT (bitarray)
986: Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
987: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
988: */
989: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
990: {\
991: PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
992: nlnk = 0;\
993: _lnkdata = idx_start;\
994: for (_k=0; _k<nidx; _k++){\
995: _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
996: if (_incrlev > level) continue;\
997: _entry = idx[_k];\
998: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
999: /* search for insertion location */\
1000: do {\
1001: _location = _lnkdata;\
1002: _lnkdata = lnk[_location];\
1003: } while (_entry > _lnkdata);\
1004: /* insertion location is found, add entry into lnk */\
1005: lnk[_location] = _entry;\
1006: lnk[_entry] = _lnkdata;\
1007: lnklvl[_entry] = _incrlev;\
1008: nlnk++;\
1009: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
1010: } else { /* existing entry: update lnklvl */\
1011: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1012: }\
1013: }\
1014: }
1016: /*
1017: Add a index set into a sorted linked list
1018: Input Parameters:
1019: nidx - number of input idx
1020: idx - interger array used for storing column indices
1021: level - level of fill, e.g., ICC(level)
1022: idxlvl - level of idx
1023: idx_start - starting index of the list
1024: lnk - linked list(an integer array) that is created
1025: lnklvl - levels of lnk
1026: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1027: output Parameters:
1028: nlnk - number of newly added idx
1029: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1030: lnklvl - levels of lnk
1031: bt - updated PetscBT (bitarray)
1032: */
1033: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
1034: {\
1035: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1036: nlnk = 0;\
1037: _lnkdata = idx_start;\
1038: for (_k=0; _k<nidx; _k++){\
1039: _incrlev = idxlvl[_k] + 1;\
1040: if (_incrlev > level) continue;\
1041: _entry = idx[_k];\
1042: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1043: /* search for insertion location */\
1044: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
1045: do {\
1046: _location = _lnkdata;\
1047: _lnkdata = lnk[_location];\
1048: } while (_entry > _lnkdata);\
1049: /* insertion location is found, add entry into lnk */\
1050: lnk[_location] = _entry;\
1051: lnk[_entry] = _lnkdata;\
1052: lnklvl[_entry] = _incrlev;\
1053: nlnk++;\
1054: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
1055: } else { /* existing entry: update lnklvl */\
1056: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1057: }\
1058: }\
1059: }
1061: /*
1062: Add a SORTED index set into a sorted linked list
1063: Input Parameters:
1064: nidx - number of input indices
1065: idx - sorted interger array used for storing column indices
1066: level - level of fill, e.g., ICC(level)
1067: idxlvl - level of idx
1068: idx_start - starting index of the list
1069: lnk - linked list(an integer array) that is created
1070: lnklvl - levels of lnk
1071: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1072: output Parameters:
1073: nlnk - number of newly added idx
1074: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1075: lnklvl - levels of lnk
1076: bt - updated PetscBT (bitarray)
1077: */
1078: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
1079: {\
1080: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1081: nlnk = 0;\
1082: _lnkdata = idx_start;\
1083: for (_k=0; _k<nidx; _k++){\
1084: _incrlev = idxlvl[_k] + 1;\
1085: if (_incrlev > level) continue;\
1086: _entry = idx[_k];\
1087: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1088: /* search for insertion location */\
1089: do {\
1090: _location = _lnkdata;\
1091: _lnkdata = lnk[_location];\
1092: } while (_entry > _lnkdata);\
1093: /* insertion location is found, add entry into lnk */\
1094: lnk[_location] = _entry;\
1095: lnk[_entry] = _lnkdata;\
1096: lnklvl[_entry] = _incrlev;\
1097: nlnk++;\
1098: _lnkdata = _entry; /* next search starts from here */\
1099: } else { /* existing entry: update lnklvl */\
1100: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1101: }\
1102: }\
1103: }
1105: /*
1106: Add a SORTED index set into a sorted linked list for ICC
1107: Input Parameters:
1108: nidx - number of input indices
1109: idx - sorted interger array used for storing column indices
1110: level - level of fill, e.g., ICC(level)
1111: idxlvl - level of idx
1112: idx_start - starting index of the list
1113: lnk - linked list(an integer array) that is created
1114: lnklvl - levels of lnk
1115: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1116: idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
1117: output Parameters:
1118: nlnk - number of newly added indices
1119: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1120: lnklvl - levels of lnk
1121: bt - updated PetscBT (bitarray)
1122: Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
1123: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
1124: */
1125: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
1126: {\
1127: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1128: nlnk = 0;\
1129: _lnkdata = idx_start;\
1130: for (_k=0; _k<nidx; _k++){\
1131: _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
1132: if (_incrlev > level) continue;\
1133: _entry = idx[_k];\
1134: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1135: /* search for insertion location */\
1136: do {\
1137: _location = _lnkdata;\
1138: _lnkdata = lnk[_location];\
1139: } while (_entry > _lnkdata);\
1140: /* insertion location is found, add entry into lnk */\
1141: lnk[_location] = _entry;\
1142: lnk[_entry] = _lnkdata;\
1143: lnklvl[_entry] = _incrlev;\
1144: nlnk++;\
1145: _lnkdata = _entry; /* next search starts from here */\
1146: } else { /* existing entry: update lnklvl */\
1147: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1148: }\
1149: }\
1150: }
1152: /*
1153: Copy data on the list into an array, then initialize the list
1154: Input Parameters:
1155: idx_start - starting index of the list
1156: lnk_max - max value of lnk indicating the end of the list
1157: nlnk - number of data on the list to be copied
1158: lnk - linked list
1159: lnklvl - level of lnk
1160: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1161: output Parameters:
1162: indices - array that contains the copied data
1163: lnk - linked list that is cleaned and initialize
1164: lnklvl - level of lnk that is reinitialized
1165: bt - PetscBT (bitarray) with all bits set to false
1166: */
1167: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
1168: {\
1169: PetscInt _j,_idx=idx_start;\
1170: for (_j=0; _j<nlnk; _j++){\
1171: _idx = lnk[_idx];\
1172: *(indices+_j) = _idx;\
1173: *(indiceslvl+_j) = lnklvl[_idx];\
1174: lnklvl[_idx] = -1;\
1175: PetscBTClear(bt,_idx);\
1176: }\
1177: lnk[idx_start] = lnk_max;\
1178: }
1179: /*
1180: Free memories used by the list
1181: */
1182: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(&(bt)))
1184: /* -------------------------------------------------------------------------------------------------------*/
1185: #include <petscbt.h>
1188: /*
1189: Create and initialize a condensed linked list -
1190: same as PetscLLCreate(), but uses a scalable array 'lnk' with size of max number of entries, not O(N).
1191: Barry suggested this approach (Dec. 6, 2011):
1192: I've thought of an alternative way of representing a linked list that is efficient but doesn't have the O(N) scaling issue
1193: (it may be faster than the O(N) even sequentially due to less crazy memory access).
1195: Instead of having some like a 2 -> 4 -> 11 -> 22 list that uses slot 2 4 11 and 22 in a big array use a small array with two slots
1196: for each entry for example [ 2 1 | 4 3 | 22 -1 | 11 2] so the first number (of the pair) is the value while the second tells you where
1197: in the list the next entry is. Inserting a new link means just append another pair at the end. For example say we want to insert 13 into the
1198: list it would then become [2 1 | 4 3 | 22 -1 | 11 4 | 13 2 ] you just add a pair at the end and fix the point for the one that points to it.
1199: That is 11 use to point to the 2 slot, after the change 11 points to the 4th slot which has the value 13. Note that values are always next
1200: to each other so memory access is much better than using the big array.
1202: Example:
1203: nlnk_max=5, lnk_max=36:
1204: Initial list: [0, 0 | 36, 2 | 0, 0 | 0, 0 | 0, 0 | 0, 0 | 0, 0]
1205: here, head_node has index 2 with value lnk[2]=lnk_max=36,
1206: 0-th entry is used to store the number of entries in the list,
1207: The initial lnk represents head -> tail(marked by 36) with number of entries = lnk[0]=0.
1209: Now adding a sorted set {2,4}, the list becomes
1210: [2, 0 | 36, 4 |2, 6 | 4, 2 | 0, 0 | 0, 0 | 0, 0 ]
1211: represents head -> 2 -> 4 -> tail with number of entries = lnk[0]=2.
1213: Then adding a sorted set {0,3,35}, the list
1214: [5, 0 | 36, 8 | 2, 10 | 4, 12 | 0, 4 | 3, 6 | 35, 2 ]
1215: represents head -> 0 -> 2 -> 3 -> 4 -> 35 -> tail with number of entries = lnk[0]=5.
1217: Input Parameters:
1218: nlnk_max - max length of the list
1219: lnk_max - max value of the entries
1220: Output Parameters:
1221: lnk - list created and initialized
1222: bt - PetscBT (bitarray) with all bits set to false. Note: bt has size lnk_max, not nln_max!
1223: */
1224: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate(PetscInt nlnk_max,PetscInt lnk_max,PetscInt **lnk,PetscBT *bt)
1225: {
1227: PetscInt *llnk;
1230: PetscMalloc1(2*(nlnk_max+2),lnk);
1231: PetscBTCreate(lnk_max,bt);
1232: llnk = *lnk;
1233: llnk[0] = 0; /* number of entries on the list */
1234: llnk[2] = lnk_max; /* value in the head node */
1235: llnk[3] = 2; /* next for the head node */
1236: return(0);
1237: }
1241: /*
1242: Add a SORTED ascending index set into a sorted linked list. See PetscLLCondensedCreate() for detailed description.
1243: Input Parameters:
1244: nidx - number of input indices
1245: indices - sorted interger array
1246: lnk - condensed linked list(an integer array) that is created
1247: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1248: output Parameters:
1249: lnk - the sorted(increasing order) linked list containing previous and newly added non-redundate indices
1250: bt - updated PetscBT (bitarray)
1251: */
1252: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted(PetscInt nidx,const PetscInt indices[],PetscInt lnk[],PetscBT bt)
1253: {
1254: PetscInt _k,_entry,_location,_next,_lnkdata,_nlnk,_newnode;
1257: _nlnk = lnk[0]; /* num of entries on the input lnk */
1258: _location = 2; /* head */
1259: for (_k=0; _k<nidx; _k++){
1260: _entry = indices[_k];
1261: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */
1262: /* search for insertion location */
1263: do {
1264: _next = _location + 1; /* link from previous node to next node */
1265: _location = lnk[_next]; /* idx of next node */
1266: _lnkdata = lnk[_location];/* value of next node */
1267: } while (_entry > _lnkdata);
1268: /* insertion location is found, add entry into lnk */
1269: _newnode = 2*(_nlnk+2); /* index for this new node */
1270: lnk[_next] = _newnode; /* connect previous node to the new node */
1271: lnk[_newnode] = _entry; /* set value of the new node */
1272: lnk[_newnode+1] = _location; /* connect new node to next node */
1273: _location = _newnode; /* next search starts from the new node */
1274: _nlnk++;
1275: } \
1276: }\
1277: lnk[0] = _nlnk; /* number of entries in the list */
1278: return(0);
1279: }
1283: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean(PetscInt lnk_max,PetscInt nidx,PetscInt *indices,PetscInt lnk[],PetscBT bt)
1284: {
1286: PetscInt _k,_next,_nlnk;
1289: _next = lnk[3]; /* head node */
1290: _nlnk = lnk[0]; /* num of entries on the list */
1291: for (_k=0; _k<_nlnk; _k++){
1292: indices[_k] = lnk[_next];
1293: _next = lnk[_next + 1];
1294: PetscBTClear(bt,indices[_k]);
1295: }
1296: lnk[0] = 0; /* num of entries on the list */
1297: lnk[2] = lnk_max; /* initialize head node */
1298: lnk[3] = 2; /* head node */
1299: return(0);
1300: }
1304: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedView(PetscInt *lnk)
1305: {
1307: PetscInt k;
1310: PetscPrintf(PETSC_COMM_SELF,"LLCondensed of size %d, (val, next)\n",lnk[0]);
1311: for (k=2; k< lnk[0]+2; k++){
1312: PetscPrintf(PETSC_COMM_SELF," %D: (%D, %D)\n",2*k,lnk[2*k],lnk[2*k+1]);
1313: }
1314: return(0);
1315: }
1319: /*
1320: Free memories used by the list
1321: */
1322: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy(PetscInt *lnk,PetscBT bt)
1323: {
1327: PetscFree(lnk);
1328: PetscBTDestroy(&bt);
1329: return(0);
1330: }
1332: /* -------------------------------------------------------------------------------------------------------*/
1335: /*
1336: Same as PetscLLCondensedCreate(), but does not use non-scalable O(lnk_max) bitarray
1337: Input Parameters:
1338: nlnk_max - max length of the list
1339: Output Parameters:
1340: lnk - list created and initialized
1341: */
1342: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate_Scalable(PetscInt nlnk_max,PetscInt **lnk)
1343: {
1345: PetscInt *llnk;
1348: PetscMalloc1(2*(nlnk_max+2),lnk);
1349: llnk = *lnk;
1350: llnk[0] = 0; /* number of entries on the list */
1351: llnk[2] = PETSC_MAX_INT; /* value in the head node */
1352: llnk[3] = 2; /* next for the head node */
1353: return(0);
1354: }
1358: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted_Scalable(PetscInt nidx,const PetscInt indices[],PetscInt lnk[])
1359: {
1360: PetscInt _k,_entry,_location,_next,_lnkdata,_nlnk,_newnode;
1361: _nlnk = lnk[0]; /* num of entries on the input lnk */
1362: _location = 2; /* head */ \
1363: for (_k=0; _k<nidx; _k++){
1364: _entry = indices[_k];
1365: /* search for insertion location */
1366: do {
1367: _next = _location + 1; /* link from previous node to next node */
1368: _location = lnk[_next]; /* idx of next node */
1369: _lnkdata = lnk[_location];/* value of next node */
1370: } while (_entry > _lnkdata);
1371: if (_entry < _lnkdata) {
1372: /* insertion location is found, add entry into lnk */
1373: _newnode = 2*(_nlnk+2); /* index for this new node */
1374: lnk[_next] = _newnode; /* connect previous node to the new node */
1375: lnk[_newnode] = _entry; /* set value of the new node */
1376: lnk[_newnode+1] = _location; /* connect new node to next node */
1377: _location = _newnode; /* next search starts from the new node */
1378: _nlnk++;
1379: }
1380: }
1381: lnk[0] = _nlnk; /* number of entries in the list */
1382: return 0;
1383: }
1387: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean_Scalable(PetscInt nidx,PetscInt *indices,PetscInt *lnk)
1388: {
1389: PetscInt _k,_next,_nlnk;
1390: _next = lnk[3]; /* head node */
1391: _nlnk = lnk[0];
1392: for (_k=0; _k<_nlnk; _k++){
1393: indices[_k] = lnk[_next];
1394: _next = lnk[_next + 1];
1395: }
1396: lnk[0] = 0; /* num of entries on the list */
1397: lnk[3] = 2; /* head node */
1398: return 0;
1399: }
1403: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy_Scalable(PetscInt *lnk)
1404: {
1405: return PetscFree(lnk);
1406: }
1408: /* -------------------------------------------------------------------------------------------------------*/
1409: /*
1410: lnk[0] number of links
1411: lnk[1] number of entries
1412: lnk[3n] value
1413: lnk[3n+1] len
1414: lnk[3n+2] link to next value
1416: The next three are always the first link
1418: lnk[3] PETSC_MIN_INT+1
1419: lnk[4] 1
1420: lnk[5] link to first real entry
1422: The next three are always the last link
1424: lnk[6] PETSC_MAX_INT - 1
1425: lnk[7] 1
1426: lnk[8] next valid link (this is the same as lnk[0] but without the decreases)
1427: */
1431: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate_fast(PetscInt nlnk_max,PetscInt **lnk)
1432: {
1434: PetscInt *llnk;
1437: PetscMalloc1(3*(nlnk_max+3),lnk);
1438: llnk = *lnk;
1439: llnk[0] = 0; /* nlnk: number of entries on the list */
1440: llnk[1] = 0; /* number of integer entries represented in list */
1441: llnk[3] = PETSC_MIN_INT+1; /* value in the first node */
1442: llnk[4] = 1; /* count for the first node */
1443: llnk[5] = 6; /* next for the first node */
1444: llnk[6] = PETSC_MAX_INT-1; /* value in the last node */
1445: llnk[7] = 1; /* count for the last node */
1446: llnk[8] = 0; /* next valid node to be used */
1447: return(0);
1448: }
1450: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted_fast(PetscInt nidx,const PetscInt indices[],PetscInt lnk[])
1451: {
1452: PetscInt k,entry,prev,next;
1453: prev = 3; /* first value */
1454: next = lnk[prev+2];
1455: for (k=0; k<nidx; k++){
1456: entry = indices[k];
1457: /* search for insertion location */
1458: while (entry >= lnk[next]) {
1459: prev = next;
1460: next = lnk[next+2];
1461: }
1462: /* entry is in range of previous list */
1463: if (entry < lnk[prev]+lnk[prev+1]) continue;
1464: lnk[1]++;
1465: /* entry is right after previous list */
1466: if (entry == lnk[prev]+lnk[prev+1]) {
1467: lnk[prev+1]++;
1468: if (lnk[next] == entry+1) { /* combine two contiquous strings */
1469: lnk[prev+1] += lnk[next+1];
1470: lnk[prev+2] = lnk[next+2];
1471: next = lnk[next+2];
1472: lnk[0]--;
1473: }
1474: continue;
1475: }
1476: /* entry is right before next list */
1477: if (entry == lnk[next]-1) {
1478: lnk[next]--;
1479: lnk[next+1]++;
1480: prev = next;
1481: next = lnk[prev+2];
1482: continue;
1483: }
1484: /* add entry into lnk */
1485: lnk[prev+2] = 3*((lnk[8]++)+3); /* connect previous node to the new node */
1486: prev = lnk[prev+2];
1487: lnk[prev] = entry; /* set value of the new node */
1488: lnk[prev+1] = 1; /* number of values in contiquous string is one to start */
1489: lnk[prev+2] = next; /* connect new node to next node */
1490: lnk[0]++;
1491: }
1492: return 0;
1493: }
1495: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean_fast(PetscInt nidx,PetscInt *indices,PetscInt *lnk)
1496: {
1497: PetscInt _k,_next,_nlnk,cnt,j;
1498: _next = lnk[5]; /* first node */
1499: _nlnk = lnk[0];
1500: cnt = 0;
1501: for (_k=0; _k<_nlnk; _k++){
1502: for (j=0; j<lnk[_next+1]; j++) {
1503: indices[cnt++] = lnk[_next] + j;
1504: }
1505: _next = lnk[_next + 2];
1506: }
1507: lnk[0] = 0; /* nlnk: number of links */
1508: lnk[1] = 0; /* number of integer entries represented in list */
1509: lnk[3] = PETSC_MIN_INT+1; /* value in the first node */
1510: lnk[4] = 1; /* count for the first node */
1511: lnk[5] = 6; /* next for the first node */
1512: lnk[6] = PETSC_MAX_INT-1; /* value in the last node */
1513: lnk[7] = 1; /* count for the last node */
1514: lnk[8] = 0; /* next valid location to make link */
1515: return 0;
1516: }
1518: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedView_fast(PetscInt *lnk)
1519: {
1520: PetscInt k,next,nlnk;
1521: next = lnk[5]; /* first node */
1522: nlnk = lnk[0];
1523: for (k=0; k<nlnk; k++){
1524: #if 0 /* Debugging code */
1525: printf("%d value %d len %d next %d\n",next,lnk[next],lnk[next+1],lnk[next+2]);
1526: #endif
1527: next = lnk[next + 2];
1528: }
1529: return 0;
1530: }
1532: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy_fast(PetscInt *lnk)
1533: {
1534: return PetscFree(lnk);
1535: }
1537: /* alias PetscSortIntWithScalarArray while MatScalar == PetscScalar */
1538: PETSC_STATIC_INLINE PetscErrorCode PetscSortIntWithMatScalarArray(PetscInt n,PetscInt *idx,PetscScalar *val)
1539: {
1540: #if !defined(PETSC_USE_REAL_MAT_SINGLE)
1541: return PetscSortIntWithScalarArray(n,idx,val);
1542: #else
1543: {
1544: MatScalar mtmp;
1545: return PetscSortIntWithDataArray(n,idx,val,sizeof(MatScalar),&mtmp);
1546: }
1547: #endif
1548: }
1550: PETSC_EXTERN PetscLogEvent MAT_Mult, MAT_MultMatrixFree, MAT_Mults, MAT_MultConstrained, MAT_MultAdd, MAT_MultTranspose;
1551: PETSC_EXTERN PetscLogEvent MAT_MultTransposeConstrained, MAT_MultTransposeAdd, MAT_Solve, MAT_Solves, MAT_SolveAdd, MAT_SolveTranspose;
1552: PETSC_EXTERN PetscLogEvent MAT_SolveTransposeAdd, MAT_SOR, MAT_ForwardSolve, MAT_BackwardSolve, MAT_LUFactor, MAT_LUFactorSymbolic;
1553: PETSC_EXTERN PetscLogEvent MAT_LUFactorNumeric, MAT_CholeskyFactor, MAT_CholeskyFactorSymbolic, MAT_CholeskyFactorNumeric, MAT_ILUFactor;
1554: PETSC_EXTERN PetscLogEvent MAT_ILUFactorSymbolic, MAT_ICCFactorSymbolic, MAT_Copy, MAT_Convert, MAT_Scale, MAT_AssemblyBegin;
1555: PETSC_EXTERN PetscLogEvent MAT_AssemblyEnd, MAT_SetValues, MAT_GetValues, MAT_GetRow, MAT_GetRowIJ, MAT_GetSubMatrices, MAT_GetColoring, MAT_GetOrdering, MAT_RedundantMat;
1556: PETSC_EXTERN PetscLogEvent MAT_IncreaseOverlap, MAT_Partitioning, MAT_Coarsen, MAT_ZeroEntries, MAT_Load, MAT_View, MAT_AXPY, MAT_FDColoringCreate, MAT_TransposeColoringCreate;
1557: PETSC_EXTERN PetscLogEvent MAT_FDColoringSetUp, MAT_FDColoringApply, MAT_Transpose, MAT_FDColoringFunction,MAT_GetSubMatrix;
1558: PETSC_EXTERN PetscLogEvent MAT_MatMult, MAT_MatSolve,MAT_MatMultSymbolic, MAT_MatMultNumeric,MAT_Getlocalmatcondensed,MAT_GetBrowsOfAcols,MAT_GetBrowsOfAocols;
1559: PETSC_EXTERN PetscLogEvent MAT_PtAP, MAT_PtAPSymbolic, MAT_PtAPNumeric,MAT_Seqstompinum,MAT_Seqstompisym,MAT_Seqstompi,MAT_Getlocalmat;
1560: PETSC_EXTERN PetscLogEvent MAT_RARt, MAT_RARtSymbolic, MAT_RARtNumeric;
1561: PETSC_EXTERN PetscLogEvent MAT_MatTransposeMult, MAT_MatTransposeMultSymbolic, MAT_MatTransposeMultNumeric;
1562: PETSC_EXTERN PetscLogEvent MAT_TransposeMatMult, MAT_TransposeMatMultSymbolic, MAT_TransposeMatMultNumeric;
1563: PETSC_EXTERN PetscLogEvent MAT_MatMatMult, MAT_MatMatMultSymbolic, MAT_MatMatMultNumeric;
1564: PETSC_EXTERN PetscLogEvent MAT_Applypapt, MAT_Applypapt_symbolic, MAT_Applypapt_numeric;
1565: PETSC_EXTERN PetscLogEvent MAT_Getsymtranspose, MAT_Transpose_SeqAIJ, MAT_Getsymtransreduced,MAT_GetSequentialNonzeroStructure;
1567: PETSC_EXTERN PetscLogEvent MATMFFD_Mult;
1568: PETSC_EXTERN PetscLogEvent MAT_GetMultiProcBlock;
1569: PETSC_EXTERN PetscLogEvent MAT_CUSPCopyToGPU, MAT_CUSPARSECopyToGPU, MAT_SetValuesBatch, MAT_SetValuesBatchI, MAT_SetValuesBatchII, MAT_SetValuesBatchIII, MAT_SetValuesBatchIV;
1570: PETSC_EXTERN PetscLogEvent MAT_ViennaCLCopyToGPU;
1571: PETSC_EXTERN PetscLogEvent MAT_Merge,MAT_Residual;
1572: PETSC_EXTERN PetscLogEvent Mat_Coloring_Apply,Mat_Coloring_Comm,Mat_Coloring_Local,Mat_Coloring_ISCreate,Mat_Coloring_SetUp,Mat_Coloring_Weights;
1574: #endif