Actual source code: gcreate.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
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->preallocated = PETSC_FALSE;
89: *A = B;
90: return(0);
91: }
95: /*@
96: MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected.
98: Logically Collective on Mat
100: Input Parameters:
101: + mat - matrix obtained from MatCreate()
102: - flg - PETSC_TRUE indicates you want the error generated
104: Level: advanced
106: .keywords: Mat, set, initial guess, nonzero
108: .seealso: PCSetErrorIfFailure()
109: @*/
110: PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg)
111: {
115: mat->erroriffailure = flg;
116: return(0);
117: }
121: /*@
122: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
124: Collective on Mat
126: Input Parameters:
127: + A - the matrix
128: . m - number of local rows (or PETSC_DECIDE)
129: . n - number of local columns (or PETSC_DECIDE)
130: . M - number of global rows (or PETSC_DETERMINE)
131: - N - number of global columns (or PETSC_DETERMINE)
133: Notes:
134: m (n) and M (N) cannot be both PETSC_DECIDE
135: If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.
137: If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
138: user must ensure that they are chosen to be compatible with the
139: vectors. To do this, one first considers the matrix-vector product
140: 'y = A x'. The 'm' that is used in the above routine must match the
141: local size used in the vector creation routine VecCreateMPI() for 'y'.
142: Likewise, the 'n' used must match that used as the local size in
143: VecCreateMPI() for 'x'.
145: You cannot change the sizes once they have been set.
147: The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.
149: Level: beginner
151: .seealso: MatGetSize(), PetscSplitOwnership()
152: @*/
153: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
154: {
159: 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);
160: 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);
161: 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);
162: 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);
163: A->rmap->n = m;
164: A->cmap->n = n;
165: A->rmap->N = M;
166: A->cmap->N = N;
167: return(0);
168: }
172: /*@
173: MatSetFromOptions - Creates a matrix where the type is determined
174: from the options database. Generates a parallel MPI matrix if the
175: communicator has more than one processor. The default matrix type is
176: AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
177: you do not select a type in the options database.
179: Collective on Mat
181: Input Parameter:
182: . A - the matrix
184: Options Database Keys:
185: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
186: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
187: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
188: . -mat_type mpidense - dense type, uses MatCreateDense()
189: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
190: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
192: Even More Options Database Keys:
193: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
194: for additional format-specific options.
196: Level: beginner
198: .keywords: matrix, create
200: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
201: MatCreateSeqDense(), MatCreateDense(),
202: MatCreateSeqBAIJ(), MatCreateBAIJ(),
203: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
204: MatConvert()
205: @*/
206: PetscErrorCode MatSetFromOptions(Mat B)
207: {
209: const char *deft = MATAIJ;
210: char type[256];
211: PetscBool flg,set;
216: PetscObjectOptionsBegin((PetscObject)B);
218: if (B->rmap->bs < 0) {
219: PetscInt newbs = -1;
220: PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
221: if (flg) {
222: PetscLayoutSetBlockSize(B->rmap,newbs);
223: PetscLayoutSetBlockSize(B->cmap,newbs);
224: }
225: }
227: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
228: if (flg) {
229: MatSetType(B,type);
230: } else if (!((PetscObject)B)->type_name) {
231: MatSetType(B,deft);
232: }
234: PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
235: PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
236: PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);
238: if (B->ops->setfromoptions) {
239: (*B->ops->setfromoptions)(PetscOptionsObject,B);
240: }
242: flg = PETSC_FALSE;
243: 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);
244: if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
245: flg = PETSC_FALSE;
246: 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);
247: if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
249: /* process any options handlers added with PetscObjectAddOptionsHandler() */
250: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);
251: PetscOptionsEnd();
252: return(0);
253: }
257: /*@
258: MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices
260: Collective on Mat
262: Input Arguments:
263: + A - matrix being preallocated
264: . bs - block size
265: . dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
266: . onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
267: . dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
268: - onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix
270: Level: beginner
272: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
273: PetscSplitOwnership()
274: @*/
275: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
276: {
278: void (*aij)(void);
281: MatSetBlockSize(A,bs);
282: MatGetBlockSize(A,&bs);
283: PetscLayoutSetUp(A->rmap);
284: PetscLayoutSetUp(A->cmap);
285: MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
286: MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
287: MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
288: MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
289: /*
290: In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
291: good before going on with it.
292: */
293: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
294: if (!aij) {
295: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
296: }
297: if (aij) {
298: if (bs == 1) {
299: MatSeqAIJSetPreallocation(A,0,dnnz);
300: MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
301: } else { /* Convert block-row precallocation to scalar-row */
302: PetscInt i,m,*sdnnz,*sonnz;
303: MatGetLocalSize(A,&m,NULL);
304: PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
305: for (i=0; i<m; i++) {
306: if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
307: if (onnz) sonnz[i] = onnz[i/bs] * bs;
308: }
309: MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
310: MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
311: PetscFree2(sdnnz,sonnz);
312: }
313: }
314: return(0);
315: }
317: /*
318: Merges some information from Cs header to A; the C object is then destroyed
320: This is somewhat different from MatHeaderReplace() it would be nice to merge the code
321: */
324: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
325: {
327: PetscInt refct;
328: PetscOps Abops;
329: struct _MatOps Aops;
330: char *mtype,*mname;
331: void *spptr;
334: /* save the parts of A we need */
335: Abops = ((PetscObject)A)->bops[0];
336: Aops = A->ops[0];
337: refct = ((PetscObject)A)->refct;
338: mtype = ((PetscObject)A)->type_name;
339: mname = ((PetscObject)A)->name;
340: spptr = A->spptr;
342: /* zero these so the destroy below does not free them */
343: ((PetscObject)A)->type_name = 0;
344: ((PetscObject)A)->name = 0;
346: /* free all the interior data structures from mat */
347: (*A->ops->destroy)(A);
349: PetscFree((*C)->spptr);
351: PetscLayoutDestroy(&A->rmap);
352: PetscLayoutDestroy(&A->cmap);
353: PetscFunctionListDestroy(&((PetscObject)A)->qlist);
354: PetscObjectListDestroy(&((PetscObject)A)->olist);
356: /* copy C over to A */
357: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
359: /* return the parts of A we saved */
360: ((PetscObject)A)->bops[0] = Abops;
361: A->ops[0] = Aops;
362: ((PetscObject)A)->refct = refct;
363: ((PetscObject)A)->type_name = mtype;
364: ((PetscObject)A)->name = mname;
365: A->spptr = spptr;
367: /* since these two are copied into A we do not want them destroyed in C */
368: ((PetscObject)*C)->qlist = 0;
369: ((PetscObject)*C)->olist = 0;
371: PetscHeaderDestroy(C);
372: return(0);
373: }
374: /*
375: Replace A's header with that of C; the C object is then destroyed
377: This is essentially code moved from MatDestroy()
379: This is somewhat different from MatHeaderMerge() it would be nice to merge the code
381: Used in DM hence is declared PETSC_EXTERN
382: */
385: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C)
386: {
387: PetscErrorCode ierr;
388: PetscInt refct;
389: PetscObjectState state;
390: struct _p_Mat buffer;
395: if (A == *C) return(0);
397: 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);
399: /* swap C and A */
400: refct = ((PetscObject)A)->refct;
401: state = ((PetscObject)A)->state;
402: PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));
403: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
404: PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));
405: ((PetscObject)A)->refct = refct;
406: ((PetscObject)A)->state = state + 1;
408: ((PetscObject)*C)->refct = 1;
409: MatDestroy(C);
410: return(0);
411: }