Actual source code: gcreate.c
petsc-3.11.4 2019-09-28
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: }