Actual source code: aijsell.c

petsc-3.13.4 2020-08-01
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,"MatMatMultSymbolic_seqdense_seqaijsell_C",NULL);
 42:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",NULL);

 44:   if (reuse == MAT_INITIAL_MATRIX) aijsell = (Mat_SeqAIJSELL*)B->spptr;

 46:   /* Clean up the Mat_SeqAIJSELL data structure.
 47:    * Note that MatDestroy() simply returns if passed a NULL value, so it's OK to call even if the shadow matrix was never constructed. */
 48:   MatDestroy(&aijsell->S);
 49:   PetscFree(B->spptr);

 51:   /* Change the type of B to MATSEQAIJ. */
 52:   PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ);

 54:   *newmat = B;
 55:   return(0);
 56: }

 58: PetscErrorCode MatDestroy_SeqAIJSELL(Mat A)
 59: {
 61:   Mat_SeqAIJSELL  *aijsell = (Mat_SeqAIJSELL*) A->spptr;


 65:   /* If MatHeaderMerge() was used, then this SeqAIJSELL matrix will not have an 
 66:    * spptr pointer. */
 67:   if (aijsell) {
 68:     /* Clean up everything in the Mat_SeqAIJSELL data structure, then free A->spptr. */
 69:     MatDestroy(&aijsell->S);
 70:     PetscFree(A->spptr);
 71:   }

 73:   /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ()
 74:    * to destroy everything that remains. */
 75:   PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ);
 76:   /* Note that I don't call MatSetType().  I believe this is because that
 77:    * is only to be called when *building* a matrix.  I could be wrong, but
 78:    * that is how things work for the SuperLU matrix class. */
 79:   MatDestroy_SeqAIJ(A);
 80:   return(0);
 81: }

 83: /* Build or update the shadow matrix if and only if needed.
 84:  * We track the ObjectState to determine when this needs to be done. */
 85: PETSC_INTERN PetscErrorCode MatSeqAIJSELL_build_shadow(Mat A)
 86: {
 87:   PetscErrorCode   ierr;
 88:   Mat_SeqAIJSELL   *aijsell = (Mat_SeqAIJSELL*) A->spptr;
 89:   PetscObjectState state;

 92:   PetscObjectStateGet((PetscObject)A,&state);
 93:   if (aijsell->S && aijsell->state == state) {
 94:     /* The existing shadow matrix is up-to-date, so simply exit. */
 95:     return(0);
 96:   }

 98:   PetscLogEventBegin(MAT_Convert,A,0,0,0);
 99:   if (aijsell->S) {
100:     MatConvert_SeqAIJ_SeqSELL(A,MATSEQSELL,MAT_REUSE_MATRIX,&aijsell->S);
101:   } else {
102:     MatConvert_SeqAIJ_SeqSELL(A,MATSEQSELL,MAT_INITIAL_MATRIX,&aijsell->S);
103:   }
104:   PetscLogEventEnd(MAT_Convert,A,0,0,0);

106:   /* Record the ObjectState so that we can tell when the shadow matrix needs updating */
107:   PetscObjectStateGet((PetscObject)A,&aijsell->state);

109:   return(0);
110: }

112: PetscErrorCode MatDuplicate_SeqAIJSELL(Mat A, MatDuplicateOption op, Mat *M)
113: {
115:   Mat_SeqAIJSELL *aijsell;
116:   Mat_SeqAIJSELL *aijsell_dest;

119:   MatDuplicate_SeqAIJ(A,op,M);
120:   aijsell      = (Mat_SeqAIJSELL*) A->spptr;
121:   aijsell_dest = (Mat_SeqAIJSELL*) (*M)->spptr;
122:   PetscArraycpy(aijsell_dest,aijsell,1);
123:   /* We don't duplicate the shadow matrix -- that will be constructed as needed. */
124:   aijsell_dest->S = NULL;
125:   if (aijsell->eager_shadow) {
126:     MatSeqAIJSELL_build_shadow(A);
127:   }
128:   return(0);
129: }

131: PetscErrorCode MatAssemblyEnd_SeqAIJSELL(Mat A, MatAssemblyType mode)
132: {
133:   PetscErrorCode  ierr;
134:   Mat_SeqAIJ      *a = (Mat_SeqAIJ*)A->data;
135:   Mat_SeqAIJSELL  *aijsell = (Mat_SeqAIJSELL*)A->spptr;

138:   if (mode == MAT_FLUSH_ASSEMBLY) return(0);

140:   /* I disable the use of the inode routines so that the AIJSELL ones will be
141:    * used instead, but I wonder if it might make sense (and is feasible) to
142:    * use some of them. */
143:   a->inode.use = PETSC_FALSE;

145:   /* Since a MATSEQAIJSELL matrix is really just a MATSEQAIJ with some
146:    * extra information and some different methods, call the AssemblyEnd 
147:    * routine for a MATSEQAIJ.
148:    * I'm not sure if this is the best way to do this, but it avoids
149:    * a lot of code duplication. */

151:   MatAssemblyEnd_SeqAIJ(A, mode);

153:   /* If the user has requested "eager" shadowing, create the SELL shadow matrix (if needed; the function checks).
154:    * (The default is to take a "lazy" approach, deferring this until something like MatMult() is called.) */
155:   if (aijsell->eager_shadow) {
156:     MatSeqAIJSELL_build_shadow(A);
157:   }

159:   return(0);
160: }

162: PetscErrorCode MatMult_SeqAIJSELL(Mat A,Vec xx,Vec yy)
163: {
164:   Mat_SeqAIJSELL    *aijsell = (Mat_SeqAIJSELL*)A->spptr;
165:   PetscErrorCode    ierr;

168:   MatSeqAIJSELL_build_shadow(A);
169:   MatMult_SeqSELL(aijsell->S,xx,yy);
170:   return(0);
171: }

173: PetscErrorCode MatMultTranspose_SeqAIJSELL(Mat A,Vec xx,Vec yy)
174: {
175:   Mat_SeqAIJSELL    *aijsell=(Mat_SeqAIJSELL*)A->spptr;
176:   PetscErrorCode    ierr;

179:   MatSeqAIJSELL_build_shadow(A);
180:   MatMultTranspose_SeqSELL(aijsell->S,xx,yy);
181:   return(0);
182: }

184: PetscErrorCode MatMultAdd_SeqAIJSELL(Mat A,Vec xx,Vec yy,Vec zz)
185: {
186:   Mat_SeqAIJSELL    *aijsell=(Mat_SeqAIJSELL*)A->spptr;
187:   PetscErrorCode    ierr;

190:   MatSeqAIJSELL_build_shadow(A);
191:   MatMultAdd_SeqSELL(aijsell->S,xx,yy,zz);
192:   return(0);
193: }

195: PetscErrorCode MatMultTransposeAdd_SeqAIJSELL(Mat A,Vec xx,Vec yy,Vec zz)
196: {
197:   Mat_SeqAIJSELL    *aijsell=(Mat_SeqAIJSELL*)A->spptr;
198:   PetscErrorCode    ierr;

201:   MatSeqAIJSELL_build_shadow(A);
202:   MatMultTransposeAdd_SeqSELL(aijsell->S,xx,yy,zz);
203:   return(0);
204: }

