Actual source code: gcreate.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  2: #include <petsc-private/matimpl.h>       /*I "petscmat.h"  I*/

  6: /*@
  7:    MatCreate - Creates a matrix where the type is determined
  8:    from either a call to MatSetType() or from the options database
  9:    with a call to MatSetFromOptions(). The default matrix type is
 10:    AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ()
 11:    if you do not set a type in the options database. If you never
 12:    call MatSetType() or MatSetFromOptions() it will generate an
 13:    error when you try to use the matrix.

 15:    Collective on MPI_Comm

 17:    Input Parameter:
 18: .  comm - MPI communicator

 20:    Output Parameter:
 21: .  A - the matrix

 23:    Options Database Keys:
 24: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
 25: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
 26: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
 27: .    -mat_type mpidense - dense type, uses MatCreateDense()
 28: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
 29: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

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

 35:    Notes:

 37:    Level: beginner

 39:    User manual sections:
 40: +   Section 3.1 Creating and Assembling Matrices
 41: -   Chapter 3 Matrices

 43: .keywords: matrix, create

 45: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
 46:           MatCreateSeqDense(), MatCreateDense(),
 47:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
 48:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
 49:           MatConvert()
 50: @*/
 51: PetscErrorCode  MatCreate(MPI_Comm comm,Mat *A)
 52: {
 53:   Mat            B;


 59:   *A = NULL;
 60:   MatInitializePackage();

 62:   PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
 63:   PetscLayoutCreate(comm,&B->rmap);
 64:   PetscLayoutCreate(comm,&B->cmap);

 66:   B->preallocated = PETSC_FALSE;
 67:   *A              = B;
 68:   return(0);
 69: }

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

 76:   Collective on Mat

 78:   Input Parameters:
 79: +  A - the matrix
 80: .  m - number of local rows (or PETSC_DECIDE)
 81: .  n - number of local columns (or PETSC_DECIDE)
 82: .  M - number of global rows (or PETSC_DETERMINE)
 83: -  N - number of global columns (or PETSC_DETERMINE)

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

 89:    If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
 90:    user must ensure that they are chosen to be compatible with the
 91:    vectors. To do this, one first considers the matrix-vector product
 92:    'y = A x'. The 'm' that is used in the above routine must match the
 93:    local size used in the vector creation routine VecCreateMPI() for 'y'.
 94:    Likewise, the 'n' used must match that used as the local size in
 95:    VecCreateMPI() for 'x'.

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

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

101:   Level: beginner

103: .seealso: MatGetSize(), PetscSplitOwnership()
104: @*/
105: PetscErrorCode  MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
106: {
111:   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);
112:   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);
113:   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);
114:   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);
115:   A->rmap->n = m;
116:   A->cmap->n = n;
117:   A->rmap->N = M;
118:   A->cmap->N = N;
119:   return(0);
120: }

124: /*@
125:    MatSetFromOptions - Creates a matrix where the type is determined
126:    from the options database. Generates a parallel MPI matrix if the
127:    communicator has more than one processor.  The default matrix type is
128:    AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
129:    you do not select a type in the options database.

131:    Collective on Mat

133:    Input Parameter:
134: .  A - the matrix

136:    Options Database Keys:
137: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
138: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
139: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
140: .    -mat_type mpidense - dense type, uses MatCreateDense()
141: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
142: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

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

148:    Level: beginner

150: .keywords: matrix, create

152: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
153:           MatCreateSeqDense(), MatCreateDense(),
154:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
155:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
156:           MatConvert()
157: @*/
158: PetscErrorCode  MatSetFromOptions(Mat B)
159: {
161:   const char     *deft = MATAIJ;
162:   char           type[256];
163:   PetscBool      flg,set;


168:   PetscObjectOptionsBegin((PetscObject)B);

170:   if (B->rmap->bs < 0) {
171:     PetscInt newbs = -1;
172:     PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
173:     if (flg) {
174:       PetscLayoutSetBlockSize(B->rmap,newbs);
175:       PetscLayoutSetBlockSize(B->cmap,newbs);
176:     }
177:   }

179:   PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
180:   if (flg) {
181:     MatSetType(B,type);
182:   } else if (!((PetscObject)B)->type_name) {
183:     MatSetType(B,deft);
184:   }

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

190:   if (B->ops->setfromoptions) {
191:     (*B->ops->setfromoptions)(B);
192:   }

194:   flg  = PETSC_FALSE;
195:   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);
196:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
197:   flg  = PETSC_FALSE;
198:   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);
199:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}

201:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
202:   PetscObjectProcessOptionsHandlers((PetscObject)B);
203:   PetscOptionsEnd();
204:   return(0);
205: }

209: /*@
210:    MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices

212:    Collective on Mat

214:    Input Arguments:
215: +  A - matrix being preallocated
216: .  bs - block size
217: .  dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
218: .  onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
219: .  dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
220: -  onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

222:    Level: beginner

224: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
225:           PetscSplitOwnership()
226: @*/
227: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
228: {
230:   void           (*aij)(void);

233:   MatSetBlockSize(A,bs);
234:   MatGetBlockSize(A,&bs);
235:   PetscLayoutSetUp(A->rmap);
236:   PetscLayoutSetUp(A->cmap);
237:   MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
238:   MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
239:   MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
240:   MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
241:   /*
242:     In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
243:     good before going on with it.
244:   */
245:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
246:   if (!aij) {
247:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
248:   }
249:   if (aij) {
250:     if (bs == 1) {
251:       MatSeqAIJSetPreallocation(A,0,dnnz);
252:       MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
253:     } else {                    /* Convert block-row precallocation to scalar-row */
254:       PetscInt i,m,*sdnnz,*sonnz;
255:       MatGetLocalSize(A,&m,NULL);
256:       PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
257:       for (i=0; i<m; i++) {
258:         if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
259:         if (onnz) sonnz[i] = onnz[i/bs] * bs;
260:       }
261:       MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
262:       MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
263:       PetscFree2(sdnnz,sonnz);
264:     }
265:   }
266:   return(0);
267: }

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

272:         This is somewhat different from MatHeaderReplace() it would be nice to merge the code
273: */
276: PetscErrorCode MatHeaderMerge(Mat A,Mat C)
277: {
279:   PetscInt       refct;
280:   PetscOps       *Abops;
281:   MatOps         Aops;
282:   char           *mtype,*mname;
283:   void           *spptr;

286:   /* save the parts of A we need */
287:   Abops = ((PetscObject)A)->bops;
288:   Aops  = A->ops;
289:   refct = ((PetscObject)A)->refct;
290:   mtype = ((PetscObject)A)->type_name;
291:   mname = ((PetscObject)A)->name;
292:   spptr = A->spptr;

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

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

301:   PetscFree(C->spptr);

303:   PetscLayoutDestroy(&A->rmap);
304:   PetscLayoutDestroy(&A->cmap);
305:   PetscFunctionListDestroy(&((PetscObject)A)->qlist);
306:   PetscObjectListDestroy(&((PetscObject)A)->olist);

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

311:   /* return the parts of A we saved */
312:   ((PetscObject)A)->bops      = Abops;
313:   A->ops                      = Aops;
314:   ((PetscObject)A)->refct     = refct;
315:   ((PetscObject)A)->type_name = mtype;
316:   ((PetscObject)A)->name      = mname;
317:   A->spptr                    = spptr;

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

323:   PetscHeaderDestroy(&C);
324:   return(0);
325: }
326: /*
327:         Replace A's header with that of C; the C object is then destroyed

329:         This is essentially code moved from MatDestroy()

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

333:         Used in DM hence is declared PETSC_EXTERN
334: */
337: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C)
338: {
339:   PetscErrorCode   ierr;
340:   PetscInt         refct;
341:   PetscObjectState state;

346:   if (A == C) return(0);
348:   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);

350:   /* free all the interior data structures from mat */
351:   (*A->ops->destroy)(A);
352:   PetscHeaderDestroy_Private((PetscObject)A);
353:   PetscFree(A->ops);
354:   PetscLayoutDestroy(&A->rmap);
355:   PetscLayoutDestroy(&A->cmap);
356:   PetscFree(A->spptr);

358:   /* copy C over to A */
359:   refct = ((PetscObject)A)->refct;
360:   state = ((PetscObject)A)->state;
361:   PetscMemcpy(A,C,sizeof(struct _p_Mat));

363:   ((PetscObject)A)->refct = refct;
364:   ((PetscObject)A)->state = state + 1;

366:   PetscFree(C);
367:   return(0);
368: }