Actual source code: aijsell.c
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'. */
21: Mat B = *newmat;
22: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
24: PetscFunctionBegin;
25: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
27: /* Reset the original function pointers. */
28: B->ops->duplicate = MatDuplicate_SeqAIJ;
29: B->ops->assemblyend = MatAssemblyEnd_SeqAIJ;
30: B->ops->destroy = MatDestroy_SeqAIJ;
31: B->ops->mult = MatMult_SeqAIJ;
32: B->ops->multtranspose = MatMultTranspose_SeqAIJ;
33: B->ops->multadd = MatMultAdd_SeqAIJ;
34: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJ;
35: B->ops->sor = MatSOR_SeqAIJ;
37: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaijsell_seqaij_C", NULL));
39: if (reuse == MAT_INITIAL_MATRIX) aijsell = (Mat_SeqAIJSELL *)B->spptr;
41: /* Clean up the Mat_SeqAIJSELL data structure.
42: * Note that MatDestroy() simply returns if passed a NULL value, so it's OK to call even if the shadow matrix was never constructed. */
43: PetscCall(MatDestroy(&aijsell->S));
44: PetscCall(PetscFree(B->spptr));
46: /* Change the type of B to MATSEQAIJ. */
47: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ));
49: *newmat = B;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode MatDestroy_SeqAIJSELL(Mat A)
54: {
55: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
57: PetscFunctionBegin;
59: /* If MatHeaderMerge() was used, then this SeqAIJSELL matrix will not have an
60: * spptr pointer. */
61: if (aijsell) {
62: /* Clean up everything in the Mat_SeqAIJSELL data structure, then free A->spptr. */
63: PetscCall(MatDestroy(&aijsell->S));
64: PetscCall(PetscFree(A->spptr));
65: }
67: /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ()
68: * to destroy everything that remains. */
69: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ));
70: /* Note that I don't call MatSetType(). I believe this is because that
71: * is only to be called when *building* a matrix. I could be wrong, but
72: * that is how things work for the SuperLU matrix class. */
73: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaijsell_seqaij_C", NULL));
74: PetscCall(MatDestroy_SeqAIJ(A));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /* Build or update the shadow matrix if and only if needed.
79: * We track the ObjectState to determine when this needs to be done. */
80: PETSC_INTERN PetscErrorCode MatSeqAIJSELL_build_shadow(Mat A)
81: {
82: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
83: PetscObjectState state;
85: PetscFunctionBegin;
86: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
87: if (aijsell->S && aijsell->state == state) {
88: /* The existing shadow matrix is up-to-date, so simply exit. */
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: PetscCall(PetscLogEventBegin(MAT_Convert, A, 0, 0, 0));
93: if (aijsell->S) {
94: PetscCall(MatConvert_SeqAIJ_SeqSELL(A, MATSEQSELL, MAT_REUSE_MATRIX, &aijsell->S));
95: } else {
96: PetscCall(MatConvert_SeqAIJ_SeqSELL(A, MATSEQSELL, MAT_INITIAL_MATRIX, &aijsell->S));
97: }
98: PetscCall(PetscLogEventEnd(MAT_Convert, A, 0, 0, 0));
100: /* Record the ObjectState so that we can tell when the shadow matrix needs updating */
101: PetscCall(PetscObjectStateGet((PetscObject)A, &aijsell->state));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode MatDuplicate_SeqAIJSELL(Mat A, MatDuplicateOption op, Mat *M)
107: {
108: Mat_SeqAIJSELL *aijsell;
109: Mat_SeqAIJSELL *aijsell_dest;
111: PetscFunctionBegin;
112: PetscCall(MatDuplicate_SeqAIJ(A, op, M));
113: aijsell = (Mat_SeqAIJSELL *)A->spptr;
114: aijsell_dest = (Mat_SeqAIJSELL *)(*M)->spptr;
115: PetscCall(PetscArraycpy(aijsell_dest, aijsell, 1));
116: /* We don't duplicate the shadow matrix -- that will be constructed as needed. */
117: aijsell_dest->S = NULL;
118: if (aijsell->eager_shadow) PetscCall(MatSeqAIJSELL_build_shadow(A));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode MatAssemblyEnd_SeqAIJSELL(Mat A, MatAssemblyType mode)
123: {
124: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
125: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
127: PetscFunctionBegin;
128: if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(PETSC_SUCCESS);
130: /* I disable the use of the inode routines so that the AIJSELL ones will be
131: * used instead, but I wonder if it might make sense (and is feasible) to
132: * use some of them. */
133: a->inode.use = PETSC_FALSE;
135: /* Since a MATSEQAIJSELL matrix is really just a MATSEQAIJ with some
136: * extra information and some different methods, call the AssemblyEnd
137: * routine for a MATSEQAIJ.
138: * I'm not sure if this is the best way to do this, but it avoids
139: * a lot of code duplication. */
141: PetscCall(MatAssemblyEnd_SeqAIJ(A, mode));
143: /* If the user has requested "eager" shadowing, create the SELL shadow matrix (if needed; the function checks).
144: * (The default is to take a "lazy" approach, deferring this until something like MatMult() is called.) */
145: if (aijsell->eager_shadow) PetscCall(MatSeqAIJSELL_build_shadow(A));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatMult_SeqAIJSELL(Mat A, Vec xx, Vec yy)
151: {
152: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
154: PetscFunctionBegin;
155: PetscCall(MatSeqAIJSELL_build_shadow(A));
156: PetscCall(MatMult_SeqSELL(aijsell->S, xx, yy));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode MatMultTranspose_SeqAIJSELL(Mat A, Vec xx, Vec yy)
161: {
162: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
164: PetscFunctionBegin;
165: PetscCall(MatSeqAIJSELL_build_shadow(A));
166: PetscCall(MatMultTranspose_SeqSELL(aijsell->S, xx, yy));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode MatMultAdd_SeqAIJSELL(Mat A, Vec xx, Vec yy, Vec zz)
171: {
172: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
174: PetscFunctionBegin;
175: PetscCall(MatSeqAIJSELL_build_shadow(A));
176: PetscCall(MatMultAdd_SeqSELL(aijsell->S, xx, yy, zz));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode MatMultTransposeAdd_SeqAIJSELL(Mat A, Vec xx, Vec yy, Vec zz)
181: {
182: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
184: PetscFunctionBegin;
185: PetscCall(MatSeqAIJSELL_build_shadow(A));
186: PetscCall(MatMultTransposeAdd_SeqSELL(aijsell->S, xx, yy, zz));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: static PetscErrorCode MatSOR_SeqAIJSELL(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
191: {
192: Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL *)A->spptr;
194: PetscFunctionBegin;
195: PetscCall(MatSeqAIJSELL_build_shadow(A));
196: PetscCall(MatSOR_SeqSELL(aijsell->S, bb, omega, flag, fshift, its, lits, xx));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /* MatConvert_SeqAIJ_SeqAIJSELL converts a SeqAIJ matrix into a
201: * SeqAIJSELL matrix. This routine is called by the MatCreate_SeqAIJSELL()
202: * routine, but can also be used to convert an assembled SeqAIJ matrix
203: * into a SeqAIJSELL one. */
204: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJSELL(Mat A, MatType type, MatReuse reuse, Mat *newmat)
205: {
206: Mat B = *newmat;
207: Mat_SeqAIJ *b;
208: Mat_SeqAIJSELL *aijsell;
209: PetscBool set;
210: PetscBool sametype;
212: PetscFunctionBegin;
213: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
215: PetscCall(PetscObjectTypeCompare((PetscObject)A, type, &sametype));
216: if (sametype) PetscFunctionReturn(PETSC_SUCCESS);
218: PetscCall(PetscNew(&aijsell));
219: b = (Mat_SeqAIJ *)B->data;
220: B->spptr = (void *)aijsell;
222: /* Disable use of the inode routines so that the AIJSELL ones will be used instead.
223: * This happens in MatAssemblyEnd_SeqAIJSELL as well, but the assembly end may not be called, so set it here, too.
224: * As noted elsewhere, I wonder if it might make sense and be feasible to use some of the inode routines. */
225: b->inode.use = PETSC_FALSE;
227: /* Set function pointers for methods that we inherit from AIJ but override.
228: * We also parse some command line options below, since those determine some of the methods we point to. */
229: B->ops->duplicate = MatDuplicate_SeqAIJSELL;
230: B->ops->assemblyend = MatAssemblyEnd_SeqAIJSELL;
231: B->ops->destroy = MatDestroy_SeqAIJSELL;
233: aijsell->S = NULL;
234: aijsell->eager_shadow = PETSC_FALSE;
236: /* Parse command line options. */
237: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "AIJSELL Options", "Mat");
238: PetscCall(PetscOptionsBool("-mat_aijsell_eager_shadow", "Eager Shadowing", "None", (PetscBool)aijsell->eager_shadow, (PetscBool *)&aijsell->eager_shadow, &set));
239: PetscOptionsEnd();
241: /* If A has already been assembled and eager shadowing is specified, build the shadow matrix. */
242: if (A->assembled && aijsell->eager_shadow) PetscCall(MatSeqAIJSELL_build_shadow(A));
244: B->ops->mult = MatMult_SeqAIJSELL;
245: B->ops->multtranspose = MatMultTranspose_SeqAIJSELL;
246: B->ops->multadd = MatMultAdd_SeqAIJSELL;
247: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJSELL;
248: B->ops->sor = MatSOR_SeqAIJSELL;
250: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaijsell_seqaij_C", MatConvert_SeqAIJSELL_SeqAIJ));
252: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJSELL));
253: *newmat = B;
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@C
258: MatCreateSeqAIJSELL - Creates a sparse matrix of type `MATSEQAIJSELL`.
260: Collective
262: Input Parameters:
263: + comm - MPI communicator, set to `PETSC_COMM_SELF`
264: . m - number of rows
265: . n - number of columns
266: . nz - number of nonzeros per row (same for all rows)
267: - nnz - array containing the number of nonzeros in the various rows
268: (possibly different for each row) or `NULL`
270: Output Parameter:
271: . A - the matrix
273: Options Database Keys:
274: . -mat_aijsell_eager_shadow - Construct shadow matrix upon matrix assembly; default is to take a "lazy" approach,
275: performing this step the first time the matrix is applied
277: Level: intermediate
279: Notes:
280: This type inherits from AIJ and is largely identical, but keeps a "shadow" copy of the matrix
281: in `MATSEQSELL` format, which is used when this format may be more suitable for a requested
282: operation. Currently, `MATSEQSELL` format is used for `MatMult()`, `MatMultTranspose()`,
283: `MatMultAdd()`, `MatMultTransposeAdd()`, and `MatSOR()` operations.
285: If `nnz` is given then `nz` is ignored
287: Because `MATSEQAIJSELL` is a subtype of `MATSEQAIJ`, the option `-mat_seqaij_type seqaijsell` can be used to make
288: sequential `MATSEQAIJ` matrices default to being instances of `MATSEQAIJSELL`.
290: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateMPIAIJSELL()`, `MatSetValues()`
291: @*/
292: PetscErrorCode MatCreateSeqAIJSELL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
293: {
294: PetscFunctionBegin;
295: PetscCall(MatCreate(comm, A));
296: PetscCall(MatSetSizes(*A, m, n, m, n));
297: PetscCall(MatSetType(*A, MATSEQAIJSELL));
298: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJSELL(Mat A)
303: {
304: PetscFunctionBegin;
305: PetscCall(MatSetType(A, MATSEQAIJ));
306: PetscCall(MatConvert_SeqAIJ_SeqAIJSELL(A, MATSEQAIJSELL, MAT_INPLACE_MATRIX, &A));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }