Actual source code: gcreate.c

petsc-3.11.2 2019-05-18
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>

  4: PETSC_INTERN PetscErrorCode MatSetBlockSizes_Default(Mat mat,PetscInt rbs, PetscInt cbs)
  5: {
  7:   if (!mat->preallocated) return(0);
  8:   if (mat->rmap->bs > 0 && mat->rmap->bs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change row block size %D to %D\n",mat->rmap->bs,rbs);
  9:   if (mat->cmap->bs > 0 && mat->cmap->bs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change column block size %D to %D\n",mat->cmap->bs,cbs);
 10:   return(0);
 11: }

 13: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y,PetscScalar a)
 14: {
 16:   PetscInt       i,start,end;
 17:   PetscScalar    alpha = a;
 18:   PetscBool      prevoption;

 21:   MatGetOption(Y,MAT_NO_OFF_PROC_ENTRIES,&prevoption);
 22:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
 23:   MatGetOwnershipRange(Y,&start,&end);
 24:   for (i=start; i<end; i++) {
 25:     if (i < Y->cmap->N) {
 26:       MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
 27:     }
 28:   }
 29:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
 30:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
 31:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);
 32:   return(0);
 33: }

 35: /*@
 36:    MatCreate - Creates a matrix where the type is determined
 37:    from either a call to MatSetType() or from the options database
 38:    with a call to MatSetFromOptions(). The default matrix type is
 39:    AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ()
 40:    if you do not set a type in the options database. If you never
 41:    call MatSetType() or MatSetFromOptions() it will generate an
 42:    error when you try to use the matrix.

 44:    Collective on MPI_Comm

 46:    Input Parameter:
 47: .  comm - MPI communicator

 49:    Output Parameter:
 50: .  A - the matrix

 52:    Options Database Keys:
 53: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
 54: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
 55: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
 56: .    -mat_type mpidense - dense type, uses MatCreateDense()
 57: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
 58: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

 60:    Even More Options Database Keys:
 61:    See the manpages for particular formats (e.g., MatCreateSeqAIJ())
 62:    for additional format-specific options.

 64:    Notes:

 66:    Level: beginner

 68:    User manual sections:
 69: +   Section 3.1 Creating and Assembling Matrices
 70: -   Chapter 3 Matrices

 72: .keywords: matrix, create

 74: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
 75:           MatCreateSeqDense(), MatCreateDense(),
 76:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
 77:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
 78:           MatConvert()
 79: @*/
 80: PetscErrorCode  MatCreate(MPI_Comm comm,Mat *A)
 81: {
 82:   Mat            B;


 88:   *A = NULL;
 89:   MatInitializePackage();

 91:   PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
 92:   PetscLayoutCreate(comm,&B->rmap);
 93:   PetscLayoutCreate(comm,&B->cmap);
 94:   PetscStrallocpy(VECSTANDARD,&B->defaultvectype);

 96:   B->congruentlayouts = PETSC_DECIDE;
 97:   B->preallocated     = PETSC_FALSE;
 98:   *A                  = B;
 99:   return(0);
100: }

102: /*@
103:    MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected.

105:    Logically Collective on Mat

107:    Input Parameters:
108: +  mat -  matrix obtained from MatCreate()
109: -  flg - PETSC_TRUE indicates you want the error generated

111:    Level: advanced

113: .keywords: Mat, set, initial guess, nonzero

115: .seealso: PCSetErrorIfFailure()
116: @*/
117: PetscErrorCode  MatSetErrorIfFailure(Mat mat,PetscBool flg)
118: {
122:   mat->erroriffailure = flg;
123:   return(0);
124: }

