Actual source code: gcreate.c

petsc-master 2016-08-24
Report Typos and Errors
 2:  #include <petsc/private/matimpl.h>

  6: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y,PetscScalar a)
  7: {
  9:   PetscInt       i,start,end;
 10:   PetscScalar    alpha = a;
 11:   PetscBool      prevoption;

 14:   MatGetOption(Y,MAT_NO_OFF_PROC_ENTRIES,&prevoption);
 15:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
 16:   MatGetOwnershipRange(Y,&start,&end);
 17:   for (i=start; i<end; i++) {
 18:     MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
 19:   }
 20:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
 21:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
 22:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);
 23:   return(0);
 24: }

 28: /*@
 29:    MatCreate - Creates a matrix where the type is determined
 30:    from either a call to MatSetType() or from the options database
 31:    with a call to MatSetFromOptions(). The default matrix type is
 32:    AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ()
 33:    if you do not set a type in the options database. If you never
 34:    call MatSetType() or MatSetFromOptions() it will generate an
 35:    error when you try to use the matrix.

 37:    Collective on MPI_Comm

 39:    Input Parameter:
 40: .  comm - MPI communicator

 42:    Output Parameter:
 43: .  A - the matrix

 45:    Options Database Keys:
 46: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
 47: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
 48: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
 49: .    -mat_type mpidense - dense type, uses MatCreateDense()
 50: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
 51: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

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

 57:    Notes:

 59:    Level: beginner

 61:    User manual sections:
 62: +   Section 3.1 Creating and Assembling Matrices
 63: -   Chapter 3 Matrices

 65: .keywords: matrix, create

 67: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
 68:           MatCreateSeqDense(), MatCreateDense(),
 69:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
 70:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
 71:           MatConvert()
 72: @*/
 73: PetscErrorCode  MatCreate(MPI_Comm comm,Mat *A)
 74: {
 75:   Mat            B;


 81:   *A = NULL;
 82:   MatInitializePackage();

 84:   PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
 85:   PetscLayoutCreate(comm,&B->rmap);
 86:   PetscLayoutCreate(comm,&B->cmap);

 88:   B->congruentlayouts = -1;
 89:   B->preallocated     = PETSC_FALSE;
 90:   *A                  = B;
 91:   return(0);
 92: }

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

 99:    Logically Collective on Mat

101:    Input Parameters:
102: +  mat -  matrix obtained from MatCreate()
103: -  flg - PETSC_TRUE indicates you want the error generated

105:    Level: advanced

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

109: .seealso: PCSetErrorIfFailure()
110: @*/
111: PetscErrorCode  MatSetErrorIfFailure(Mat mat,PetscBool flg)
112: {
116:   mat->erroriffailure = flg;
117:   return(0);
118: }

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

125:   Collective on Mat

127:   Input Parameters:
128: +  A - the matrix
129: .  m - number of local rows (or PETSC_DECIDE)
130: .  n - number of local columns (or PETSC_DECIDE)
131: .  M - number of global rows (or PETSC_DETERMINE)
132: -  N - number of global columns (or PETSC_DETERMINE)

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

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

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

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

150:   Level: beginner

152: .seealso: MatGetSize(), PetscSplitOwnership()
153: @*/
154: PetscErrorCode  MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
155: {
160:   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);
161:   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);
162:   if ((A->rmap->n >= 0 && A->rmap->N >= 0) && (A->rmap->n != m || 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);
163:   if ((A->cmap->n >= 0 && A->cmap->N >= 0) && (A->cmap->n != n || 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);
164:   A->rmap->n = m;
165:   A->cmap->n = n;
166:   A->rmap->N = M;
167:   A->cmap->N = N;
168:   return(0);
169: }

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

180:    Collective on Mat

182:    Input Parameter:
183: .  A - the matrix

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

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

197:    Level: beginner

199: .keywords: matrix, create

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


217:   PetscObjectOptionsBegin((PetscObject)B);

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

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

235:   PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
236:   PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
237:   PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);

239:   if (B->ops->setfromoptions) {
240:     (*B->ops->setfromoptions)(PetscOptionsObject,B);
241:   }

243:   flg  = PETSC_FALSE;
244:   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);
245:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
246:   flg  = PETSC_FALSE;
247:   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);
248:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}

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

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

261:    Collective on Mat

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

271:    Level: beginner

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

283:   MatSetBlockSize(A,bs);
284:   MatGetBlockSize(A,&bs);
285:   PetscLayoutSetUp(A->rmap);
286:   PetscLayoutSetUp(A->cmap);
287:   MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
288:   MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
289:   MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
290:   MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
291:   /*
292:     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
293:     good before going on with it.
294:   */
295:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
296:   PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
297:   if (!aij && !is) {
298:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
299:   }
300:   if (aij || is) {
301:     if (bs == 1) {
302:       MatSeqAIJSetPreallocation(A,0,dnnz);
303:       MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
304:       MatISSetPreallocation(A,0,dnnz,0,onnz);
305:     } else {                    /* Convert block-row precallocation to scalar-row */
306:       PetscInt i,m,*sdnnz,*sonnz;
307:       MatGetLocalSize(A,&m,NULL);
308:       PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
309:       for (i=0; i<m; i++) {
310:         if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
311:         if (onnz) sonnz[i] = onnz[i/bs] * bs;
312:       }
313:       MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
314:       MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
315:       MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
316:       PetscFree2(sdnnz,sonnz);
317:     }
318:   }
319:   return(0);
320: }

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

325:         This is somewhat different from MatHeaderReplace() it would be nice to merge the code
326: */
329: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
330: {
332:   PetscInt       refct;
333:   PetscOps       Abops;
334:   struct _MatOps Aops;
335:   char           *mtype,*mname;
336:   void           *spptr;

339:   /* save the parts of A we need */
340:   Abops = ((PetscObject)A)->bops[0];
341:   Aops  = A->ops[0];
342:   refct = ((PetscObject)A)->refct;
343:   mtype = ((PetscObject)A)->type_name;
344:   mname = ((PetscObject)A)->name;
345:   spptr = A->spptr;

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

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

354:   PetscFree((*C)->spptr);

356:   PetscLayoutDestroy(&A->rmap);
357:   PetscLayoutDestroy(&A->cmap);
358:   PetscFunctionListDestroy(&((PetscObject)A)->qlist);
359:   PetscObjectListDestroy(&((PetscObject)A)->olist);

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

364:   /* return the parts of A we saved */
365:   ((PetscObject)A)->bops[0]   = Abops;
366:   A->ops[0]                   = Aops;
367:   ((PetscObject)A)->refct     = refct;
368:   ((PetscObject)A)->type_name = mtype;
369:   ((PetscObject)A)->name      = mname;
370:   A->spptr                    = spptr;

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

376:   PetscHeaderDestroy(C);
377:   return(0);
378: }
379: /*
380:         Replace A's header with that of C; the C object is then destroyed

382:         This is essentially code moved from MatDestroy()

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

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

400:   if (A == *C) return(0);
402:   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);

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

413:   ((PetscObject)*C)->refct = 1;
414:   MatDestroy(C);
415:   return(0);
416: }