Actual source code: aijsell.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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: }