Actual source code: gcreate.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/
4: #if 0
7: static PetscErrorCode MatPublish_Base(PetscObject obj)
8: {
10: return(0);
11: }
12: #endif
16: /*@
17: MatCreate - Creates a matrix where the type is determined
18: from either a call to MatSetType() or from the options database
19: with a call to MatSetFromOptions(). The default matrix type is
20: AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ()
21: if you do not set a type in the options database. If you never
22: call MatSetType() or MatSetFromOptions() it will generate an
23: error when you try to use the matrix.
25: Collective on MPI_Comm
27: Input Parameter:
28: . comm - MPI communicator
29:
30: Output Parameter:
31: . A - the matrix
33: Options Database Keys:
34: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
35: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
36: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
37: . -mat_type mpidense - dense type, uses MatCreateDense()
38: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
39: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
41: Even More Options Database Keys:
42: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
43: for additional format-specific options.
45: Notes:
47: Level: beginner
49: User manual sections:
50: + Section 3.1 Creating and Assembling Matrices
51: - Chapter 3 Matrices
53: .keywords: matrix, create
55: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
56: MatCreateSeqDense(), MatCreateDense(),
57: MatCreateSeqBAIJ(), MatCreateBAIJ(),
58: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
59: MatConvert()
60: @*/
61: PetscErrorCode MatCreate(MPI_Comm comm,Mat *A)
62: {
63: Mat B;
69: *A = PETSC_NULL;
70: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
71: MatInitializePackage(PETSC_NULL);
72: #endif
74: PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_CLASSID,0,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
75: PetscLayoutCreate(comm,&B->rmap);
76: PetscLayoutCreate(comm,&B->cmap);
77: B->preallocated = PETSC_FALSE;
78: *A = B;
79: return(0);
80: }
84: /*@
85: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
87: Collective on Mat
89: Input Parameters:
90: + A - the matrix
91: . m - number of local rows (or PETSC_DECIDE)
92: . n - number of local columns (or PETSC_DECIDE)
93: . M - number of global rows (or PETSC_DETERMINE)
94: - N - number of global columns (or PETSC_DETERMINE)
96: Notes:
97: m (n) and M (N) cannot be both PETSC_DECIDE
98: If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.
100: If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
101: user must ensure that they are chosen to be compatible with the
102: vectors. To do this, one first considers the matrix-vector product
103: 'y = A x'. The 'm' that is used in the above routine must match the
104: local size used in the vector creation routine VecCreateMPI() for 'y'.
105: Likewise, the 'n' used must match that used as the local size in
106: VecCreateMPI() for 'x'.
108: You cannot change the sizes once they have been set.
110: The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.
112: Level: beginner
114: .seealso: MatGetSize(), PetscSplitOwnership()
115: @*/
116: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
117: {
120: 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);
121: 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);
122: 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);
123: 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);
124: A->rmap->n = m;
125: A->cmap->n = n;
126: A->rmap->N = M;
127: A->cmap->N = N;
128: return(0);
129: }
133: /*@
134: MatSetFromOptions - Creates a matrix where the type is determined
135: from the options database. Generates a parallel MPI matrix if the
136: communicator has more than one processor. The default matrix type is
137: AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
138: you do not select a type in the options database.
140: Collective on Mat
142: Input Parameter:
143: . A - the matrix
145: Options Database Keys:
146: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
147: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
148: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
149: . -mat_type mpidense - dense type, uses MatCreateDense()
150: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
151: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
153: Even More Options Database Keys:
154: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
155: for additional format-specific options.
157: Level: beginner
159: .keywords: matrix, create
161: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
162: MatCreateSeqDense(), MatCreateDense(),
163: MatCreateSeqBAIJ(), MatCreateBAIJ(),
164: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
165: MatConvert()
166: @*/
167: PetscErrorCode MatSetFromOptions(Mat B)
168: {
170: const char *deft = MATAIJ;
171: char type[256];
172: PetscBool flg,set;
177: PetscObjectOptionsBegin((PetscObject)B);
179: if (B->rmap->bs < 0) {
180: PetscInt newbs = -1;
181: PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
182: if (flg) {
183: PetscLayoutSetBlockSize(B->rmap,newbs);
184: PetscLayoutSetBlockSize(B->cmap,newbs);
185: }
186: }
188: PetscOptionsList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
189: if (flg) {
190: MatSetType(B,type);
191: } else if (!((PetscObject)B)->type_name) {
192: MatSetType(B,deft);
193: }
195: if (B->ops->setfromoptions) {
196: (*B->ops->setfromoptions)(B);
197: }
199: flg = PETSC_FALSE;
200: 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);
201: if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
202: flg = PETSC_FALSE;
203: 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);
204: if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
206: /* process any options handlers added with PetscObjectAddOptionsHandler() */
207: PetscObjectProcessOptionsHandlers((PetscObject)B);
208: PetscOptionsEnd();
210: return(0);
211: }
215: /*@
216: MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices
218: Collective on Mat
220: Input Arguments:
221: + A - matrix being preallocated
222: . bs - block size
223: . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
224: . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
225: . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
226: - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix
228: Level: beginner
230: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation()
231: @*/
232: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt *dnnz,const PetscInt *onnz,const PetscInt *dnnzu,const PetscInt *onnzu)
233: {
235: void (*aij)(void);
238: MatSetBlockSize(A,bs);
239: PetscLayoutSetUp(A->rmap);
240: PetscLayoutSetUp(A->cmap);
241: MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
242: MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
243: MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
244: MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
245: /*
246: In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
247: good before going on with it.
248: */
249: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
250: if (!aij) {
251: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
252: }
253: if (aij) {
254: if (bs == 1) {
255: MatSeqAIJSetPreallocation(A,0,dnnz);
256: MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
257: } else { /* Convert block-row precallocation to scalar-row */
258: PetscInt i,m,*sdnnz,*sonnz;
259: MatGetLocalSize(A,&m,PETSC_NULL);
260: PetscMalloc2((!!dnnz)*m,PetscInt,&sdnnz,(!!onnz)*m,PetscInt,&sonnz);
261: for (i=0; i<m; i++) {
262: if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
263: if (onnz) sonnz[i] = onnz[i/bs] * bs;
264: }
265: MatSeqAIJSetPreallocation(A,0,dnnz?sdnnz:PETSC_NULL);
266: MatMPIAIJSetPreallocation(A,0,dnnz?sdnnz:PETSC_NULL,0,onnz?sonnz:PETSC_NULL);
267: PetscFree2(sdnnz,sonnz);
268: }
269: }
270: return(0);
271: }
273: /*
274: Merges some information from Cs header to A; the C object is then destroyed
276: This is somewhat different from MatHeaderReplace() it would be nice to merge the code
277: */
280: PetscErrorCode MatHeaderMerge(Mat A,Mat C)
281: {
282: PetscErrorCode ierr;
283: PetscInt refct;
284: PetscOps *Abops;
285: MatOps Aops;
286: char *mtype,*mname;
287: void *spptr;
290: /* save the parts of A we need */
291: Abops = ((PetscObject)A)->bops;
292: Aops = A->ops;
293: refct = ((PetscObject)A)->refct;
294: mtype = ((PetscObject)A)->type_name;
295: mname = ((PetscObject)A)->name;
296: spptr = A->spptr;
298: /* zero these so the destroy below does not free them */
299: ((PetscObject)A)->type_name = 0;
300: ((PetscObject)A)->name = 0;
302: /* free all the interior data structures from mat */
303: (*A->ops->destroy)(A);
305: PetscFree(C->spptr);
307: PetscLayoutDestroy(&A->rmap);
308: PetscLayoutDestroy(&A->cmap);
309: PetscFListDestroy(&((PetscObject)A)->qlist);
310: PetscOListDestroy(&((PetscObject)A)->olist);
312: /* copy C over to A */
313: PetscMemcpy(A,C,sizeof(struct _p_Mat));
315: /* return the parts of A we saved */
316: ((PetscObject)A)->bops = Abops;
317: A->ops = Aops;
318: ((PetscObject)A)->refct = refct;
319: ((PetscObject)A)->type_name = mtype;
320: ((PetscObject)A)->name = mname;
321: A->spptr = spptr;
323: /* since these two are copied into A we do not want them destroyed in C */
324: ((PetscObject)C)->qlist = 0;
325: ((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
335: */
338: PetscErrorCode MatHeaderReplace(Mat A,Mat C)
339: {
341: PetscInt refct;
346: if (A == C) return(0);
348: if (((PetscObject)C)->refct != 1) SETERRQ1(((PetscObject)C)->comm,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: PetscMemcpy(A,C,sizeof(struct _p_Mat));
361: ((PetscObject)A)->refct = refct;
362: PetscLogObjectDestroy((PetscObject)C);
363: PetscFree(C);
364: return(0);
365: }