126: /*@
127:   MatSetSizes - Sets the local and global sizes, and checks to determine compatibility

129:   Collective on Mat

131:   Input Parameters:
132: +  A - the matrix
133: .  m - number of local rows (or PETSC_DECIDE)
134: .  n - number of local columns (or PETSC_DECIDE)
135: .  M - number of global rows (or PETSC_DETERMINE)
136: -  N - number of global columns (or PETSC_DETERMINE)

138:    Notes:
139:    m (n) and M (N) cannot be both PETSC_DECIDE
140:    If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.

142:    If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
143:    user must ensure that they are chosen to be compatible with the
144:    vectors. To do this, one first considers the matrix-vector product
145:    'y = A x'. The 'm' that is used in the above routine must match the
146:    local size used in the vector creation routine VecCreateMPI() for 'y'.
147:    Likewise, the 'n' used must match that used as the local size in
148:    VecCreateMPI() for 'x'.

150:    You cannot change the sizes once they have been set.

152:    The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.

154:   Level: beginner

156: .seealso: MatGetSize(), PetscSplitOwnership()
157: @*/
158: PetscErrorCode  MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
159: {
164:   if (M > 0 && m > M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D cannot be larger than global column size %D",m,M);
165:   if (N > 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D cannot be larger than global row size %D",n,N);
166:   if ((A->rmap->n >= 0 && A->rmap->N >= 0) && (A->rmap->n != m || (M > 0 && A->rmap->N != M))) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset row sizes to %D local %D global after previously setting them to %D local %D global",m,M,A->rmap->n,A->rmap->N);
167:   if ((A->cmap->n >= 0 && A->cmap->N >= 0) && (A->cmap->n != n || (N > 0 && A->cmap->N != N))) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset column sizes to %D local %D global after previously setting them to %D local %D global",n,N,A->cmap->n,A->cmap->N);
168:   A->rmap->n = m;
169:   A->cmap->n = n;
170:   A->rmap->N = M > -1 ? M : A->rmap->N;
171:   A->cmap->N = N > -1 ? N : A->cmap->N;
172:   return(0);
173: }

175: /*@
176:    MatSetFromOptions - Creates a matrix where the type is determined
177:    from the options database. Generates a parallel MPI matrix if the
178:    communicator has more than one processor.  The default matrix type is
179:    AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
180:    you do not select a type in the options database.

182:    Collective on Mat

184:    Input Parameter:
185: .  A - the matrix

187:    Options Database Keys:
188: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
189: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
190: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
191: .    -mat_type mpidense - dense type, uses MatCreateDense()
192: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
193: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

195:    Even More Options Database Keys:
196:    See the manpages for particular formats (e.g., MatCreateSeqAIJ())
197:    for additional format-specific options.

199:    Level: beginner

201: .keywords: matrix, create

203: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
204:           MatCreateSeqDense(), MatCreateDense(),
205:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
206:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
207:           MatConvert()
208: @*/
209: PetscErrorCode  MatSetFromOptions(Mat B)
210: {
212:   const char     *deft = MATAIJ;
213:   char           type[256];
214:   PetscBool      flg,set;


219:   PetscObjectOptionsBegin((PetscObject)B);

221:   if (B->rmap->bs < 0) {
222:     PetscInt newbs = -1;
223:     PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
224:     if (flg) {
225:       PetscLayoutSetBlockSize(B->rmap,newbs);
226:       PetscLayoutSetBlockSize(B->cmap,newbs);
227:     }
228:   }

230:   PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
231:   if (flg) {
232:     MatSetType(B,type);
233:   } else if (!((PetscObject)B)->type_name) {
234:     MatSetType(B,deft);
235:   }

237:   PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
238:   PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
239:   PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);
240:   PetscOptionsBool("-mat_error_if_failure","Generate an error if an error occurs when factoring the matrix","MatSetErrorIfFailure",B->erroriffailure,&B->erroriffailure,NULL);

242:   if (B->ops->setfromoptions) {
243:     (*B->ops->setfromoptions)(PetscOptionsObject,B);
244:   }

246:   flg  = PETSC_FALSE;
247:   PetscOptionsBool("-mat_new_nonzero_location_err","Generate an error if new nonzeros are created in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);
248:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
249:   flg  = PETSC_FALSE;
250:   PetscOptionsBool("-mat_new_nonzero_allocation_err","Generate an error if new nonzeros are allocated in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);
251:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}

253:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
254:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);
255:   PetscOptionsEnd();
256:   return(0);
257: }

259: /*@C
260:    MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions.

262:    Collective on Mat

264:    Input Arguments:
265: +  A - matrix being preallocated
266: .  bs - block size
267: .  dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
268: .  onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
269: .  dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
270: -  onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

272:    Level: beginner

274: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
275:           PetscSplitOwnership()
276: @*/
277: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
278: {
280:   void           (*aij)(void);
281:   void           (*is)(void);
282:   void           (*hyp)(void) = NULL;

285:   MatSetBlockSize(A,bs);
286:   MatGetBlockSize(A,&bs);
287:   PetscLayoutSetUp(A->rmap);
288:   PetscLayoutSetUp(A->cmap);
289:   MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
290:   MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
291:   MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
292:   MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
293:   /*
294:     In general, we have to do extra work to preallocate for scalar (AIJ) or unassembled (IS) matrices so we check whether it will do any
295:     good before going on with it.
296:   */
297:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
298:   PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
299: #if defined(PETSC_HAVE_HYPRE)
300:   PetscObjectQueryFunction((PetscObject)A,"MatHYPRESetPreallocation_C",&hyp);
301: #endif
302:   if (!aij && !is && !hyp) {
303:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
304:   }
305:   if (aij || is || hyp) {
306:     if (bs == 1) {
307:       MatSeqAIJSetPreallocation(A,0,dnnz);
308:       MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
309:       MatISSetPreallocation(A,0,dnnz,0,onnz);
310: #if defined(PETSC_HAVE_HYPRE)
311:       MatHYPRESetPreallocation(A,0,dnnz,0,onnz);
312: #endif
313:     } else {                    /* Convert block-row precallocation to scalar-row */
314:       PetscInt i,m,*sdnnz,*sonnz;
315:       MatGetLocalSize(A,&m,NULL);
316:       PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
317:       for (i=0; i<m; i++) {
318:         if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
319:         if (onnz) sonnz[i] = onnz[i/bs] * bs;
320:       }
321:       MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
322:       MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
323:       MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
324: #if defined(PETSC_HAVE_HYPRE)
325:       MatHYPRESetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
326: #endif
327:       PetscFree2(sdnnz,sonnz);
328:     }
329:   }
330:   return(0);
331: }

333: /*
334:         Merges some information from Cs header to A; the C object is then destroyed

336:         This is somewhat different from MatHeaderReplace() it would be nice to merge the code
337: */
338: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
339: {
341:   PetscInt       refct;
342:   PetscOps       Abops;
343:   struct _MatOps Aops;
344:   char           *mtype,*mname;

347:   /* save the parts of A we need */
348:   Abops = ((PetscObject)A)->bops[0];
349:   Aops  = A->ops[0];
350:   refct = ((PetscObject)A)->refct;
351:   mtype = ((PetscObject)A)->type_name;
352:   mname = ((PetscObject)A)->name;

354:   /* zero these so the destroy below does not free them */
355:   ((PetscObject)A)->type_name = 0;
356:   ((PetscObject)A)->name      = 0;

358:   /* free all the interior data structures from mat */
359:   (*A->ops->destroy)(A);

361:   PetscFree(A->defaultvectype);
362:   PetscLayoutDestroy(&A->rmap);
363:   PetscLayoutDestroy(&A->cmap);
364:   PetscFunctionListDestroy(&((PetscObject)A)->qlist);
365:   PetscObjectListDestroy(&((PetscObject)A)->olist);

367:   /* copy C over to A */
368:   PetscMemcpy(A,*C,sizeof(struct _p_Mat));

370:   /* return the parts of A we saved */
371:   ((PetscObject)A)->bops[0]   = Abops;
372:   A->ops[0]                   = Aops;
373:   ((PetscObject)A)->refct     = refct;
374:   ((PetscObject)A)->type_name = mtype;
375:   ((PetscObject)A)->name      = mname;

377:   /* since these two are copied into A we do not want them destroyed in C */
378:   ((PetscObject)*C)->qlist = 0;
379:   ((PetscObject)*C)->olist = 0;

381:   PetscHeaderDestroy(C);
382:   return(0);
383: }
384: /*
385:         Replace A's header with that of C; the C object is then destroyed

387:         This is essentially code moved from MatDestroy()

389:         This is somewhat different from MatHeaderMerge() it would be nice to merge the code

391:         Used in DM hence is declared PETSC_EXTERN
392: */
393: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C)
394: {
395:   PetscErrorCode   ierr;
396:   PetscInt         refct;
397:   PetscObjectState state;
398:   struct _p_Mat    buffer;

403:   if (A == *C) return(0);
405:   if (((PetscObject)*C)->refct != 1) SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Object C has refct %D > 1, would leave hanging reference",((PetscObject)*C)->refct);

407:   /* swap C and A */
408:   refct = ((PetscObject)A)->refct;
409:   state = ((PetscObject)A)->state;
410:   PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));
411:   PetscMemcpy(A,*C,sizeof(struct _p_Mat));
412:   PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));
413:   ((PetscObject)A)->refct = refct;
414:   ((PetscObject)A)->state = state + 1;

416:   ((PetscObject)*C)->refct = 1;
417:   MatDestroy(C);
418:   return(0);
419: }