Actual source code: aijsell.c
petsc-3.12.5 2020-03-29
1: /*
2: Defines basic operations for the MATSEQAIJSELL matrix class.
3: This class is derived from the MATAIJCLASS, but maintains a "shadow" copy
4: of the matrix stored in MATSEQSELL format, which is used as appropriate for
5: performing operations for which this format is more suitable.
6: */
8: #include <../src/mat/impls/aij/seq/aij.h>
9: #include <../src/mat/impls/sell/seq/sell.h>
11: typedef struct {
12: Mat S; /* The SELL formatted "shadow" matrix. */
13: PetscBool eager_shadow;
14: PetscObjectState state; /* State of the matrix when shadow matrix was last constructed. */
15: } Mat_SeqAIJSELL;
17: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJSELL_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
18: {
19: /* This routine is only called to convert a MATAIJSELL to its base PETSc type, */
20: /* so we will ignore 'MatType type'. */
22: Mat B = *newmat;
23: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*) A->spptr;
26: if (reuse == MAT_INITIAL_MATRIX) {
27: MatDuplicate(A,MAT_COPY_VALUES,&B);
28: }
30: /* Reset the original function pointers. */
31: B->ops->duplicate = MatDuplicate_SeqAIJ;
32: B->ops->assemblyend = MatAssemblyEnd_SeqAIJ;
33: B->ops->destroy = MatDestroy_SeqAIJ;
34: B->ops->mult = MatMult_SeqAIJ;
35: B->ops->multtranspose = MatMultTranspose_SeqAIJ;
36: B->ops->multadd = MatMultAdd_SeqAIJ;
37: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJ;
38: B->ops->sor = MatSOR_SeqAIJ;
40: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijsell_seqaij_C",NULL);
41: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaijsell_C",NULL);
42: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaijsell_C",NULL);
43: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",NULL);
44: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaijsell_C",NULL);
46: if (reuse == MAT_INITIAL_MATRIX) aijsell = (Mat_SeqAIJSELL*)B->spptr;
48: /* Clean up the Mat_SeqAIJSELL data structure.
49: * Note that MatDestroy() simply returns if passed a NULL value, so it's OK to call even if the shadow matrix was never constructed. */
50: MatDestroy(&aijsell->S);
51: PetscFree(B->spptr);
53: /* Change the type of B to MATSEQAIJ. */
54: PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ);
56: *newmat = B;
57: return(0);
58: }
60: PetscErrorCode MatDestroy_SeqAIJSELL(Mat A)
61: {
63: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*) A->spptr;
67: /* If MatHeaderMerge() was used, then this SeqAIJSELL matrix will not have an
68: * spptr pointer. */
69: if (aijsell) {
70: /* Clean up everything in the Mat_SeqAIJSELL data structure, then free A->spptr. */
71: MatDestroy(&aijsell->S);
72: PetscFree(A->spptr);
73: }
75: /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ()
76: * to destroy everything that remains. */
77: PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ);
78: /* Note that I don't call MatSetType(). I believe this is because that
79: * is only to be called when *building* a matrix. I could be wrong, but
80: * that is how things work for the SuperLU matrix class. */
81: MatDestroy_SeqAIJ(A);
82: return(0);
83: }
85: /* Build or update the shadow matrix if and only if needed.
86: * We track the ObjectState to determine when this needs to be done. */
87: PETSC_INTERN PetscErrorCode MatSeqAIJSELL_build_shadow(Mat A)
88: {
89: PetscErrorCode ierr;
90: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*) A->spptr;
91: PetscObjectState state;
94: PetscObjectStateGet((PetscObject)A,&state);
95: if (aijsell->S && aijsell->state == state) {
96: /* The existing shadow matrix is up-to-date, so simply exit. */
97: return(0);
98: }
100: PetscLogEventBegin(MAT_Convert,A,0,0,0);
101: if (aijsell->S) {
102: MatConvert_SeqAIJ_SeqSELL(A,MATSEQSELL,MAT_REUSE_MATRIX,&aijsell->S);
103: } else {
104: MatConvert_SeqAIJ_SeqSELL(A,MATSEQSELL,MAT_INITIAL_MATRIX,&aijsell->S);
105: }
106: PetscLogEventEnd(MAT_Convert,A,0,0,0);
108: /* Record the ObjectState so that we can tell when the shadow matrix needs updating */
109: PetscObjectStateGet((PetscObject)A,&aijsell->state);
111: return(0);
112: }
114: PetscErrorCode MatDuplicate_SeqAIJSELL(Mat A, MatDuplicateOption op, Mat *M)
115: {
117: Mat_SeqAIJSELL *aijsell;
118: Mat_SeqAIJSELL *aijsell_dest;
121: MatDuplicate_SeqAIJ(A,op,M);
122: aijsell = (Mat_SeqAIJSELL*) A->spptr;
123: aijsell_dest = (Mat_SeqAIJSELL*) (*M)->spptr;
124: PetscArraycpy(aijsell_dest,aijsell,1);
125: /* We don't duplicate the shadow matrix -- that will be constructed as needed. */
126: aijsell_dest->S = NULL;
127: if (aijsell->eager_shadow) {
128: MatSeqAIJSELL_build_shadow(A);
129: }
130: return(0);
131: }
133: PetscErrorCode MatAssemblyEnd_SeqAIJSELL(Mat A, MatAssemblyType mode)
134: {
135: PetscErrorCode ierr;
136: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
137: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*)A->spptr;
140: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
142: /* I disable the use of the inode routines so that the AIJSELL ones will be
143: * used instead, but I wonder if it might make sense (and is feasible) to
144: * use some of them. */
145: a->inode.use = PETSC_FALSE;
147: /* Since a MATSEQAIJSELL matrix is really just a MATSEQAIJ with some
148: * extra information and some different methods, call the AssemblyEnd
149: * routine for a MATSEQAIJ.
150: * I'm not sure if this is the best way to do this, but it avoids
151: * a lot of code duplication. */
153: MatAssemblyEnd_SeqAIJ(A, mode);
155: /* If the user has requested "eager" shadowing, create the SELL shadow matrix (if needed; the function checks).
156: * (The default is to take a "lazy" approach, deferring this until something like MatMult() is called.) */
157: if (aijsell->eager_shadow) {
158: MatSeqAIJSELL_build_shadow(A);
159: }
161: return(0);
162: }
164: PetscErrorCode MatMult_SeqAIJSELL(Mat A,Vec xx,Vec yy)
165: {
166: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*)A->spptr;
167: PetscErrorCode ierr;
170: MatSeqAIJSELL_build_shadow(A);
171: MatMult_SeqSELL(aijsell->S,xx,yy);
172: return(0);
173: }
175: PetscErrorCode MatMultTranspose_SeqAIJSELL(Mat A,Vec xx,Vec yy)
176: {
177: Mat_SeqAIJSELL *aijsell=(Mat_SeqAIJSELL*)A->spptr;
178: PetscErrorCode ierr;
181: MatSeqAIJSELL_build_shadow(A);
182: MatMultTranspose_SeqSELL(aijsell->S,xx,yy);
183: return(0);
184: }
186: PetscErrorCode MatMultAdd_SeqAIJSELL(Mat A,Vec xx,Vec yy,Vec zz)
187: {
188: Mat_SeqAIJSELL *aijsell=(Mat_SeqAIJSELL*)A->spptr;
189: PetscErrorCode ierr;
192: MatSeqAIJSELL_build_shadow(A);
193: MatMultAdd_SeqSELL(aijsell->S,xx,yy,zz);
194: return(0);
195: }
197: PetscErrorCode MatMultTransposeAdd_SeqAIJSELL(Mat A,Vec xx,Vec yy,Vec zz)
198: {
199: Mat_SeqAIJSELL *aijsell=(Mat_SeqAIJSELL*)A->spptr;
200: PetscErrorCode ierr;
203: MatSeqAIJSELL_build_shadow(A);
204: MatMultTransposeAdd_SeqSELL(aijsell->S,xx,yy,zz);
205: return(0);
206: }
208: PetscErrorCode MatSOR_SeqAIJSELL(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
209: {
210: Mat_SeqAIJSELL *aijsell=(Mat_SeqAIJSELL*)A->spptr;
211: PetscErrorCode ierr;
214: MatSeqAIJSELL_build_shadow(A);
215: MatSOR_SeqSELL(aijsell->S,bb,omega,flag,fshift,its,lits,xx);
216: return(0);
217: }
219: /* This function prototype is needed in MatConvert_SeqAIJ_SeqAIJSELL(), below. */
220: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*);
222: /* MatConvert_SeqAIJ_SeqAIJSELL converts a SeqAIJ matrix into a
223: * SeqAIJSELL matrix. This routine is called by the MatCreate_SeqAIJSELL()
224: * routine, but can also be used to convert an assembled SeqAIJ matrix
225: * into a SeqAIJSELL one. */
226: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJSELL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
227: {
229: Mat B = *newmat;
230: Mat_SeqAIJ *b;
231: Mat_SeqAIJSELL *aijsell;
232: PetscBool set;
233: PetscBool sametype;
236: if (reuse == MAT_INITIAL_MATRIX) {
237: MatDuplicate(A,MAT_COPY_VALUES,&B);
238: }
240: PetscObjectTypeCompare((PetscObject)A,type,&sametype);
241: if (sametype) return(0);
243: PetscNewLog(B,&aijsell);
244: b = (Mat_SeqAIJ*) B->data;
245: B->spptr = (void*) aijsell;
247: /* Disable use of the inode routines so that the AIJSELL ones will be used instead.
248: * This happens in MatAssemblyEnd_SeqAIJSELL as well, but the assembly end may not be called, so set it here, too.
249: * As noted elsewhere, I wonder if it might make sense and be feasible to use some of the inode routines. */
250: b->inode.use = PETSC_FALSE;
252: /* Set function pointers for methods that we inherit from AIJ but override.
253: * We also parse some command line options below, since those determine some of the methods we point to. */
254: B->ops->duplicate = MatDuplicate_SeqAIJSELL;
255: B->ops->assemblyend = MatAssemblyEnd_SeqAIJSELL;
256: B->ops->destroy = MatDestroy_SeqAIJSELL;
258: aijsell->S = NULL;
259: aijsell->eager_shadow = PETSC_FALSE;
261: /* Parse command line options. */
262: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"AIJSELL Options","Mat");
263: PetscOptionsBool("-mat_aijsell_eager_shadow","Eager Shadowing","None",(PetscBool)aijsell->eager_shadow,(PetscBool*)&aijsell->eager_shadow,&set);
264: PetscOptionsEnd();
266: /* If A has already been assembled and eager shadowing is specified, build the shadow matrix. */
267: if (A->assembled && aijsell->eager_shadow) {
268: MatSeqAIJSELL_build_shadow(A);
269: }
271: B->ops->mult = MatMult_SeqAIJSELL;
272: B->ops->multtranspose = MatMultTranspose_SeqAIJSELL;
273: B->ops->multadd = MatMultAdd_SeqAIJSELL;
274: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJSELL;
275: B->ops->sor = MatSOR_SeqAIJSELL;
277: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijsell_seqaij_C",MatConvert_SeqAIJSELL_SeqAIJ);
278: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaijsell_C",MatMatMult_SeqDense_SeqAIJ);
279: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaijsell_C",MatMatMultSymbolic_SeqDense_SeqAIJ);
280: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",MatMatMultNumeric_SeqDense_SeqAIJ);
281: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaijsell_C",MatPtAP_IS_XAIJ);
283: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJSELL);
284: *newmat = B;
285: return(0);
286: }
288: /*@C
289: MatCreateSeqAIJSELL - Creates a sparse matrix of type SEQAIJSELL.
290: This type inherits from AIJ and is largely identical, but keeps a "shadow"
291: copy of the matrix in SEQSELL format, which is used when this format
292: may be more suitable for a requested operation. Currently, SEQSELL format
293: is used for MatMult, MatMultTranspose, MatMultAdd, MatMultTransposeAdd,
294: and MatSOR operations.
295: Because SEQAIJSELL is a subtype of SEQAIJ, the option "-mat_seqaij_type seqaijsell" can be used to make
296: sequential AIJ matrices default to being instances of MATSEQAIJSELL.
298: Collective
300: Input Parameters:
301: + comm - MPI communicator, set to PETSC_COMM_SELF
302: . m - number of rows
303: . n - number of columns
304: . nz - number of nonzeros per row (same for all rows)
305: - nnz - array containing the number of nonzeros in the various rows
306: (possibly different for each row) or NULL
308: Output Parameter:
309: . A - the matrix
311: Options Database Keys:
312: . -mat_aijsell_eager_shadow - Construct shadow matrix upon matrix assembly; default is to take a "lazy" approach, performing this step the first time the matrix is applied
314: Notes:
315: If nnz is given then nz is ignored
317: Level: intermediate
319: .seealso: MatCreate(), MatCreateMPIAIJSELL(), MatSetValues()
320: @*/
321: PetscErrorCode MatCreateSeqAIJSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
322: {
326: MatCreate(comm,A);
327: MatSetSizes(*A,m,n,m,n);
328: MatSetType(*A,MATSEQAIJSELL);
329: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
330: return(0);
331: }
333: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJSELL(Mat A)
334: {
338: MatSetType(A,MATSEQAIJ);
339: MatConvert_SeqAIJ_SeqAIJSELL(A,MATSEQAIJSELL,MAT_INPLACE_MATRIX,&A);
340: return(0);
341: }