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