Actual source code: matimpl.h
petsc-3.15.0 2021-03-30
2: #ifndef __MATIMPL_H
5: #include <petscmat.h>
6: #include <petscmatcoarsen.h>
7: #include <petsc/private/petscimpl.h>
9: PETSC_EXTERN PetscBool MatRegisterAllCalled;
10: PETSC_EXTERN PetscBool MatSeqAIJRegisterAllCalled;
11: PETSC_EXTERN PetscBool MatOrderingRegisterAllCalled;
12: PETSC_EXTERN PetscBool MatColoringRegisterAllCalled;
13: PETSC_EXTERN PetscBool MatPartitioningRegisterAllCalled;
14: PETSC_EXTERN PetscBool MatCoarsenRegisterAllCalled;
15: PETSC_EXTERN PetscErrorCode MatRegisterAll(void);
16: PETSC_EXTERN PetscErrorCode MatOrderingRegisterAll(void);
17: PETSC_EXTERN PetscErrorCode MatColoringRegisterAll(void);
18: PETSC_EXTERN PetscErrorCode MatPartitioningRegisterAll(void);
19: PETSC_EXTERN PetscErrorCode MatCoarsenRegisterAll(void);
20: PETSC_EXTERN PetscErrorCode MatSeqAIJRegisterAll(void);
22: /*
23: This file defines the parts of the matrix data structure that are
24: shared by all matrix types.
25: */
27: /*
28: If you add entries here also add them to the MATOP enum
29: in include/petscmat.h and src/mat/f90-mod/petscmat.h
30: */
31: typedef struct _MatOps *MatOps;
32: struct _MatOps {
33: /* 0*/
34: PetscErrorCode (*setvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
35: PetscErrorCode (*getrow)(Mat,PetscInt,PetscInt *,PetscInt*[],PetscScalar*[]);
36: PetscErrorCode (*restorerow)(Mat,PetscInt,PetscInt *,PetscInt *[],PetscScalar *[]);
37: PetscErrorCode (*mult)(Mat,Vec,Vec);
38: PetscErrorCode (*multadd)(Mat,Vec,Vec,Vec);
39: /* 5*/
40: PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
41: PetscErrorCode (*multtransposeadd)(Mat,Vec,Vec,Vec);
42: PetscErrorCode (*solve)(Mat,Vec,Vec);
43: PetscErrorCode (*solveadd)(Mat,Vec,Vec,Vec);
44: PetscErrorCode (*solvetranspose)(Mat,Vec,Vec);
45: /*10*/
46: PetscErrorCode (*solvetransposeadd)(Mat,Vec,Vec,Vec);
47: PetscErrorCode (*lufactor)(Mat,IS,IS,const MatFactorInfo*);
48: PetscErrorCode (*choleskyfactor)(Mat,IS,const MatFactorInfo*);
49: PetscErrorCode (*sor)(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
50: PetscErrorCode (*transpose)(Mat,MatReuse,Mat*);
51: /*15*/
52: PetscErrorCode (*getinfo)(Mat,MatInfoType,MatInfo*);
53: PetscErrorCode (*equal)(Mat,Mat,PetscBool*);
54: PetscErrorCode (*getdiagonal)(Mat,Vec);
55: PetscErrorCode (*diagonalscale)(Mat,Vec,Vec);
56: PetscErrorCode (*norm)(Mat,NormType,PetscReal*);
57: /*20*/
58: PetscErrorCode (*assemblybegin)(Mat,MatAssemblyType);
59: PetscErrorCode (*assemblyend)(Mat,MatAssemblyType);
60: PetscErrorCode (*setoption)(Mat,MatOption,PetscBool);
61: PetscErrorCode (*zeroentries)(Mat);
62: /*24*/
63: PetscErrorCode (*zerorows)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
64: PetscErrorCode (*lufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
65: PetscErrorCode (*lufactornumeric)(Mat,Mat,const MatFactorInfo*);
66: PetscErrorCode (*choleskyfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
67: PetscErrorCode (*choleskyfactornumeric)(Mat,Mat,const MatFactorInfo*);
68: /*29*/
69: PetscErrorCode (*setup)(Mat);
70: PetscErrorCode (*ilufactorsymbolic)(Mat,Mat,IS,IS,const MatFactorInfo*);
71: PetscErrorCode (*iccfactorsymbolic)(Mat,Mat,IS,const MatFactorInfo*);
72: PetscErrorCode (*getdiagonalblock)(Mat,Mat*);
73: PetscErrorCode (*setinf)(Mat);
74: /*34*/
75: PetscErrorCode (*duplicate)(Mat,MatDuplicateOption,Mat*);
76: PetscErrorCode (*forwardsolve)(Mat,Vec,Vec);
77: PetscErrorCode (*backwardsolve)(Mat,Vec,Vec);
78: PetscErrorCode (*ilufactor)(Mat,IS,IS,const MatFactorInfo*);
79: PetscErrorCode (*iccfactor)(Mat,IS,const MatFactorInfo*);
80: /*39*/
81: PetscErrorCode (*axpy)(Mat,PetscScalar,Mat,MatStructure);
82: PetscErrorCode (*createsubmatrices)(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat *[]);
83: PetscErrorCode (*increaseoverlap)(Mat,PetscInt,IS[],PetscInt);
84: PetscErrorCode (*getvalues)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscScalar []);
85: PetscErrorCode (*copy)(Mat,Mat,MatStructure);
86: /*44*/
87: PetscErrorCode (*getrowmax)(Mat,Vec,PetscInt[]);
88: PetscErrorCode (*scale)(Mat,PetscScalar);
89: PetscErrorCode (*shift)(Mat,PetscScalar);
90: PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
91: PetscErrorCode (*zerorowscolumns)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
92: /*49*/
93: PetscErrorCode (*setrandom)(Mat,PetscRandom);
94: PetscErrorCode (*getrowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
95: PetscErrorCode (*restorerowij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt *,const PetscInt *[],const PetscInt *[],PetscBool *);
96: PetscErrorCode (*getcolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
97: PetscErrorCode (*restorecolumnij)(Mat,PetscInt,PetscBool ,PetscBool ,PetscInt*,const PetscInt *[],const PetscInt *[],PetscBool *);
98: /*54*/
99: PetscErrorCode (*fdcoloringcreate)(Mat,ISColoring,MatFDColoring);
100: PetscErrorCode (*coloringpatch)(Mat,PetscInt,PetscInt,ISColoringValue[],ISColoring*);
101: PetscErrorCode (*setunfactored)(Mat);
102: PetscErrorCode (*permute)(Mat,IS,IS,Mat*);
103: PetscErrorCode (*setvaluesblocked)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
104: /*59*/
105: PetscErrorCode (*createsubmatrix)(Mat,IS,IS,MatReuse,Mat*);
106: PetscErrorCode (*destroy)(Mat);
107: PetscErrorCode (*view)(Mat,PetscViewer);
108: PetscErrorCode (*convertfrom)(Mat,MatType,MatReuse,Mat*);
109: PetscErrorCode (*placeholder_63)(void);
110: /*64*/
111: PetscErrorCode (*matmatmultsymbolic)(Mat,Mat,Mat,PetscReal,Mat);
112: PetscErrorCode (*matmatmultnumeric)(Mat,Mat,Mat,Mat);
113: PetscErrorCode (*setlocaltoglobalmapping)(Mat,ISLocalToGlobalMapping,ISLocalToGlobalMapping);
114: PetscErrorCode (*setvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
115: PetscErrorCode (*zerorowslocal)(Mat,PetscInt,const PetscInt[],PetscScalar,Vec,Vec);
116: /*69*/
117: PetscErrorCode (*getrowmaxabs)(Mat,Vec,PetscInt[]);
118: PetscErrorCode (*getrowminabs)(Mat,Vec,PetscInt[]);
119: PetscErrorCode (*convert)(Mat, MatType,MatReuse,Mat*);
120: PetscErrorCode (*hasoperation)(Mat,MatOperation,PetscBool*);
121: PetscErrorCode (*placeholder_73)(void);
122: /*74*/
123: PetscErrorCode (*setvaluesadifor)(Mat,PetscInt,void*);
124: PetscErrorCode (*fdcoloringapply)(Mat,MatFDColoring,Vec,void*);
125: PetscErrorCode (*setfromoptions)(PetscOptionItems*,Mat);
126: PetscErrorCode (*multconstrained)(Mat,Vec,Vec);
127: PetscErrorCode (*multtransposeconstrained)(Mat,Vec,Vec);
128: /*79*/
129: PetscErrorCode (*findzerodiagonals)(Mat,IS*);
130: PetscErrorCode (*mults)(Mat,Vecs,Vecs);
131: PetscErrorCode (*solves)(Mat,Vecs,Vecs);
132: PetscErrorCode (*getinertia)(Mat,PetscInt*,PetscInt*,PetscInt*);
133: PetscErrorCode (*load)(Mat,PetscViewer);
134: /*84*/
135: PetscErrorCode (*issymmetric)(Mat,PetscReal,PetscBool*);
136: PetscErrorCode (*ishermitian)(Mat,PetscReal,PetscBool*);
137: PetscErrorCode (*isstructurallysymmetric)(Mat,PetscBool *);
138: PetscErrorCode (*setvaluesblockedlocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
139: PetscErrorCode (*getvecs)(Mat,Vec*,Vec*);
140: /*89*/
141: PetscErrorCode (*placeholder_89)(void);
142: PetscErrorCode (*matmultsymbolic)(Mat,Mat,PetscReal,Mat);
143: PetscErrorCode (*matmultnumeric)(Mat,Mat,Mat);
144: PetscErrorCode (*placeholder_92)(void);
145: PetscErrorCode (*ptapsymbolic)(Mat,Mat,PetscReal,Mat); /* double dispatch wrapper routine */
146: /*94*/
147: PetscErrorCode (*ptapnumeric)(Mat,Mat,Mat); /* double dispatch wrapper routine */
148: PetscErrorCode (*placeholder_95)(void);
149: PetscErrorCode (*mattransposemultsymbolic)(Mat,Mat,PetscReal,Mat);
150: PetscErrorCode (*mattransposemultnumeric)(Mat,Mat,Mat);
151: PetscErrorCode (*bindtocpu)(Mat,PetscBool);
152: /*99*/
153: PetscErrorCode (*productsetfromoptions)(Mat);
154: PetscErrorCode (*productsymbolic)(Mat);
155: PetscErrorCode (*productnumeric)(Mat);
156: PetscErrorCode (*conjugate)(Mat); /* complex conjugate */
157: PetscErrorCode (*viewnative)(Mat,PetscViewer);
158: /*104*/
159: PetscErrorCode (*setvaluesrow)(Mat,PetscInt,const PetscScalar[]);
160: PetscErrorCode (*realpart)(Mat);
161: PetscErrorCode (*imaginarypart)(Mat);
162: PetscErrorCode (*getrowuppertriangular)(Mat);
163: PetscErrorCode (*restorerowuppertriangular)(Mat);
164: /*109*/
165: PetscErrorCode (*matsolve)(Mat,Mat,Mat);
166: PetscErrorCode (*matsolvetranspose)(Mat,Mat,Mat);
167: PetscErrorCode (*getrowmin)(Mat,Vec,PetscInt[]);
168: PetscErrorCode (*getcolumnvector)(Mat,Vec,PetscInt);
169: PetscErrorCode (*missingdiagonal)(Mat,PetscBool *,PetscInt*);
170: /*114*/
171: PetscErrorCode (*getseqnonzerostructure)(Mat,Mat *);
172: PetscErrorCode (*create)(Mat);
173: PetscErrorCode (*getghosts)(Mat,PetscInt*,const PetscInt *[]);
174: PetscErrorCode (*getlocalsubmatrix)(Mat,IS,IS,Mat*);
175: PetscErrorCode (*restorelocalsubmatrix)(Mat,IS,IS,Mat*);
176: /*119*/
177: PetscErrorCode (*multdiagonalblock)(Mat,Vec,Vec);
178: PetscErrorCode (*hermitiantranspose)(Mat,MatReuse,Mat*);
179: PetscErrorCode (*multhermitiantranspose)(Mat,Vec,Vec);
180: PetscErrorCode (*multhermitiantransposeadd)(Mat,Vec,Vec,Vec);
181: PetscErrorCode (*getmultiprocblock)(Mat,MPI_Comm,MatReuse,Mat*);
182: /*124*/
183: PetscErrorCode (*findnonzerorows)(Mat,IS*);
184: PetscErrorCode (*getcolumnnorms)(Mat,NormType,PetscReal*);
185: PetscErrorCode (*invertblockdiagonal)(Mat,const PetscScalar**);
186: PetscErrorCode (*invertvariableblockdiagonal)(Mat,PetscInt,const PetscInt*,PetscScalar*);
187: PetscErrorCode (*createsubmatricesmpi)(Mat,PetscInt,const IS[], const IS[], MatReuse, Mat**);
188: /*129*/
189: PetscErrorCode (*setvaluesbatch)(Mat,PetscInt,PetscInt,PetscInt*,const PetscScalar*);
190: PetscErrorCode (*placeholder_130)(void);
191: PetscErrorCode (*transposematmultsymbolic)(Mat,Mat,PetscReal,Mat);
192: PetscErrorCode (*transposematmultnumeric)(Mat,Mat,Mat);
193: PetscErrorCode (*transposecoloringcreate)(Mat,ISColoring,MatTransposeColoring);
194: /*134*/
195: PetscErrorCode (*transcoloringapplysptoden)(MatTransposeColoring,Mat,Mat);
196: PetscErrorCode (*transcoloringapplydentosp)(MatTransposeColoring,Mat,Mat);
197: PetscErrorCode (*placeholder_136)(void);
198: PetscErrorCode (*rartsymbolic)(Mat,Mat,PetscReal,Mat); /* double dispatch wrapper routine */
199: PetscErrorCode (*rartnumeric)(Mat,Mat,Mat); /* double dispatch wrapper routine */
200: /*139*/
201: PetscErrorCode (*setblocksizes)(Mat,PetscInt,PetscInt);
202: PetscErrorCode (*aypx)(Mat,PetscScalar,Mat,MatStructure);
203: PetscErrorCode (*residual)(Mat,Vec,Vec,Vec);
204: PetscErrorCode (*fdcoloringsetup)(Mat,ISColoring,MatFDColoring);
205: PetscErrorCode (*findoffblockdiagonalentries)(Mat,IS*);
206: PetscErrorCode (*creatempimatconcatenateseqmat)(MPI_Comm,Mat,PetscInt,MatReuse,Mat*);
207: /*145*/
208: PetscErrorCode (*destroysubmatrices)(PetscInt,Mat*[]);
209: PetscErrorCode (*mattransposesolve)(Mat,Mat,Mat);
210: PetscErrorCode (*getvalueslocal)(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscScalar[]);
211: };
212: /*
213: If you add MatOps entries above also add them to the MATOP enum
214: in include/petscmat.h and src/mat/f90-mod/petscmat.h
215: */
217: #include <petscsys.h>
218: PETSC_EXTERN PetscErrorCode MatRegisterOp(MPI_Comm, const char[], PetscVoidFunction, const char[], PetscInt, ...);
219: PETSC_EXTERN PetscErrorCode MatQueryOp(MPI_Comm, PetscVoidFunction*, const char[], PetscInt, ...);
221: typedef struct _p_MatRootName* MatRootName;
222: struct _p_MatRootName {
223: char *rname,*sname,*mname;
224: MatRootName next;
225: };
227: PETSC_EXTERN MatRootName MatRootNameList;
229: /*
230: Utility private matrix routines
231: */
232: PETSC_INTERN PetscErrorCode MatFindNonzeroRowsOrCols_Basic(Mat,PetscBool,PetscReal,IS*);
233: PETSC_INTERN PetscErrorCode MatConvert_Basic(Mat,MatType,MatReuse,Mat*);
234: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat,MatType,MatReuse,Mat*);
235: PETSC_INTERN PetscErrorCode MatConvertFrom_Shell(Mat,MatType,MatReuse,Mat*);
236: PETSC_INTERN PetscErrorCode MatCopy_Basic(Mat,Mat,MatStructure);
237: PETSC_INTERN PetscErrorCode MatDiagonalSet_Default(Mat,Vec,InsertMode);
238: #if defined(PETSC_HAVE_SCALAPACK)
239: PETSC_INTERN PetscErrorCode MatConvert_Dense_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
240: #endif
241: PETSC_INTERN PetscErrorCode MatSetPreallocationCOO_Basic(Mat,PetscInt,const PetscInt[],const PetscInt[]);
242: PETSC_INTERN PetscErrorCode MatSetValuesCOO_Basic(Mat,const PetscScalar[],InsertMode);
244: /* these callbacks rely on the old matrix function pointers for
245: matmat operations. They are unsafe, and should be removed.
246: However, the amount of work needed to clean up all the
247: implementations is not negligible */
248: PETSC_INTERN PetscErrorCode MatProductSymbolic_AB(Mat);
249: PETSC_INTERN PetscErrorCode MatProductNumeric_AB(Mat);
250: PETSC_INTERN PetscErrorCode MatProductSymbolic_AtB(Mat);
251: PETSC_INTERN PetscErrorCode MatProductNumeric_AtB(Mat);
252: PETSC_INTERN PetscErrorCode MatProductSymbolic_ABt(Mat);
253: PETSC_INTERN PetscErrorCode MatProductNumeric_ABt(Mat);
254: PETSC_INTERN PetscErrorCode MatProductNumeric_PtAP(Mat);
255: PETSC_INTERN PetscErrorCode MatProductNumeric_RARt(Mat);
256: PETSC_INTERN PetscErrorCode MatProductSymbolic_ABC(Mat);
257: PETSC_INTERN PetscErrorCode MatProductNumeric_ABC(Mat);
259: PETSC_INTERN PetscErrorCode MatProductCreate_Private(Mat,Mat,Mat,Mat);
260: /* this callback handles all the different triple products and
261: does not rely on the function pointers; used by cuSPARSE and KOKKOS-KERNELS */
262: PETSC_INTERN PetscErrorCode MatProductSymbolic_ABC_Basic(Mat);
265: #if defined(PETSC_USE_DEBUG)
266: # define MatCheckPreallocated(A,arg) do { \
267: 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); \
268: } while (0)
269: #else
270: # define MatCheckPreallocated(A,arg) do {} while (0)
271: #endif
273: #if defined(PETSC_USE_DEBUG)
274: # define MatCheckProduct(A,arg) do { \
275: if (PetscUnlikely(!(A)->product)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Argument %D \"%s\" is not a matrix obtained from MatProductCreate()",(arg),#A); \
276: } while (0)
277: #else
278: # define MatCheckProduct(A,arg) do {} while (0)
279: #endif
281: /*
282: The stash is used to temporarily store inserted matrix values that
283: belong to another processor. During the assembly phase the stashed
284: values are moved to the correct processor and
285: */
287: typedef struct _MatStashSpace *PetscMatStashSpace;
289: struct _MatStashSpace {
290: PetscMatStashSpace next;
291: PetscScalar *space_head,*val;
292: PetscInt *idx,*idy;
293: PetscInt total_space_size;
294: PetscInt local_used;
295: PetscInt local_remaining;
296: };
298: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceGet(PetscInt,PetscInt,PetscMatStashSpace *);
299: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceContiguous(PetscInt,PetscMatStashSpace *,PetscScalar *,PetscInt *,PetscInt *);
300: PETSC_EXTERN PetscErrorCode PetscMatStashSpaceDestroy(PetscMatStashSpace*);
302: typedef struct {
303: PetscInt count;
304: } MatStashHeader;
306: typedef struct {
307: void *buffer; /* Of type blocktype, dynamically constructed */
308: PetscInt count;
309: char pending;
310: } MatStashFrame;
312: typedef struct _MatStash MatStash;
313: struct _MatStash {
314: PetscInt nmax; /* maximum stash size */
315: PetscInt umax; /* user specified max-size */
316: PetscInt oldnmax; /* the nmax value used previously */
317: PetscInt n; /* stash size */
318: PetscInt bs; /* block size of the stash */
319: PetscInt reallocs; /* preserve the no of mallocs invoked */
320: PetscMatStashSpace space_head,space; /* linked list to hold stashed global row/column numbers and matrix values */
322: PetscErrorCode (*ScatterBegin)(Mat,MatStash*,PetscInt*);
323: PetscErrorCode (*ScatterGetMesg)(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,PetscScalar**,PetscInt*);
324: PetscErrorCode (*ScatterEnd)(MatStash*);
325: PetscErrorCode (*ScatterDestroy)(MatStash*);
327: /* The following variables are used for communication */
328: MPI_Comm comm;
329: PetscMPIInt size,rank;
330: PetscMPIInt tag1,tag2;
331: MPI_Request *send_waits; /* array of send requests */
332: MPI_Request *recv_waits; /* array of receive requests */
333: MPI_Status *send_status; /* array of send status */
334: PetscInt nsends,nrecvs; /* numbers of sends and receives */
335: PetscScalar *svalues; /* sending data */
336: PetscInt *sindices;
337: PetscScalar **rvalues; /* receiving data (values) */
338: PetscInt **rindices; /* receiving data (indices) */
339: PetscInt nprocessed; /* number of messages already processed */
340: PetscMPIInt *flg_v; /* indicates what messages have arrived so far and from whom */
341: PetscBool reproduce;
342: PetscInt reproduce_count;
344: /* The following variables are used for BTS communication */
345: PetscBool first_assembly_done; /* Is the first time matrix assembly done? */
346: PetscBool use_status; /* Use MPI_Status to determine number of items in each message */
347: PetscMPIInt nsendranks;
348: PetscMPIInt nrecvranks;
349: PetscMPIInt *sendranks;
350: PetscMPIInt *recvranks;
351: MatStashHeader *sendhdr,*recvhdr;
352: MatStashFrame *sendframes; /* pointers to the main messages */
353: MatStashFrame *recvframes;
354: MatStashFrame *recvframe_active;
355: PetscInt recvframe_i; /* index of block within active frame */
356: PetscMPIInt recvframe_count; /* Count actually sent for current frame */
357: PetscInt recvcount; /* Number of receives processed so far */
358: PetscMPIInt *some_indices; /* From last call to MPI_Waitsome */
359: MPI_Status *some_statuses; /* Statuses from last call to MPI_Waitsome */
360: PetscMPIInt some_count; /* Number of requests completed in last call to MPI_Waitsome */
361: PetscMPIInt some_i; /* Index of request currently being processed */
362: MPI_Request *sendreqs;
363: MPI_Request *recvreqs;
364: PetscSegBuffer segsendblocks;
365: PetscSegBuffer segrecvframe;
366: PetscSegBuffer segrecvblocks;
367: MPI_Datatype blocktype;
368: size_t blocktype_size;
369: InsertMode *insertmode; /* Pointer to check mat->insertmode and set upon message arrival in case no local values have been set. */
370: };
372: #if !defined(PETSC_HAVE_MPIUNI)
373: PETSC_INTERN PetscErrorCode MatStashScatterDestroy_BTS(MatStash*);
374: #endif
375: PETSC_INTERN PetscErrorCode MatStashCreate_Private(MPI_Comm,PetscInt,MatStash*);
376: PETSC_INTERN PetscErrorCode MatStashDestroy_Private(MatStash*);
377: PETSC_INTERN PetscErrorCode MatStashScatterEnd_Private(MatStash*);
378: PETSC_INTERN PetscErrorCode MatStashSetInitialSize_Private(MatStash*,PetscInt);
379: PETSC_INTERN PetscErrorCode MatStashGetInfo_Private(MatStash*,PetscInt*,PetscInt*);
380: PETSC_INTERN PetscErrorCode MatStashValuesRow_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscBool);
381: PETSC_INTERN PetscErrorCode MatStashValuesCol_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscBool);
382: PETSC_INTERN PetscErrorCode MatStashValuesRowBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
383: PETSC_INTERN PetscErrorCode MatStashValuesColBlocked_Private(MatStash*,PetscInt,PetscInt,const PetscInt[],const PetscScalar[],PetscInt,PetscInt,PetscInt);
384: PETSC_INTERN PetscErrorCode MatStashScatterBegin_Private(Mat,MatStash*,PetscInt*);
385: PETSC_INTERN PetscErrorCode MatStashScatterGetMesg_Private(MatStash*,PetscMPIInt*,PetscInt**,PetscInt**,PetscScalar**,PetscInt*);
386: PETSC_INTERN PetscErrorCode MatGetInfo_External(Mat,MatInfoType,MatInfo*);
388: typedef struct {
389: PetscInt dim;
390: PetscInt dims[4];
391: PetscInt starts[4];
392: PetscBool noc; /* this is a single component problem, hence user will not set MatStencil.c */
393: } MatStencilInfo;
395: /* Info about using compressed row format */
396: typedef struct {
397: PetscBool use; /* indicates compressed rows have been checked and will be used */
398: PetscInt nrows; /* number of non-zero rows */
399: PetscInt *i; /* compressed row pointer */
400: PetscInt *rindex; /* compressed row index */
401: } Mat_CompressedRow;
402: PETSC_EXTERN PetscErrorCode MatCheckCompressedRow(Mat,PetscInt,Mat_CompressedRow*,PetscInt*,PetscInt,PetscReal);
404: typedef struct { /* used by MatCreateRedundantMatrix() for reusing matredundant */
405: PetscInt nzlocal,nsends,nrecvs;
406: PetscMPIInt *send_rank,*recv_rank;
407: PetscInt *sbuf_nz,*rbuf_nz,*sbuf_j,**rbuf_j;
408: PetscScalar *sbuf_a,**rbuf_a;
409: MPI_Comm subcomm; /* when user does not provide a subcomm */
410: IS isrow,iscol;
411: Mat *matseq;
412: } Mat_Redundant;
414: typedef struct { /* used by MatProduct() */
415: MatProductType type;
416: char *alg;
417: Mat A,B,C,Dwork;
418: PetscReal fill;
419: PetscBool api_user; /* used to distinguish command line options and to indicate the matrix values are ready to be consumed at symbolic phase if needed */
421: /* Some products may display the information on the algorithm used */
422: PetscErrorCode (*view)(Mat,PetscViewer);
424: /* many products have intermediate data structures, each specific to Mat types and product type */
425: PetscBool clear; /* whether or not to clear the data structures after MatProductNumeric has been called */
426: void *data; /* where to stash those structures */
427: PetscErrorCode (*destroy)(void*); /* destroy routine */
428: } Mat_Product;
430: #define CSRDataStructure(datatype) \
431: int *i; \
432: int *j; \
433: datatype *a;\
434: PetscInt n;\
435: PetscInt ignorezeroentries;
437: typedef struct {
438: CSRDataStructure(PetscScalar)
439: } PetscCSRDataStructure;
441: struct _p_SplitCSRMat {
442: PetscInt cstart,cend,rstart,rend;
443: PetscCSRDataStructure diag,offdiag;
444: PetscInt *colmap;
445: PetscMPIInt rank;
446: };
448: struct _p_Mat {
449: PETSCHEADER(struct _MatOps);
450: PetscLayout rmap,cmap;
451: void *data; /* implementation-specific data */
452: MatFactorType factortype; /* MAT_FACTOR_LU, ILU, CHOLESKY or ICC */
453: PetscBool useordering; /* factorization using ordering provide to routine (most PETSc implementations) */
454: PetscBool assembled; /* is the matrix assembled? */
455: PetscBool was_assembled; /* new values inserted into assembled mat */
456: PetscInt num_ass; /* number of times matrix has been assembled */
457: PetscObjectState nonzerostate; /* each time new nonzeros locations are introduced into the matrix this is updated */
458: PetscObjectState ass_nonzerostate; /* nonzero state at last assembly */
459: MatInfo info; /* matrix information */
460: InsertMode insertmode; /* have values been inserted in matrix or added? */
461: MatStash stash,bstash; /* used for assembling off-proc mat emements */
462: MatNullSpace nullsp; /* null space (operator is singular) */
463: MatNullSpace transnullsp; /* null space of transpose of operator */
464: MatNullSpace nearnullsp; /* near null space to be used by multigrid methods */
465: PetscInt congruentlayouts; /* are the rows and columns layouts congruent? */
466: PetscBool preallocated;
467: MatStencilInfo stencil; /* information for structured grid */
468: PetscBool symmetric,hermitian,structurally_symmetric,spd;
469: PetscBool symmetric_set,hermitian_set,structurally_symmetric_set,spd_set; /* if true, then corresponding flag is correct*/
470: PetscBool symmetric_eternal;
471: PetscBool nooffprocentries,nooffproczerorows;
472: PetscBool assembly_subset; /* set by MAT_SUBSET_OFF_PROC_ENTRIES */
473: PetscBool submat_singleis; /* for efficient PCSetUp_ASM() */
474: PetscBool structure_only;
475: PetscBool sortedfull; /* full, sorted rows are inserted */
476: PetscBool force_diagonals; /* set by MAT_FORCE_DIAGONAL_ENTRIES */
477: #if defined(PETSC_HAVE_DEVICE)
478: PetscOffloadMask offloadmask; /* a mask which indicates where the valid matrix data is (GPU, CPU or both) */
479: PetscBool boundtocpu;
480: #endif
481: void *spptr; /* pointer for special library like SuperLU */
482: char *solvertype;
483: PetscBool checksymmetryonassembly,checknullspaceonassembly;
484: PetscReal checksymmetrytol;
485: Mat schur; /* Schur complement matrix */
486: MatFactorSchurStatus schur_status; /* status of the Schur complement matrix */
487: Mat_Redundant *redundant; /* used by MatCreateRedundantMatrix() */
488: PetscBool erroriffailure; /* Generate an error if detected (for example a zero pivot) instead of returning */
489: MatFactorError factorerrortype; /* type of error in factorization */
490: PetscReal factorerror_zeropivot_value; /* If numerical zero pivot was detected this is the computed value */
491: PetscInt factorerror_zeropivot_row; /* Row where zero pivot was detected */
492: PetscInt nblocks,*bsizes; /* support for MatSetVariableBlockSizes() */
493: char *defaultvectype;
494: Mat_Product *product;
495: PetscBool form_explicit_transpose; /* hint to generate an explicit mat tranpsose for operations like MatMultTranspose() */
496: PetscBool transupdated; /* whether or not the explicitly generated transpose is up-to-date */
497: };
499: PETSC_INTERN PetscErrorCode MatAXPY_Basic(Mat,PetscScalar,Mat,MatStructure);
500: PETSC_INTERN PetscErrorCode MatAXPY_BasicWithPreallocation(Mat,Mat,PetscScalar,Mat,MatStructure);
501: PETSC_INTERN PetscErrorCode MatAXPY_Basic_Preallocate(Mat,Mat,Mat*);
502: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat,PetscScalar,Mat);
504: /*
505: Utility for MatFactor (Schur complement)
506: */
507: PETSC_INTERN PetscErrorCode MatFactorFactorizeSchurComplement_Private(Mat);
508: PETSC_INTERN PetscErrorCode MatFactorInvertSchurComplement_Private(Mat);
509: PETSC_INTERN PetscErrorCode MatFactorUpdateSchurStatus_Private(Mat);
510: PETSC_INTERN PetscErrorCode MatFactorSetUpInPlaceSchur_Private(Mat);
512: /*
513: Utility for MatZeroRows
514: */
515: PETSC_INTERN PetscErrorCode MatZeroRowsMapLocal_Private(Mat,PetscInt,const PetscInt*,PetscInt*,PetscInt**);
517: /*
518: Utility for MatView/MatLoad
519: */
520: PETSC_INTERN PetscErrorCode MatView_Binary_BlockSizes(Mat,PetscViewer);
521: PETSC_INTERN PetscErrorCode MatLoad_Binary_BlockSizes(Mat,PetscViewer);
524: /*
525: Object for partitioning graphs
526: */
528: typedef struct _MatPartitioningOps *MatPartitioningOps;
529: struct _MatPartitioningOps {
530: PetscErrorCode (*apply)(MatPartitioning,IS*);
531: PetscErrorCode (*applynd)(MatPartitioning,IS*);
532: PetscErrorCode (*setfromoptions)(PetscOptionItems*,MatPartitioning);
533: PetscErrorCode (*destroy)(MatPartitioning);
534: PetscErrorCode (*view)(MatPartitioning,PetscViewer);
535: PetscErrorCode (*improve)(MatPartitioning,IS*);
536: };
538: struct _p_MatPartitioning {
539: PETSCHEADER(struct _MatPartitioningOps);
540: Mat adj;
541: PetscInt *vertex_weights;
542: PetscReal *part_weights;
543: PetscInt n; /* number of partitions */
544: void *data;
545: PetscInt setupcalled;
546: PetscBool use_edge_weights; /* A flag indicates whether or not to use edge weights */
547: };
549: /* needed for parallel nested dissection by ParMetis and PTSCOTCH */
550: PETSC_INTERN PetscErrorCode MatPartitioningSizesToSep_Private(PetscInt,PetscInt[],PetscInt[],PetscInt[]);
552: /*
553: Object for coarsen graphs
554: */
555: typedef struct _MatCoarsenOps *MatCoarsenOps;
556: struct _MatCoarsenOps {
557: PetscErrorCode (*apply)(MatCoarsen);
558: PetscErrorCode (*setfromoptions)(PetscOptionItems*,MatCoarsen);
559: PetscErrorCode (*destroy)(MatCoarsen);
560: PetscErrorCode (*view)(MatCoarsen,PetscViewer);
561: };
563: struct _p_MatCoarsen {
564: PETSCHEADER(struct _MatCoarsenOps);
565: Mat graph;
566: void *subctx;
567: /* */
568: PetscBool strict_aggs;
569: IS perm;
570: PetscCoarsenData *agg_lists;
571: };
573: /*
574: MatFDColoring is used to compute Jacobian matrices efficiently
575: via coloring. The data structure is explained below in an example.
577: Color = 0 1 0 2 | 2 3 0
578: ---------------------------------------------------
579: 00 01 | 05
580: 10 11 | 14 15 Processor 0
581: 22 23 | 25
582: 32 33 |
583: ===================================================
584: | 44 45 46
585: 50 | 55 Processor 1
586: | 64 66
587: ---------------------------------------------------
589: ncolors = 4;
591: ncolumns = {2,1,1,0}
592: columns = {{0,2},{1},{3},{}}
593: nrows = {4,2,3,3}
594: rows = {{0,1,2,3},{0,1},{1,2,3},{0,1,2}}
595: vwscale = {dx(0),dx(1),dx(2),dx(3)} MPI Vec
596: vscale = {dx(0),dx(1),dx(2),dx(3),dx(4),dx(5)} Seq Vec
598: ncolumns = {1,0,1,1}
599: columns = {{6},{},{4},{5}}
600: nrows = {3,0,2,2}
601: rows = {{0,1,2},{},{1,2},{1,2}}
602: vwscale = {dx(4),dx(5),dx(6)} MPI Vec
603: vscale = {dx(0),dx(4),dx(5),dx(6)} Seq Vec
605: See the routine MatFDColoringApply() for how this data is used
606: to compute the Jacobian.
608: */
609: typedef struct {
610: PetscInt row;
611: PetscInt col;
612: PetscScalar *valaddr; /* address of value */
613: } MatEntry;
615: typedef struct {
616: PetscInt row;
617: PetscScalar *valaddr; /* address of value */
618: } MatEntry2;
620: struct _p_MatFDColoring{
621: PETSCHEADER(int);
622: PetscInt M,N,m; /* total rows, columns; local rows */
623: PetscInt rstart; /* first row owned by local processor */
624: PetscInt ncolors; /* number of colors */
625: PetscInt *ncolumns; /* number of local columns for a color */
626: PetscInt **columns; /* lists the local columns of each color (using global column numbering) */
627: IS *isa; /* these are the IS that contain the column values given in columns */
628: PetscInt *nrows; /* number of local rows for each color */
629: MatEntry *matentry; /* holds (row, column, address of value) for Jacobian matrix entry */
630: MatEntry2 *matentry2; /* holds (row, address of value) for Jacobian matrix entry */
631: PetscScalar *dy; /* store a block of F(x+dx)-F(x) when J is in BAIJ format */
632: PetscReal error_rel; /* square root of relative error in computing function */
633: PetscReal umin; /* minimum allowable u'dx value */
634: Vec w1,w2,w3; /* work vectors used in computing Jacobian */
635: PetscBool fset; /* indicates that the initial function value F(X) is set */
636: PetscErrorCode (*f)(void); /* function that defines Jacobian */
637: void *fctx; /* optional user-defined context for use by the function f */
638: Vec vscale; /* holds FD scaling, i.e. 1/dx for each perturbed column */
639: PetscInt currentcolor; /* color for which function evaluation is being done now */
640: const char *htype; /* "wp" or "ds" */
641: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_LOCAL */
642: PetscInt brows,bcols; /* number of block rows or columns for speedup inserting the dense matrix into sparse Jacobian */
643: PetscBool setupcalled; /* true if setup has been called */
644: PetscBool viewed; /* true if the -mat_fd_coloring_view has been triggered already */
645: void (*ftn_func_pointer)(void),*ftn_func_cntx; /* serve the same purpose as *fortran_func_pointers in PETSc objects */
646: PetscObjectId matid; /* matrix this object was created with, must always be the same */
647: };
649: typedef struct _MatColoringOps *MatColoringOps;
650: struct _MatColoringOps {
651: PetscErrorCode (*destroy)(MatColoring);
652: PetscErrorCode (*setfromoptions)(PetscOptionItems*,MatColoring);
653: PetscErrorCode (*view)(MatColoring,PetscViewer);
654: PetscErrorCode (*apply)(MatColoring,ISColoring*);
655: PetscErrorCode (*weights)(MatColoring,PetscReal**,PetscInt**);
656: };
658: struct _p_MatColoring {
659: PETSCHEADER(struct _MatColoringOps);
660: Mat mat;
661: PetscInt dist; /* distance of the coloring */
662: PetscInt maxcolors; /* the maximum number of colors returned, maxcolors=1 for MIS */
663: void *data; /* inner context */
664: PetscBool valid; /* check to see if what is produced is a valid coloring */
665: MatColoringWeightType weight_type; /* type of weight computation to be performed */
666: PetscReal *user_weights; /* custom weights and permutation */
667: PetscInt *user_lperm;
668: PetscBool valid_iscoloring; /* check to see if matcoloring is produced a valid iscoloring */
669: };
671: struct _p_MatTransposeColoring{
672: PETSCHEADER(int);
673: PetscInt M,N,m; /* total rows, columns; local rows */
674: PetscInt rstart; /* first row owned by local processor */
675: PetscInt ncolors; /* number of colors */
676: PetscInt *ncolumns; /* number of local columns for a color */
677: PetscInt *nrows; /* number of local rows for each color */
678: PetscInt currentcolor; /* color for which function evaluation is being done now */
679: ISColoringType ctype; /* IS_COLORING_GLOBAL or IS_COLORING_LOCAL */
681: PetscInt *colorforrow,*colorforcol; /* pointer to rows and columns */
682: PetscInt *rows; /* lists the local rows for each color (using the local row numbering) */
683: PetscInt *den2sp; /* maps (row,color) in the dense matrix to index of sparse matrix array a->a */
684: PetscInt *columns; /* lists the local columns of each color (using global column numbering) */
685: PetscInt brows; /* number of rows for efficient implementation of MatTransColoringApplyDenToSp() */
686: PetscInt *lstart; /* array used for loop over row blocks of Csparse */
687: };
689: /*
690: Null space context for preconditioner/operators
691: */
692: struct _p_MatNullSpace {
693: PETSCHEADER(int);
694: PetscBool has_cnst;
695: PetscInt n;
696: Vec* vecs;
697: PetscScalar* alpha; /* for projections */
698: PetscErrorCode (*remove)(MatNullSpace,Vec,void*); /* for user provided removal function */
699: void* rmctx; /* context for remove() function */
700: };
702: /*
703: Checking zero pivot for LU, ILU preconditioners.
704: */
705: typedef struct {
706: PetscInt nshift,nshift_max;
707: PetscReal shift_amount,shift_lo,shift_hi,shift_top,shift_fraction;
708: PetscBool newshift;
709: PetscReal rs; /* active row sum of abs(offdiagonals) */
710: PetscScalar pv; /* pivot of the active row */
711: } FactorShiftCtx;
713: /*
714: Used by MatCreateSubMatrices_MPIXAIJ_Local()
715: */
716: #include <petscctable.h>
717: typedef struct { /* used by MatCreateSubMatrices_MPIAIJ_SingleIS_Local() and MatCreateSubMatrices_MPIAIJ_Local */
718: PetscInt id; /* index of submats, only submats[0] is responsible for deleting some arrays below */
719: PetscInt nrqs,nrqr;
720: PetscInt **rbuf1,**rbuf2,**rbuf3,**sbuf1,**sbuf2;
721: PetscInt **ptr;
722: PetscInt *tmp;
723: PetscInt *ctr;
724: PetscInt *pa; /* proc array */
725: PetscInt *req_size,*req_source1,*req_source2;
726: PetscBool allcolumns,allrows;
727: PetscBool singleis;
728: PetscInt *row2proc; /* row to proc map */
729: PetscInt nstages;
730: #if defined(PETSC_USE_CTABLE)
731: PetscTable cmap,rmap;
732: PetscInt *cmap_loc,*rmap_loc;
733: #else
734: PetscInt *cmap,*rmap;
735: #endif
737: PetscErrorCode (*destroy)(Mat);
738: } Mat_SubSppt;
740: PETSC_EXTERN PetscErrorCode MatFactorDumpMatrix(Mat);
741: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat,PetscScalar);
742: PETSC_INTERN PetscErrorCode MatSetBlockSizes_Default(Mat,PetscInt,PetscInt);
744: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_nz(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
745: {
746: PetscReal _rs = sctx->rs;
747: PetscReal _zero = info->zeropivot*_rs;
750: if (PetscAbsScalar(sctx->pv) <= _zero && !PetscIsNanScalar(sctx->pv)){
751: /* force |diag| > zeropivot*rs */
752: if (!sctx->nshift) sctx->shift_amount = info->shiftamount;
753: else sctx->shift_amount *= 2.0;
754: sctx->newshift = PETSC_TRUE;
755: (sctx->nshift)++;
756: } else {
757: sctx->newshift = PETSC_FALSE;
758: }
759: return(0);
760: }
762: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_pd(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
763: {
764: PetscReal _rs = sctx->rs;
765: PetscReal _zero = info->zeropivot*_rs;
768: if (PetscRealPart(sctx->pv) <= _zero && !PetscIsNanScalar(sctx->pv)){
769: /* force matfactor to be diagonally dominant */
770: if (sctx->nshift == sctx->nshift_max) {
771: sctx->shift_fraction = sctx->shift_hi;
772: } else {
773: sctx->shift_lo = sctx->shift_fraction;
774: sctx->shift_fraction = (sctx->shift_hi+sctx->shift_lo)/2.;
775: }
776: sctx->shift_amount = sctx->shift_fraction * sctx->shift_top;
777: sctx->nshift++;
778: sctx->newshift = PETSC_TRUE;
779: } else {
780: sctx->newshift = PETSC_FALSE;
781: }
782: return(0);
783: }
785: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_inblocks(Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
786: {
787: PetscReal _zero = info->zeropivot;
790: if (PetscAbsScalar(sctx->pv) <= _zero && !PetscIsNanScalar(sctx->pv)){
791: sctx->pv += info->shiftamount;
792: sctx->shift_amount = 0.0;
793: sctx->nshift++;
794: }
795: sctx->newshift = PETSC_FALSE;
796: return(0);
797: }
799: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck_none(Mat fact,Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
800: {
801: PetscReal _zero = info->zeropivot;
805: sctx->newshift = PETSC_FALSE;
806: if (PetscAbsScalar(sctx->pv) <= _zero && !PetscIsNanScalar(sctx->pv)) {
807: if (!mat->erroriffailure) {
808: PetscInfo3(mat,"Detected zero pivot in factorization in row %D value %g tolerance %g\n",row,(double)PetscAbsScalar(sctx->pv),(double)_zero);
809: fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
810: fact->factorerror_zeropivot_value = PetscAbsScalar(sctx->pv);
811: fact->factorerror_zeropivot_row = row;
812: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot row %D value %g tolerance %g\n",row,(double)PetscAbsScalar(sctx->pv),(double)_zero);
813: }
814: return(0);
815: }
817: PETSC_STATIC_INLINE PetscErrorCode MatPivotCheck(Mat fact,Mat mat,const MatFactorInfo *info,FactorShiftCtx *sctx,PetscInt row)
818: {
822: if (info->shifttype == (PetscReal) MAT_SHIFT_NONZERO){
823: MatPivotCheck_nz(mat,info,sctx,row);
824: } else if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE){
825: MatPivotCheck_pd(mat,info,sctx,row);
826: } else if (info->shifttype == (PetscReal) MAT_SHIFT_INBLOCKS){
827: MatPivotCheck_inblocks(mat,info,sctx,row);
828: } else {
829: MatPivotCheck_none(fact,mat,info,sctx,row);
830: }
831: return(0);
832: }
834: /*
835: Create and initialize a linked list
836: Input Parameters:
837: idx_start - starting index of the list
838: lnk_max - max value of lnk indicating the end of the list
839: nlnk - max length of the list
840: Output Parameters:
841: lnk - list initialized
842: bt - PetscBT (bitarray) with all bits set to false
843: lnk_empty - flg indicating the list is empty
844: */
845: #define PetscLLCreate(idx_start,lnk_max,nlnk,lnk,bt) \
846: (PetscMalloc1(nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk[idx_start] = lnk_max,0))
848: #define PetscLLCreate_new(idx_start,lnk_max,nlnk,lnk,bt,lnk_empty)\
849: (PetscMalloc1(nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk_empty = PETSC_TRUE,0) ||(lnk[idx_start] = lnk_max,0))
851: /*
852: Add an index set into a sorted linked list
853: Input Parameters:
854: nidx - number of input indices
855: indices - integer array
856: idx_start - starting index of the list
857: lnk - linked list(an integer array) that is created
858: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
859: output Parameters:
860: nlnk - number of newly added indices
861: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
862: bt - updated PetscBT (bitarray)
863: */
864: #define PetscLLAdd(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
865: {\
866: PetscInt _k,_entry,_location,_lnkdata;\
867: nlnk = 0;\
868: _lnkdata = idx_start;\
869: for (_k=0; _k<nidx; _k++){\
870: _entry = indices[_k];\
871: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
872: /* search for insertion location */\
873: /* start from the beginning if _entry < previous _entry */\
874: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
875: do {\
876: _location = _lnkdata;\
877: _lnkdata = lnk[_location];\
878: } while (_entry > _lnkdata);\
879: /* insertion location is found, add entry into lnk */\
880: lnk[_location] = _entry;\
881: lnk[_entry] = _lnkdata;\
882: nlnk++;\
883: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
884: }\
885: }\
886: }
888: /*
889: Add a permuted index set into a sorted linked list
890: Input Parameters:
891: nidx - number of input indices
892: indices - integer array
893: perm - permutation of indices
894: idx_start - starting index of the list
895: lnk - linked list(an integer array) that is created
896: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
897: output Parameters:
898: nlnk - number of newly added indices
899: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
900: bt - updated PetscBT (bitarray)
901: */
902: #define PetscLLAddPerm(nidx,indices,perm,idx_start,nlnk,lnk,bt) 0;\
903: {\
904: PetscInt _k,_entry,_location,_lnkdata;\
905: nlnk = 0;\
906: _lnkdata = idx_start;\
907: for (_k=0; _k<nidx; _k++){\
908: _entry = perm[indices[_k]];\
909: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
910: /* search for insertion location */\
911: /* start from the beginning if _entry < previous _entry */\
912: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
913: do {\
914: _location = _lnkdata;\
915: _lnkdata = lnk[_location];\
916: } while (_entry > _lnkdata);\
917: /* insertion location is found, add entry into lnk */\
918: lnk[_location] = _entry;\
919: lnk[_entry] = _lnkdata;\
920: nlnk++;\
921: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
922: }\
923: }\
924: }
926: /*
927: Add a SORTED ascending index set into a sorted linked list - same as PetscLLAdd() bus skip 'if (_k && _entry < _lnkdata) _lnkdata = idx_start;'
928: Input Parameters:
929: nidx - number of input indices
930: indices - sorted integer array
931: idx_start - starting index of the list
932: lnk - linked list(an integer array) that is created
933: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
934: output Parameters:
935: nlnk - number of newly added indices
936: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
937: bt - updated PetscBT (bitarray)
938: */
939: #define PetscLLAddSorted(nidx,indices,idx_start,nlnk,lnk,bt) 0;\
940: {\
941: PetscInt _k,_entry,_location,_lnkdata;\
942: nlnk = 0;\
943: _lnkdata = idx_start;\
944: for (_k=0; _k<nidx; _k++){\
945: _entry = indices[_k];\
946: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
947: /* search for insertion location */\
948: do {\
949: _location = _lnkdata;\
950: _lnkdata = lnk[_location];\
951: } while (_entry > _lnkdata);\
952: /* insertion location is found, add entry into lnk */\
953: lnk[_location] = _entry;\
954: lnk[_entry] = _lnkdata;\
955: nlnk++;\
956: _lnkdata = _entry; /* next search starts from here */\
957: }\
958: }\
959: }
961: #define PetscLLAddSorted_new(nidx,indices,idx_start,lnk_empty,nlnk,lnk,bt) 0; \
962: {\
963: PetscInt _k,_entry,_location,_lnkdata;\
964: if (lnk_empty){\
965: _lnkdata = idx_start; \
966: for (_k=0; _k<nidx; _k++){ \
967: _entry = indices[_k]; \
968: PetscBTSet(bt,_entry); /* mark the new entry */ \
969: _location = _lnkdata; \
970: _lnkdata = lnk[_location]; \
971: /* insertion location is found, add entry into lnk */ \
972: lnk[_location] = _entry; \
973: lnk[_entry] = _lnkdata; \
974: _lnkdata = _entry; /* next search starts from here */ \
975: } \
976: /*\
977: lnk[indices[nidx-1]] = lnk[idx_start];\
978: lnk[idx_start] = indices[0];\
979: PetscBTSet(bt,indices[0]); \
980: for (_k=1; _k<nidx; _k++){ \
981: PetscBTSet(bt,indices[_k]); \
982: lnk[indices[_k-1]] = indices[_k]; \
983: } \
984: */\
985: nlnk = nidx;\
986: lnk_empty = PETSC_FALSE;\
987: } else {\
988: nlnk = 0; \
989: _lnkdata = idx_start; \
990: for (_k=0; _k<nidx; _k++){ \
991: _entry = indices[_k]; \
992: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */ \
993: /* search for insertion location */ \
994: do { \
995: _location = _lnkdata; \
996: _lnkdata = lnk[_location]; \
997: } while (_entry > _lnkdata); \
998: /* insertion location is found, add entry into lnk */ \
999: lnk[_location] = _entry; \
1000: lnk[_entry] = _lnkdata; \
1001: nlnk++; \
1002: _lnkdata = _entry; /* next search starts from here */ \
1003: } \
1004: } \
1005: } \
1006: }
1008: /*
1009: Add a SORTED index set into a sorted linked list used for LUFactorSymbolic()
1010: Same as PetscLLAddSorted() with an additional operation:
1011: count the number of input indices that are no larger than 'diag'
1012: Input Parameters:
1013: indices - sorted integer array
1014: idx_start - starting index of the list, index of pivot row
1015: lnk - linked list(an integer array) that is created
1016: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1017: diag - index of the active row in LUFactorSymbolic
1018: nzbd - number of input indices with indices <= idx_start
1019: im - im[idx_start] is initialized as num of nonzero entries in row=idx_start
1020: output Parameters:
1021: nlnk - number of newly added indices
1022: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from indices
1023: bt - updated PetscBT (bitarray)
1024: im - im[idx_start]: unchanged if diag is not an entry
1025: : num of entries with indices <= diag if diag is an entry
1026: */
1027: #define PetscLLAddSortedLU(indices,idx_start,nlnk,lnk,bt,diag,nzbd,im) 0;\
1028: {\
1029: PetscInt _k,_entry,_location,_lnkdata,_nidx;\
1030: nlnk = 0;\
1031: _lnkdata = idx_start;\
1032: _nidx = im[idx_start] - nzbd; /* num of entries with idx_start < index <= diag */\
1033: for (_k=0; _k<_nidx; _k++){\
1034: _entry = indices[_k];\
1035: nzbd++;\
1036: if (_entry== diag) im[idx_start] = nzbd;\
1037: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1038: /* search for insertion location */\
1039: do {\
1040: _location = _lnkdata;\
1041: _lnkdata = lnk[_location];\
1042: } while (_entry > _lnkdata);\
1043: /* insertion location is found, add entry into lnk */\
1044: lnk[_location] = _entry;\
1045: lnk[_entry] = _lnkdata;\
1046: nlnk++;\
1047: _lnkdata = _entry; /* next search starts from here */\
1048: }\
1049: }\
1050: }
1052: /*
1053: Copy data on the list into an array, then initialize the list
1054: Input Parameters:
1055: idx_start - starting index of the list
1056: lnk_max - max value of lnk indicating the end of the list
1057: nlnk - number of data on the list to be copied
1058: lnk - linked list
1059: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1060: output Parameters:
1061: indices - array that contains the copied data
1062: lnk - linked list that is cleaned and initialize
1063: bt - PetscBT (bitarray) with all bits set to false
1064: */
1065: #define PetscLLClean(idx_start,lnk_max,nlnk,lnk,indices,bt) 0;\
1066: {\
1067: PetscInt _j,_idx=idx_start;\
1068: for (_j=0; _j<nlnk; _j++){\
1069: _idx = lnk[_idx];\
1070: indices[_j] = _idx;\
1071: PetscBTClear(bt,_idx);\
1072: }\
1073: lnk[idx_start] = lnk_max;\
1074: }
1075: /*
1076: Free memories used by the list
1077: */
1078: #define PetscLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(&(bt)))
1080: /* Routines below are used for incomplete matrix factorization */
1081: /*
1082: Create and initialize a linked list and its levels
1083: Input Parameters:
1084: idx_start - starting index of the list
1085: lnk_max - max value of lnk indicating the end of the list
1086: nlnk - max length of the list
1087: Output Parameters:
1088: lnk - list initialized
1089: lnk_lvl - array of size nlnk for storing levels of lnk
1090: bt - PetscBT (bitarray) with all bits set to false
1091: */
1092: #define PetscIncompleteLLCreate(idx_start,lnk_max,nlnk,lnk,lnk_lvl,bt)\
1093: (PetscIntMultError(2,nlnk,NULL) || PetscMalloc1(2*nlnk,&lnk) || PetscBTCreate(nlnk,&(bt)) || (lnk[idx_start] = lnk_max,lnk_lvl = lnk + nlnk,0))
1095: /*
1096: Initialize a sorted linked list used for ILU and ICC
1097: Input Parameters:
1098: nidx - number of input idx
1099: idx - integer array used for storing column indices
1100: idx_start - starting index of the list
1101: perm - indices of an IS
1102: lnk - linked list(an integer array) that is created
1103: lnklvl - levels of lnk
1104: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1105: output Parameters:
1106: nlnk - number of newly added idx
1107: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1108: lnklvl - levels of lnk
1109: bt - updated PetscBT (bitarray)
1110: */
1111: #define PetscIncompleteLLInit(nidx,idx,idx_start,perm,nlnk,lnk,lnklvl,bt) 0;\
1112: {\
1113: PetscInt _k,_entry,_location,_lnkdata;\
1114: nlnk = 0;\
1115: _lnkdata = idx_start;\
1116: for (_k=0; _k<nidx; _k++){\
1117: _entry = perm[idx[_k]];\
1118: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1119: /* search for insertion location */\
1120: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
1121: do {\
1122: _location = _lnkdata;\
1123: _lnkdata = lnk[_location];\
1124: } while (_entry > _lnkdata);\
1125: /* insertion location is found, add entry into lnk */\
1126: lnk[_location] = _entry;\
1127: lnk[_entry] = _lnkdata;\
1128: lnklvl[_entry] = 0;\
1129: nlnk++;\
1130: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
1131: }\
1132: }\
1133: }
1135: /*
1136: Add a SORTED index set into a sorted linked list for ILU
1137: Input Parameters:
1138: nidx - number of input indices
1139: idx - sorted integer array used for storing column indices
1140: level - level of fill, e.g., ICC(level)
1141: idxlvl - level of idx
1142: idx_start - starting index of the list
1143: lnk - linked list(an integer array) that is created
1144: lnklvl - levels of lnk
1145: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1146: prow - the row number of idx
1147: output Parameters:
1148: nlnk - number of newly added idx
1149: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1150: lnklvl - levels of lnk
1151: bt - updated PetscBT (bitarray)
1153: Note: the level of factor(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(i,prow)+lvl(prow,j)+1)
1154: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
1155: */
1156: #define PetscILULLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,lnklvl_prow) 0;\
1157: {\
1158: PetscInt _k,_entry,_location,_lnkdata,_incrlev,_lnklvl_prow=lnklvl[prow];\
1159: nlnk = 0;\
1160: _lnkdata = idx_start;\
1161: for (_k=0; _k<nidx; _k++){\
1162: _incrlev = idxlvl[_k] + _lnklvl_prow + 1;\
1163: if (_incrlev > level) continue;\
1164: _entry = idx[_k];\
1165: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1166: /* search for insertion location */\
1167: do {\
1168: _location = _lnkdata;\
1169: _lnkdata = lnk[_location];\
1170: } while (_entry > _lnkdata);\
1171: /* insertion location is found, add entry into lnk */\
1172: lnk[_location] = _entry;\
1173: lnk[_entry] = _lnkdata;\
1174: lnklvl[_entry] = _incrlev;\
1175: nlnk++;\
1176: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
1177: } else { /* existing entry: update lnklvl */\
1178: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1179: }\
1180: }\
1181: }
1183: /*
1184: Add a index set into a sorted linked list
1185: Input Parameters:
1186: nidx - number of input idx
1187: idx - integer array used for storing column indices
1188: level - level of fill, e.g., ICC(level)
1189: idxlvl - level of idx
1190: idx_start - starting index of the list
1191: lnk - linked list(an integer array) that is created
1192: lnklvl - levels of lnk
1193: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1194: output Parameters:
1195: nlnk - number of newly added idx
1196: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1197: lnklvl - levels of lnk
1198: bt - updated PetscBT (bitarray)
1199: */
1200: #define PetscIncompleteLLAdd(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
1201: {\
1202: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1203: nlnk = 0;\
1204: _lnkdata = idx_start;\
1205: for (_k=0; _k<nidx; _k++){\
1206: _incrlev = idxlvl[_k] + 1;\
1207: if (_incrlev > level) continue;\
1208: _entry = idx[_k];\
1209: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1210: /* search for insertion location */\
1211: if (_k && _entry < _lnkdata) _lnkdata = idx_start;\
1212: do {\
1213: _location = _lnkdata;\
1214: _lnkdata = lnk[_location];\
1215: } while (_entry > _lnkdata);\
1216: /* insertion location is found, add entry into lnk */\
1217: lnk[_location] = _entry;\
1218: lnk[_entry] = _lnkdata;\
1219: lnklvl[_entry] = _incrlev;\
1220: nlnk++;\
1221: _lnkdata = _entry; /* next search starts from here if next_entry > _entry */\
1222: } else { /* existing entry: update lnklvl */\
1223: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1224: }\
1225: }\
1226: }
1228: /*
1229: Add a SORTED index set into a sorted linked list
1230: Input Parameters:
1231: nidx - number of input indices
1232: idx - sorted integer array used for storing column indices
1233: level - level of fill, e.g., ICC(level)
1234: idxlvl - level of idx
1235: idx_start - starting index of the list
1236: lnk - linked list(an integer array) that is created
1237: lnklvl - levels of lnk
1238: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1239: output Parameters:
1240: nlnk - number of newly added idx
1241: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1242: lnklvl - levels of lnk
1243: bt - updated PetscBT (bitarray)
1244: */
1245: #define PetscIncompleteLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt) 0;\
1246: {\
1247: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1248: nlnk = 0;\
1249: _lnkdata = idx_start;\
1250: for (_k=0; _k<nidx; _k++){\
1251: _incrlev = idxlvl[_k] + 1;\
1252: if (_incrlev > level) continue;\
1253: _entry = idx[_k];\
1254: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1255: /* search for insertion location */\
1256: do {\
1257: _location = _lnkdata;\
1258: _lnkdata = lnk[_location];\
1259: } while (_entry > _lnkdata);\
1260: /* insertion location is found, add entry into lnk */\
1261: lnk[_location] = _entry;\
1262: lnk[_entry] = _lnkdata;\
1263: lnklvl[_entry] = _incrlev;\
1264: nlnk++;\
1265: _lnkdata = _entry; /* next search starts from here */\
1266: } else { /* existing entry: update lnklvl */\
1267: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1268: }\
1269: }\
1270: }
1272: /*
1273: Add a SORTED index set into a sorted linked list for ICC
1274: Input Parameters:
1275: nidx - number of input indices
1276: idx - sorted integer array used for storing column indices
1277: level - level of fill, e.g., ICC(level)
1278: idxlvl - level of idx
1279: idx_start - starting index of the list
1280: lnk - linked list(an integer array) that is created
1281: lnklvl - levels of lnk
1282: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1283: idxlvl_prow - idxlvl[prow], where prow is the row number of the idx
1284: output Parameters:
1285: nlnk - number of newly added indices
1286: lnk - the sorted(increasing order) linked list containing new and non-redundate entries from idx
1287: lnklvl - levels of lnk
1288: bt - updated PetscBT (bitarray)
1289: Note: the level of U(i,j) is set as lvl(i,j) = min{ lvl(i,j), lvl(prow,i)+lvl(prow,j)+1)
1290: where idx = non-zero columns of U(prow,prow+1:n-1), prow<i
1291: */
1292: #define PetscICCLLAddSorted(nidx,idx,level,idxlvl,idx_start,nlnk,lnk,lnklvl,bt,idxlvl_prow) 0;\
1293: {\
1294: PetscInt _k,_entry,_location,_lnkdata,_incrlev;\
1295: nlnk = 0;\
1296: _lnkdata = idx_start;\
1297: for (_k=0; _k<nidx; _k++){\
1298: _incrlev = idxlvl[_k] + idxlvl_prow + 1;\
1299: if (_incrlev > level) continue;\
1300: _entry = idx[_k];\
1301: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */\
1302: /* search for insertion location */\
1303: do {\
1304: _location = _lnkdata;\
1305: _lnkdata = lnk[_location];\
1306: } while (_entry > _lnkdata);\
1307: /* insertion location is found, add entry into lnk */\
1308: lnk[_location] = _entry;\
1309: lnk[_entry] = _lnkdata;\
1310: lnklvl[_entry] = _incrlev;\
1311: nlnk++;\
1312: _lnkdata = _entry; /* next search starts from here */\
1313: } else { /* existing entry: update lnklvl */\
1314: if (lnklvl[_entry] > _incrlev) lnklvl[_entry] = _incrlev;\
1315: }\
1316: }\
1317: }
1319: /*
1320: Copy data on the list into an array, then initialize the list
1321: Input Parameters:
1322: idx_start - starting index of the list
1323: lnk_max - max value of lnk indicating the end of the list
1324: nlnk - number of data on the list to be copied
1325: lnk - linked list
1326: lnklvl - level of lnk
1327: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1328: output Parameters:
1329: indices - array that contains the copied data
1330: lnk - linked list that is cleaned and initialize
1331: lnklvl - level of lnk that is reinitialized
1332: bt - PetscBT (bitarray) with all bits set to false
1333: */
1334: #define PetscIncompleteLLClean(idx_start,lnk_max,nlnk,lnk,lnklvl,indices,indiceslvl,bt) 0;\
1335: do {\
1336: PetscInt _j,_idx=idx_start;\
1337: for (_j=0; _j<nlnk; _j++){\
1338: _idx = lnk[_idx];\
1339: *(indices+_j) = _idx;\
1340: *(indiceslvl+_j) = lnklvl[_idx];\
1341: lnklvl[_idx] = -1;\
1342: PetscBTClear(bt,_idx);\
1343: }\
1344: lnk[idx_start] = lnk_max;\
1345: } while (0)
1346: /*
1347: Free memories used by the list
1348: */
1349: #define PetscIncompleteLLDestroy(lnk,bt) (PetscFree(lnk) || PetscBTDestroy(&(bt)))
1351: #define MatCheckSameLocalSize(A,ar1,B,ar2) do { \
1353: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n)) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible matrix local sizes: parameter # %d (%D x %D) != parameter # %d (%D x %D)",ar1,A->rmap->n,A->cmap->n,ar2,B->rmap->n,B->cmap->n);} while (0)
1355: #define MatCheckSameSize(A,ar1,B,ar2) do { \
1356: if ((A->rmap->N != B->rmap->N) || (A->cmap->N != B->cmap->N)) SETERRQ6(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible matrix global sizes: parameter # %d (%D x %D) != parameter # %d (%D x %D)",ar1,A->rmap->N,A->cmap->N,ar2,B->rmap->N,B->cmap->N);\
1357: MatCheckSameLocalSize(A,ar1,B,ar2);} while (0)
1359: #define VecCheckMatCompatible(M,x,ar1,b,ar2) do { \
1360: if (M->cmap->N != x->map->N) SETERRQ3(PetscObjectComm((PetscObject)M),PETSC_ERR_ARG_SIZ,"Vector global length incompatible with matrix: parameter # %d global size %D != matrix column global size %D",ar1,x->map->N,M->cmap->N); \
1361: if (M->rmap->N != b->map->N) SETERRQ3(PetscObjectComm((PetscObject)M),PETSC_ERR_ARG_SIZ,"Vector global length incompatible with matrix: parameter # %d global size %D != matrix row global size %D",ar2,b->map->N,M->rmap->N);} while (0)
1363: /* -------------------------------------------------------------------------------------------------------*/
1364: #include <petscbt.h>
1365: /*
1366: Create and initialize a condensed linked list -
1367: same as PetscLLCreate(), but uses a scalable array 'lnk' with size of max number of entries, not O(N).
1368: Barry suggested this approach (Dec. 6, 2011):
1369: I've thought of an alternative way of representing a linked list that is efficient but doesn't have the O(N) scaling issue
1370: (it may be faster than the O(N) even sequentially due to less crazy memory access).
1372: 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
1373: 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
1374: 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
1375: 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.
1376: 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
1377: to each other so memory access is much better than using the big array.
1379: Example:
1380: nlnk_max=5, lnk_max=36:
1381: Initial list: [0, 0 | 36, 2 | 0, 0 | 0, 0 | 0, 0 | 0, 0 | 0, 0]
1382: here, head_node has index 2 with value lnk[2]=lnk_max=36,
1383: 0-th entry is used to store the number of entries in the list,
1384: The initial lnk represents head -> tail(marked by 36) with number of entries = lnk[0]=0.
1386: Now adding a sorted set {2,4}, the list becomes
1387: [2, 0 | 36, 4 |2, 6 | 4, 2 | 0, 0 | 0, 0 | 0, 0 ]
1388: represents head -> 2 -> 4 -> tail with number of entries = lnk[0]=2.
1390: Then adding a sorted set {0,3,35}, the list
1391: [5, 0 | 36, 8 | 2, 10 | 4, 12 | 0, 4 | 3, 6 | 35, 2 ]
1392: represents head -> 0 -> 2 -> 3 -> 4 -> 35 -> tail with number of entries = lnk[0]=5.
1394: Input Parameters:
1395: nlnk_max - max length of the list
1396: lnk_max - max value of the entries
1397: Output Parameters:
1398: lnk - list created and initialized
1399: bt - PetscBT (bitarray) with all bits set to false. Note: bt has size lnk_max, not nln_max!
1400: */
1401: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate(PetscInt nlnk_max,PetscInt lnk_max,PetscInt **lnk,PetscBT *bt)
1402: {
1404: PetscInt *llnk,lsize = 0;
1407: PetscIntMultError(2,nlnk_max+2,&lsize);
1408: PetscMalloc1(lsize,lnk);
1409: PetscBTCreate(lnk_max,bt);
1410: llnk = *lnk;
1411: llnk[0] = 0; /* number of entries on the list */
1412: llnk[2] = lnk_max; /* value in the head node */
1413: llnk[3] = 2; /* next for the head node */
1414: return(0);
1415: }
1417: /*
1418: Add a SORTED ascending index set into a sorted linked list. See PetscLLCondensedCreate() for detailed description.
1419: Input Parameters:
1420: nidx - number of input indices
1421: indices - sorted integer array
1422: lnk - condensed linked list(an integer array) that is created
1423: bt - PetscBT (bitarray), bt[idx]=true marks idx is in lnk
1424: output Parameters:
1425: lnk - the sorted(increasing order) linked list containing previous and newly added non-redundate indices
1426: bt - updated PetscBT (bitarray)
1427: */
1428: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted(PetscInt nidx,const PetscInt indices[],PetscInt lnk[],PetscBT bt)
1429: {
1430: PetscInt _k,_entry,_location,_next,_lnkdata,_nlnk,_newnode;
1433: _nlnk = lnk[0]; /* num of entries on the input lnk */
1434: _location = 2; /* head */
1435: for (_k=0; _k<nidx; _k++){
1436: _entry = indices[_k];
1437: if (!PetscBTLookupSet(bt,_entry)){ /* new entry */
1438: /* search for insertion location */
1439: do {
1440: _next = _location + 1; /* link from previous node to next node */
1441: _location = lnk[_next]; /* idx of next node */
1442: _lnkdata = lnk[_location];/* value of next node */
1443: } while (_entry > _lnkdata);
1444: /* insertion location is found, add entry into lnk */
1445: _newnode = 2*(_nlnk+2); /* index for this new node */
1446: lnk[_next] = _newnode; /* connect previous node to the new node */
1447: lnk[_newnode] = _entry; /* set value of the new node */
1448: lnk[_newnode+1] = _location; /* connect new node to next node */
1449: _location = _newnode; /* next search starts from the new node */
1450: _nlnk++;
1451: } \
1452: }\
1453: lnk[0] = _nlnk; /* number of entries in the list */
1454: return(0);
1455: }
1457: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean(PetscInt lnk_max,PetscInt nidx,PetscInt *indices,PetscInt lnk[],PetscBT bt)
1458: {
1460: PetscInt _k,_next,_nlnk;
1463: _next = lnk[3]; /* head node */
1464: _nlnk = lnk[0]; /* num of entries on the list */
1465: for (_k=0; _k<_nlnk; _k++){
1466: indices[_k] = lnk[_next];
1467: _next = lnk[_next + 1];
1468: PetscBTClear(bt,indices[_k]);
1469: }
1470: lnk[0] = 0; /* num of entries on the list */
1471: lnk[2] = lnk_max; /* initialize head node */
1472: lnk[3] = 2; /* head node */
1473: return(0);
1474: }
1476: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedView(PetscInt *lnk)
1477: {
1479: PetscInt k;
1482: PetscPrintf(PETSC_COMM_SELF,"LLCondensed of size %D, (val, next)\n",lnk[0]);
1483: for (k=2; k< lnk[0]+2; k++){
1484: PetscPrintf(PETSC_COMM_SELF," %D: (%D, %D)\n",2*k,lnk[2*k],lnk[2*k+1]);
1485: }
1486: return(0);
1487: }
1489: /*
1490: Free memories used by the list
1491: */
1492: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy(PetscInt *lnk,PetscBT bt)
1493: {
1497: PetscFree(lnk);
1498: PetscBTDestroy(&bt);
1499: return(0);
1500: }
1502: /* -------------------------------------------------------------------------------------------------------*/
1503: /*
1504: Same as PetscLLCondensedCreate(), but does not use non-scalable O(lnk_max) bitarray
1505: Input Parameters:
1506: nlnk_max - max length of the list
1507: Output Parameters:
1508: lnk - list created and initialized
1509: */
1510: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate_Scalable(PetscInt nlnk_max,PetscInt **lnk)
1511: {
1513: PetscInt *llnk,lsize = 0;
1516: PetscIntMultError(2,nlnk_max+2,&lsize);
1517: PetscMalloc1(lsize,lnk);
1518: llnk = *lnk;
1519: llnk[0] = 0; /* number of entries on the list */
1520: llnk[2] = PETSC_MAX_INT; /* value in the head node */
1521: llnk[3] = 2; /* next for the head node */
1522: return(0);
1523: }
1525: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedExpand_Scalable(PetscInt nlnk_max,PetscInt **lnk)
1526: {
1528: PetscInt lsize = 0;
1531: PetscIntMultError(2,nlnk_max+2,&lsize);
1532: PetscRealloc(lsize*sizeof(PetscInt),lnk);
1533: return(0);
1534: }
1536: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted_Scalable(PetscInt nidx,const PetscInt indices[],PetscInt lnk[])
1537: {
1538: PetscInt _k,_entry,_location,_next,_lnkdata,_nlnk,_newnode;
1539: _nlnk = lnk[0]; /* num of entries on the input lnk */
1540: _location = 2; /* head */ \
1541: for (_k=0; _k<nidx; _k++){
1542: _entry = indices[_k];
1543: /* search for insertion location */
1544: do {
1545: _next = _location + 1; /* link from previous node to next node */
1546: _location = lnk[_next]; /* idx of next node */
1547: _lnkdata = lnk[_location];/* value of next node */
1548: } while (_entry > _lnkdata);
1549: if (_entry < _lnkdata) {
1550: /* insertion location is found, add entry into lnk */
1551: _newnode = 2*(_nlnk+2); /* index for this new node */
1552: lnk[_next] = _newnode; /* connect previous node to the new node */
1553: lnk[_newnode] = _entry; /* set value of the new node */
1554: lnk[_newnode+1] = _location; /* connect new node to next node */
1555: _location = _newnode; /* next search starts from the new node */
1556: _nlnk++;
1557: }
1558: }
1559: lnk[0] = _nlnk; /* number of entries in the list */
1560: return 0;
1561: }
1563: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean_Scalable(PetscInt nidx,PetscInt *indices,PetscInt *lnk)
1564: {
1565: PetscInt _k,_next,_nlnk;
1566: _next = lnk[3]; /* head node */
1567: _nlnk = lnk[0];
1568: for (_k=0; _k<_nlnk; _k++){
1569: indices[_k] = lnk[_next];
1570: _next = lnk[_next + 1];
1571: }
1572: lnk[0] = 0; /* num of entries on the list */
1573: lnk[3] = 2; /* head node */
1574: return 0;
1575: }
1577: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy_Scalable(PetscInt *lnk)
1578: {
1579: return PetscFree(lnk);
1580: }
1582: /* -------------------------------------------------------------------------------------------------------*/
1583: /*
1584: lnk[0] number of links
1585: lnk[1] number of entries
1586: lnk[3n] value
1587: lnk[3n+1] len
1588: lnk[3n+2] link to next value
1590: The next three are always the first link
1592: lnk[3] PETSC_MIN_INT+1
1593: lnk[4] 1
1594: lnk[5] link to first real entry
1596: The next three are always the last link
1598: lnk[6] PETSC_MAX_INT - 1
1599: lnk[7] 1
1600: lnk[8] next valid link (this is the same as lnk[0] but without the decreases)
1601: */
1603: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedCreate_fast(PetscInt nlnk_max,PetscInt **lnk)
1604: {
1606: PetscInt *llnk,lsize = 0;
1609: PetscIntMultError(3,nlnk_max+3,&lsize);
1610: PetscMalloc1(lsize,lnk);
1611: llnk = *lnk;
1612: llnk[0] = 0; /* nlnk: number of entries on the list */
1613: llnk[1] = 0; /* number of integer entries represented in list */
1614: llnk[3] = PETSC_MIN_INT+1; /* value in the first node */
1615: llnk[4] = 1; /* count for the first node */
1616: llnk[5] = 6; /* next for the first node */
1617: llnk[6] = PETSC_MAX_INT-1; /* value in the last node */
1618: llnk[7] = 1; /* count for the last node */
1619: llnk[8] = 0; /* next valid node to be used */
1620: return(0);
1621: }
1623: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedAddSorted_fast(PetscInt nidx,const PetscInt indices[],PetscInt lnk[])
1624: {
1625: PetscInt k,entry,prev,next;
1626: prev = 3; /* first value */
1627: next = lnk[prev+2];
1628: for (k=0; k<nidx; k++){
1629: entry = indices[k];
1630: /* search for insertion location */
1631: while (entry >= lnk[next]) {
1632: prev = next;
1633: next = lnk[next+2];
1634: }
1635: /* entry is in range of previous list */
1636: if (entry < lnk[prev]+lnk[prev+1]) continue;
1637: lnk[1]++;
1638: /* entry is right after previous list */
1639: if (entry == lnk[prev]+lnk[prev+1]) {
1640: lnk[prev+1]++;
1641: if (lnk[next] == entry+1) { /* combine two contiguous strings */
1642: lnk[prev+1] += lnk[next+1];
1643: lnk[prev+2] = lnk[next+2];
1644: next = lnk[next+2];
1645: lnk[0]--;
1646: }
1647: continue;
1648: }
1649: /* entry is right before next list */
1650: if (entry == lnk[next]-1) {
1651: lnk[next]--;
1652: lnk[next+1]++;
1653: prev = next;
1654: next = lnk[prev+2];
1655: continue;
1656: }
1657: /* add entry into lnk */
1658: lnk[prev+2] = 3*((lnk[8]++)+3); /* connect previous node to the new node */
1659: prev = lnk[prev+2];
1660: lnk[prev] = entry; /* set value of the new node */
1661: lnk[prev+1] = 1; /* number of values in contiguous string is one to start */
1662: lnk[prev+2] = next; /* connect new node to next node */
1663: lnk[0]++;
1664: }
1665: return 0;
1666: }
1668: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedClean_fast(PetscInt nidx,PetscInt *indices,PetscInt *lnk)
1669: {
1670: PetscInt _k,_next,_nlnk,cnt,j;
1671: _next = lnk[5]; /* first node */
1672: _nlnk = lnk[0];
1673: cnt = 0;
1674: for (_k=0; _k<_nlnk; _k++){
1675: for (j=0; j<lnk[_next+1]; j++) {
1676: indices[cnt++] = lnk[_next] + j;
1677: }
1678: _next = lnk[_next + 2];
1679: }
1680: lnk[0] = 0; /* nlnk: number of links */
1681: lnk[1] = 0; /* number of integer entries represented in list */
1682: lnk[3] = PETSC_MIN_INT+1; /* value in the first node */
1683: lnk[4] = 1; /* count for the first node */
1684: lnk[5] = 6; /* next for the first node */
1685: lnk[6] = PETSC_MAX_INT-1; /* value in the last node */
1686: lnk[7] = 1; /* count for the last node */
1687: lnk[8] = 0; /* next valid location to make link */
1688: return 0;
1689: }
1691: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedView_fast(PetscInt *lnk)
1692: {
1693: PetscInt k,next,nlnk;
1694: next = lnk[5]; /* first node */
1695: nlnk = lnk[0];
1696: for (k=0; k<nlnk; k++){
1697: #if 0 /* Debugging code */
1698: printf("%d value %d len %d next %d\n",next,lnk[next],lnk[next+1],lnk[next+2]);
1699: #endif
1700: next = lnk[next + 2];
1701: }
1702: return 0;
1703: }
1705: PETSC_STATIC_INLINE PetscErrorCode PetscLLCondensedDestroy_fast(PetscInt *lnk)
1706: {
1707: return PetscFree(lnk);
1708: }
1710: /* this is extern because it is used in MatFDColoringUseDM() which is in the DM library */
1711: PETSC_EXTERN PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,void*);
1713: PETSC_EXTERN PetscLogEvent MAT_Mult;
1714: PETSC_EXTERN PetscLogEvent MAT_MultMatrixFree;
1715: PETSC_EXTERN PetscLogEvent MAT_Mults;
1716: PETSC_EXTERN PetscLogEvent MAT_MultConstrained;
1717: PETSC_EXTERN PetscLogEvent MAT_MultAdd;
1718: PETSC_EXTERN PetscLogEvent MAT_MultTranspose;
1719: PETSC_EXTERN PetscLogEvent MAT_MultTransposeConstrained;
1720: PETSC_EXTERN PetscLogEvent MAT_MultTransposeAdd;
1721: PETSC_EXTERN PetscLogEvent MAT_Solve;
1722: PETSC_EXTERN PetscLogEvent MAT_Solves;
1723: PETSC_EXTERN PetscLogEvent MAT_SolveAdd;
1724: PETSC_EXTERN PetscLogEvent MAT_SolveTranspose;
1725: PETSC_EXTERN PetscLogEvent MAT_SolveTransposeAdd;
1726: PETSC_EXTERN PetscLogEvent MAT_SOR;
1727: PETSC_EXTERN PetscLogEvent MAT_ForwardSolve;
1728: PETSC_EXTERN PetscLogEvent MAT_BackwardSolve;
1729: PETSC_EXTERN PetscLogEvent MAT_LUFactor;
1730: PETSC_EXTERN PetscLogEvent MAT_LUFactorSymbolic;
1731: PETSC_EXTERN PetscLogEvent MAT_LUFactorNumeric;
1732: PETSC_EXTERN PetscLogEvent MAT_CholeskyFactor;
1733: PETSC_EXTERN PetscLogEvent MAT_CholeskyFactorSymbolic;
1734: PETSC_EXTERN PetscLogEvent MAT_CholeskyFactorNumeric;
1735: PETSC_EXTERN PetscLogEvent MAT_ILUFactor;
1736: PETSC_EXTERN PetscLogEvent MAT_ILUFactorSymbolic;
1737: PETSC_EXTERN PetscLogEvent MAT_ICCFactorSymbolic;
1738: PETSC_EXTERN PetscLogEvent MAT_Copy;
1739: PETSC_EXTERN PetscLogEvent MAT_Convert;
1740: PETSC_EXTERN PetscLogEvent MAT_Scale;
1741: PETSC_EXTERN PetscLogEvent MAT_AssemblyBegin;
1742: PETSC_EXTERN PetscLogEvent MAT_AssemblyEnd;
1743: PETSC_EXTERN PetscLogEvent MAT_SetValues;
1744: PETSC_EXTERN PetscLogEvent MAT_GetValues;
1745: PETSC_EXTERN PetscLogEvent MAT_GetRow;
1746: PETSC_EXTERN PetscLogEvent MAT_GetRowIJ;
1747: PETSC_EXTERN PetscLogEvent MAT_CreateSubMats;
1748: PETSC_EXTERN PetscLogEvent MAT_GetColoring;
1749: PETSC_EXTERN PetscLogEvent MAT_GetOrdering;
1750: PETSC_EXTERN PetscLogEvent MAT_RedundantMat;
1751: PETSC_EXTERN PetscLogEvent MAT_IncreaseOverlap;
1752: PETSC_EXTERN PetscLogEvent MAT_Partitioning;
1753: PETSC_EXTERN PetscLogEvent MAT_PartitioningND;
1754: PETSC_EXTERN PetscLogEvent MAT_Coarsen;
1755: PETSC_EXTERN PetscLogEvent MAT_ZeroEntries;
1756: PETSC_EXTERN PetscLogEvent MAT_Load;
1757: PETSC_EXTERN PetscLogEvent MAT_View;
1758: PETSC_EXTERN PetscLogEvent MAT_AXPY;
1759: PETSC_EXTERN PetscLogEvent MAT_FDColoringCreate;
1760: PETSC_EXTERN PetscLogEvent MAT_TransposeColoringCreate;
1761: PETSC_EXTERN PetscLogEvent MAT_FDColoringSetUp;
1762: PETSC_EXTERN PetscLogEvent MAT_FDColoringApply;
1763: PETSC_EXTERN PetscLogEvent MAT_Transpose;
1764: PETSC_EXTERN PetscLogEvent MAT_FDColoringFunction;
1765: PETSC_EXTERN PetscLogEvent MAT_CreateSubMat;
1766: PETSC_EXTERN PetscLogEvent MAT_MatSolve;
1767: PETSC_EXTERN PetscLogEvent MAT_MatTrSolve;
1768: PETSC_EXTERN PetscLogEvent MAT_MatMultSymbolic;
1769: PETSC_EXTERN PetscLogEvent MAT_MatMultNumeric;
1770: PETSC_EXTERN PetscLogEvent MAT_Getlocalmatcondensed;
1771: PETSC_EXTERN PetscLogEvent MAT_GetBrowsOfAcols;
1772: PETSC_EXTERN PetscLogEvent MAT_GetBrowsOfAocols;
1773: PETSC_EXTERN PetscLogEvent MAT_PtAPSymbolic;
1774: PETSC_EXTERN PetscLogEvent MAT_PtAPNumeric;
1775: PETSC_EXTERN PetscLogEvent MAT_Seqstompinum;
1776: PETSC_EXTERN PetscLogEvent MAT_Seqstompisym;
1777: PETSC_EXTERN PetscLogEvent MAT_Seqstompi;
1778: PETSC_EXTERN PetscLogEvent MAT_Getlocalmat;
1779: PETSC_EXTERN PetscLogEvent MAT_RARtSymbolic;
1780: PETSC_EXTERN PetscLogEvent MAT_RARtNumeric;
1781: PETSC_EXTERN PetscLogEvent MAT_MatTransposeMultSymbolic;
1782: PETSC_EXTERN PetscLogEvent MAT_MatTransposeMultNumeric;
1783: PETSC_EXTERN PetscLogEvent MAT_TransposeMatMultSymbolic;
1784: PETSC_EXTERN PetscLogEvent MAT_TransposeMatMultNumeric;
1785: PETSC_EXTERN PetscLogEvent MAT_MatMatMultSymbolic;
1786: PETSC_EXTERN PetscLogEvent MAT_MatMatMultNumeric;
1787: PETSC_EXTERN PetscLogEvent MAT_Applypapt;
1788: PETSC_EXTERN PetscLogEvent MAT_Applypapt_symbolic;
1789: PETSC_EXTERN PetscLogEvent MAT_Applypapt_numeric;
1790: PETSC_EXTERN PetscLogEvent MAT_Getsymtranspose;
1791: PETSC_EXTERN PetscLogEvent MAT_Getsymtransreduced;
1792: PETSC_EXTERN PetscLogEvent MAT_GetSequentialNonzeroStructure;
1793: PETSC_EXTERN PetscLogEvent MATMFFD_Mult;
1794: PETSC_EXTERN PetscLogEvent MAT_GetMultiProcBlock;
1795: PETSC_EXTERN PetscLogEvent MAT_CUSPARSECopyToGPU;
1796: PETSC_EXTERN PetscLogEvent MAT_CUSPARSECopyFromGPU;
1797: PETSC_EXTERN PetscLogEvent MAT_CUSPARSEGenerateTranspose;
1798: PETSC_EXTERN PetscLogEvent MAT_CUSPARSESolveAnalysis;
1799: PETSC_EXTERN PetscLogEvent MAT_SetValuesBatch;
1800: PETSC_EXTERN PetscLogEvent MAT_ViennaCLCopyToGPU;
1801: PETSC_EXTERN PetscLogEvent MAT_DenseCopyToGPU;
1802: PETSC_EXTERN PetscLogEvent MAT_DenseCopyFromGPU;
1803: PETSC_EXTERN PetscLogEvent MAT_Merge;
1804: PETSC_EXTERN PetscLogEvent MAT_Residual;
1805: PETSC_EXTERN PetscLogEvent MAT_SetRandom;
1806: PETSC_EXTERN PetscLogEvent MAT_FactorFactS;
1807: PETSC_EXTERN PetscLogEvent MAT_FactorInvS;
1808: PETSC_EXTERN PetscLogEvent MAT_PreallCOO;
1809: PETSC_EXTERN PetscLogEvent MAT_SetVCOO;
1810: PETSC_EXTERN PetscLogEvent MATCOLORING_Apply;
1811: PETSC_EXTERN PetscLogEvent MATCOLORING_Comm;
1812: PETSC_EXTERN PetscLogEvent MATCOLORING_Local;
1813: PETSC_EXTERN PetscLogEvent MATCOLORING_ISCreate;
1814: PETSC_EXTERN PetscLogEvent MATCOLORING_SetUp;
1815: PETSC_EXTERN PetscLogEvent MATCOLORING_Weights;
1817: #endif