Actual source code: gcreate.c
petsc-3.4.5 2014-06-29
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: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
61: MatInitializePackage();
62: #endif
64: PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
65: PetscLayoutCreate(comm,&B->rmap);
66: PetscLayoutCreate(comm,&B->cmap);
68: B->preallocated = PETSC_FALSE;
69: *A = B;
70: return(0);
71: }
75: /*@
76: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
78: Collective on Mat
80: Input Parameters:
81: + A - the matrix
82: . m - number of local rows (or PETSC_DECIDE)
83: . n - number of local columns (or PETSC_DECIDE)
84: . M - number of global rows (or PETSC_DETERMINE)
85: - N - number of global columns (or PETSC_DETERMINE)
87: Notes:
88: m (n) and M (N) cannot be both PETSC_DECIDE
89: If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.
91: If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
92: user must ensure that they are chosen to be compatible with the
93: vectors. To do this, one first considers the matrix-vector product
94: 'y = A x'. The 'm' that is used in the above routine must match the
95: local size used in the vector creation routine VecCreateMPI() for 'y'.
96: Likewise, the 'n' used must match that used as the local size in
97: VecCreateMPI() for 'x'.
99: You cannot change the sizes once they have been set.
101: The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.
103: Level: beginner
105: .seealso: MatGetSize(), PetscSplitOwnership()
106: @*/
107: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
108: {
113: 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);
114: 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);
115: 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);
116: 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);
117: A->rmap->n = m;
118: A->cmap->n = n;
119: A->rmap->N = M;
120: A->cmap->N = N;
121: return(0);
122: }
126: /*@
127: MatSetFromOptions - Creates a matrix where the type is determined
128: from the options database. Generates a parallel MPI matrix if the
129: communicator has more than one processor. The default matrix type is
130: AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
131: you do not select a type in the options database.
133: Collective on Mat
135: Input Parameter:
136: . A - the matrix
138: Options Database Keys:
139: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
140: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
141: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
142: . -mat_type mpidense - dense type, uses MatCreateDense()
143: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
144: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
146: Even More Options Database Keys:
147: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
148: for additional format-specific options.
150: Level: beginner
152: .keywords: matrix, create
154: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
155: MatCreateSeqDense(), MatCreateDense(),
156: MatCreateSeqBAIJ(), MatCreateBAIJ(),
157: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
158: MatConvert()
159: @*/
160: PetscErrorCode MatSetFromOptions(Mat B)
161: {
163: const char *deft = MATAIJ;
164: char type[256];
165: PetscBool flg,set;
170: PetscObjectOptionsBegin((PetscObject)B);
172: if (B->rmap->bs < 0) {
173: PetscInt newbs = -1;
174: PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
175: if (flg) {
176: PetscLayoutSetBlockSize(B->rmap,newbs);
177: PetscLayoutSetBlockSize(B->cmap,newbs);
178: }
179: }
181: PetscOptionsList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
182: if (flg) {
183: MatSetType(B,type);
184: } else if (!((PetscObject)B)->type_name) {
185: MatSetType(B,deft);
186: }
188: PetscViewerDestroy(&B->viewonassembly);
189: PetscOptionsViewer("-mat_view","Display mat with the viewer on MatAssemblyEnd()","MatView",&B->viewonassembly,&B->viewformatonassembly,NULL);
190: PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
191: PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",0.0,&B->checksymmetrytol,NULL);
192: PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",PETSC_FALSE,&B->checknullspaceonassembly,NULL);
194: if (B->ops->setfromoptions) {
195: (*B->ops->setfromoptions)(B);
196: }
198: flg = PETSC_FALSE;
199: 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);
200: if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
201: flg = PETSC_FALSE;
202: 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);
203: if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
205: /* process any options handlers added with PetscObjectAddOptionsHandler() */
206: PetscObjectProcessOptionsHandlers((PetscObject)B);
207: PetscOptionsEnd();
208: return(0);
209: }
213: /*@
214: MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices
216: Collective on Mat
218: Input Arguments:
219: + A - matrix being preallocated
220: . bs - block size
221: . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
222: . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
223: . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
224: - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix
226: Level: beginner
228: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
229: PetscSplitOwnership()
230: @*/
231: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
232: {
234: void (*aij)(void);
237: MatSetBlockSize(A,bs);
238: PetscLayoutSetUp(A->rmap);
239: PetscLayoutSetUp(A->cmap);
240: MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
241: MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
242: MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
243: MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
244: /*
245: In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
246: good before going on with it.
247: */
248: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
249: if (!aij) {
250: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
251: }
252: if (aij) {
253: if (bs == 1) {
254: MatSeqAIJSetPreallocation(A,0,dnnz);
255: MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
256: } else { /* Convert block-row precallocation to scalar-row */
257: PetscInt i,m,*sdnnz,*sonnz;
258: MatGetLocalSize(A,&m,NULL);
259: PetscMalloc2((!!dnnz)*m,PetscInt,&sdnnz,(!!onnz)*m,PetscInt,&sonnz);
260: for (i=0; i<m; i++) {
261: if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
262: if (onnz) sonnz[i] = onnz[i/bs] * bs;
263: }
264: MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
265: MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
266: PetscFree2(sdnnz,sonnz);
267: }
268: }
269: return(0);
270: }
272: /*
273: Merges some information from Cs header to A; the C object is then destroyed
275: This is somewhat different from MatHeaderReplace() it would be nice to merge the code
276: */
279: PetscErrorCode MatHeaderMerge(Mat A,Mat C)
280: {
282: PetscInt refct;
283: PetscOps *Abops;
284: MatOps Aops;
285: char *mtype,*mname;
286: void *spptr;
289: /* save the parts of A we need */
290: Abops = ((PetscObject)A)->bops;
291: Aops = A->ops;
292: refct = ((PetscObject)A)->refct;
293: mtype = ((PetscObject)A)->type_name;
294: mname = ((PetscObject)A)->name;
295: spptr = A->spptr;
297: /* zero these so the destroy below does not free them */
298: ((PetscObject)A)->type_name = 0;
299: ((PetscObject)A)->name = 0;
301: /* free all the interior data structures from mat */
302: (*A->ops->destroy)(A);
304: PetscFree(C->spptr);
306: PetscLayoutDestroy(&A->rmap);
307: PetscLayoutDestroy(&A->cmap);
308: PetscFunctionListDestroy(&((PetscObject)A)->qlist);
309: PetscObjectListDestroy(&((PetscObject)A)->olist);
311: /* copy C over to A */
312: PetscMemcpy(A,C,sizeof(struct _p_Mat));
314: /* return the parts of A we saved */
315: ((PetscObject)A)->bops = Abops;
316: A->ops = Aops;
317: ((PetscObject)A)->refct = refct;
318: ((PetscObject)A)->type_name = mtype;
319: ((PetscObject)A)->name = mname;
320: A->spptr = spptr;
322: /* since these two are copied into A we do not want them destroyed in C */
323: ((PetscObject)C)->qlist = 0;
324: ((PetscObject)C)->olist = 0;
326: PetscHeaderDestroy(&C);
327: return(0);
328: }
329: /*
330: Replace A's header with that of C; the C object is then destroyed
332: This is essentially code moved from MatDestroy()
334: This is somewhat different from MatHeaderMerge() it would be nice to merge the code
336: Used in DM hence is declared PETSC_EXTERN
337: */
340: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C)
341: {
343: PetscInt refct;
348: if (A == C) return(0);
350: 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);
352: /* free all the interior data structures from mat */
353: (*A->ops->destroy)(A);
354: PetscHeaderDestroy_Private((PetscObject)A);
355: PetscFree(A->ops);
356: PetscLayoutDestroy(&A->rmap);
357: PetscLayoutDestroy(&A->cmap);
358: PetscFree(A->spptr);
360: /* copy C over to A */
361: refct = ((PetscObject)A)->refct;
362: PetscMemcpy(A,C,sizeof(struct _p_Mat));
364: ((PetscObject)A)->refct = refct;
366: PetscLogObjectDestroy((PetscObject)C);
367: PetscFree(C);
368: return(0);
369: }