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