206: PetscErrorCode MatSOR_SeqAIJSELL(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
207: {
208:   Mat_SeqAIJSELL    *aijsell=(Mat_SeqAIJSELL*)A->spptr;
209:   PetscErrorCode    ierr;

212:   MatSeqAIJSELL_build_shadow(A);
213:   MatSOR_SeqSELL(aijsell->S,bb,omega,flag,fshift,its,lits,xx);
214:   return(0);
215: }

217: /* MatConvert_SeqAIJ_SeqAIJSELL converts a SeqAIJ matrix into a
218:  * SeqAIJSELL matrix.  This routine is called by the MatCreate_SeqAIJSELL()
219:  * routine, but can also be used to convert an assembled SeqAIJ matrix
220:  * into a SeqAIJSELL one. */
221: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJSELL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
222: {
224:   Mat            B = *newmat;
225:   Mat_SeqAIJ     *b;
226:   Mat_SeqAIJSELL *aijsell;
227:   PetscBool      set;
228:   PetscBool      sametype;

231:   if (reuse == MAT_INITIAL_MATRIX) {
232:     MatDuplicate(A,MAT_COPY_VALUES,&B);
233:   }

235:   PetscObjectTypeCompare((PetscObject)A,type,&sametype);
236:   if (sametype) return(0);

238:   PetscNewLog(B,&aijsell);
239:   b        = (Mat_SeqAIJ*) B->data;
240:   B->spptr = (void*) aijsell;

242:   /* Disable use of the inode routines so that the AIJSELL ones will be used instead.
243:    * This happens in MatAssemblyEnd_SeqAIJSELL as well, but the assembly end may not be called, so set it here, too.
244:    * As noted elsewhere, I wonder if it might make sense and be feasible to use some of the inode routines. */
245:   b->inode.use = PETSC_FALSE;

247:   /* Set function pointers for methods that we inherit from AIJ but override. 
248:    * We also parse some command line options below, since those determine some of the methods we point to. */
249:   B->ops->duplicate        = MatDuplicate_SeqAIJSELL;
250:   B->ops->assemblyend      = MatAssemblyEnd_SeqAIJSELL;
251:   B->ops->destroy          = MatDestroy_SeqAIJSELL;

253:   aijsell->S = NULL;
254:   aijsell->eager_shadow = PETSC_FALSE;

256:   /* Parse command line options. */
257:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"AIJSELL Options","Mat");
258:   PetscOptionsBool("-mat_aijsell_eager_shadow","Eager Shadowing","None",(PetscBool)aijsell->eager_shadow,(PetscBool*)&aijsell->eager_shadow,&set);
259:   PetscOptionsEnd();

261:   /* If A has already been assembled and eager shadowing is specified, build the shadow matrix. */
262:   if (A->assembled && aijsell->eager_shadow) {
263:     MatSeqAIJSELL_build_shadow(A);
264:   }

266:   B->ops->mult             = MatMult_SeqAIJSELL;
267:   B->ops->multtranspose    = MatMultTranspose_SeqAIJSELL;
268:   B->ops->multadd          = MatMultAdd_SeqAIJSELL;
269:   B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJSELL;
270:   B->ops->sor              = MatSOR_SeqAIJSELL;

272:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijsell_seqaij_C",MatConvert_SeqAIJSELL_SeqAIJ);
273:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaijsell_C",MatMatMultSymbolic_SeqDense_SeqAIJ);
274:   PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",MatMatMultNumeric_SeqDense_SeqAIJ);

276:   PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJSELL);
277:   *newmat = B;
278:   return(0);
279: }

281: /*@C
282:    MatCreateSeqAIJSELL - Creates a sparse matrix of type SEQAIJSELL.
283:    This type inherits from AIJ and is largely identical, but keeps a "shadow"
284:    copy of the matrix in SEQSELL format, which is used when this format
285:    may be more suitable for a requested operation. Currently, SEQSELL format
286:    is used for MatMult, MatMultTranspose, MatMultAdd, MatMultTransposeAdd,
287:    and MatSOR operations.
288:    Because SEQAIJSELL is a subtype of SEQAIJ, the option "-mat_seqaij_type seqaijsell" can be used to make
289:    sequential AIJ matrices default to being instances of MATSEQAIJSELL.

291:    Collective

293:    Input Parameters:
294: +  comm - MPI communicator, set to PETSC_COMM_SELF
295: .  m - number of rows
296: .  n - number of columns
297: .  nz - number of nonzeros per row (same for all rows)
298: -  nnz - array containing the number of nonzeros in the various rows
299:          (possibly different for each row) or NULL

301:    Output Parameter:
302: .  A - the matrix

304:    Options Database Keys:
305: .  -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

307:    Notes:
308:    If nnz is given then nz is ignored

310:    Level: intermediate

312: .seealso: MatCreate(), MatCreateMPIAIJSELL(), MatSetValues()
313: @*/
314: PetscErrorCode  MatCreateSeqAIJSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
315: {

319:   MatCreate(comm,A);
320:   MatSetSizes(*A,m,n,m,n);
321:   MatSetType(*A,MATSEQAIJSELL);
322:   MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
323:   return(0);
324: }

326: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJSELL(Mat A)
327: {

331:   MatSetType(A,MATSEQAIJ);
332:   MatConvert_SeqAIJ_SeqAIJSELL(A,MATSEQAIJSELL,MAT_INPLACE_MATRIX,&A);
333:   return(0);
334: }