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: }