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