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