Actual source code: gcreate.c
petsc-master 2019-11-16
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
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: Level: beginner
66: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
67: MatCreateSeqDense(), MatCreateDense(),
68: MatCreateSeqBAIJ(), MatCreateBAIJ(),
69: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
70: MatConvert()
71: @*/
72: PetscErrorCode MatCreate(MPI_Comm comm,Mat *A)
73: {
74: Mat B;
80: *A = NULL;
81: MatInitializePackage();
83: PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
84: PetscLayoutCreate(comm,&B->rmap);
85: PetscLayoutCreate(comm,&B->cmap);
86: PetscStrallocpy(VECSTANDARD,&B->defaultvectype);
88: B->congruentlayouts = PETSC_DECIDE;
89: B->preallocated = PETSC_FALSE;
90: *A = B;
91: return(0);
92: }
94: /*@
95: MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected.
97: Logically Collective on Mat
99: Input Parameters:
100: + mat - matrix obtained from MatCreate()
101: - flg - PETSC_TRUE indicates you want the error generated
103: Level: advanced
105: .seealso: PCSetErrorIfFailure()
106: @*/
107: PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg)
108: {
112: mat->erroriffailure = flg;
113: return(0);
114: }
116: /*@
117: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
119: Collective on Mat
121: Input Parameters:
122: + A - the matrix
123: . m - number of local rows (or PETSC_DECIDE)
124: . n - number of local columns (or PETSC_DECIDE)
125: . M - number of global rows (or PETSC_DETERMINE)
126: - N - number of global columns (or PETSC_DETERMINE)
128: Notes:
129: m (n) and M (N) cannot be both PETSC_DECIDE
130: If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.
132: If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
133: user must ensure that they are chosen to be compatible with the
134: vectors. To do this, one first considers the matrix-vector product
135: 'y = A x'. The 'm' that is used in the above routine must match the
136: local size used in the vector creation routine VecCreateMPI() for 'y'.
137: Likewise, the 'n' used must match that used as the local size in
138: VecCreateMPI() for 'x'.
140: You cannot change the sizes once they have been set.
142: The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.
144: Level: beginner
146: .seealso: MatGetSize(), PetscSplitOwnership()
147: @*/
148: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
149: {
154: if (M > 0 && m > M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D cannot be larger than global row size %D",m,M);
155: if (N > 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D cannot be larger than global column size %D",n,N);
156: 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);
157: 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);
158: A->rmap->n = m;
159: A->cmap->n = n;
160: A->rmap->N = M > -1 ? M : A->rmap->N;
161: A->cmap->N = N > -1 ? N : A->cmap->N;
162: return(0);
163: }
165: /*@
166: MatSetFromOptions - Creates a matrix where the type is determined
167: from the options database. Generates a parallel MPI matrix if the
168: communicator has more than one processor. The default matrix type is
169: AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
170: you do not select a type in the options database.
172: Collective on Mat
174: Input Parameter:
175: . A - the matrix
177: Options Database Keys:
178: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
179: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
180: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
181: . -mat_type mpidense - dense type, uses MatCreateDense()
182: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
183: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
185: Even More Options Database Keys:
186: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
187: for additional format-specific options.
189: Level: beginner
191: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
192: MatCreateSeqDense(), MatCreateDense(),
193: MatCreateSeqBAIJ(), MatCreateBAIJ(),
194: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
195: MatConvert()
196: @*/
197: PetscErrorCode MatSetFromOptions(Mat B)
198: {
200: const char *deft = MATAIJ;
201: char type[256];
202: PetscBool flg,set;
207: PetscObjectOptionsBegin((PetscObject)B);
209: if (B->rmap->bs < 0) {
210: PetscInt newbs = -1;
211: PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
212: if (flg) {
213: PetscLayoutSetBlockSize(B->rmap,newbs);
214: PetscLayoutSetBlockSize(B->cmap,newbs);
215: }
216: }
218: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
219: if (flg) {
220: MatSetType(B,type);
221: } else if (!((PetscObject)B)->type_name) {
222: MatSetType(B,deft);
223: }
225: PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
226: PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
227: PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);
228: PetscOptionsBool("-mat_error_if_failure","Generate an error if an error occurs when factoring the matrix","MatSetErrorIfFailure",B->erroriffailure,&B->erroriffailure,NULL);
230: if (B->ops->setfromoptions) {
231: (*B->ops->setfromoptions)(PetscOptionsObject,B);
232: }
234: flg = PETSC_FALSE;
235: 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);
236: if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
237: flg = PETSC_FALSE;
238: 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);
239: if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
241: /* process any options handlers added with PetscObjectAddOptionsHandler() */
242: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);
243: PetscOptionsEnd();
244: return(0);
245: }
247: /*@C
248: MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions.
250: Collective on Mat
252: Input Arguments:
253: + A - matrix being preallocated
254: . bs - block size
255: . dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix
256: . onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
257: . dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
258: - onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix
260: Level: beginner
262: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
263: PetscSplitOwnership()
264: @*/
265: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
266: {
268: PetscInt cbs;
269: void (*aij)(void);
270: void (*is)(void);
271: void (*hyp)(void) = NULL;
274: if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
275: MatSetBlockSize(A,bs);
276: }
277: PetscLayoutSetUp(A->rmap);
278: PetscLayoutSetUp(A->cmap);
279: MatGetBlockSizes(A,&bs,&cbs);
280: /* these routines assumes bs == cbs, this should be checked somehow */
281: MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
282: MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
283: MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
284: MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
285: /*
286: 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
287: good before going on with it.
288: */
289: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
290: PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
291: #if defined(PETSC_HAVE_HYPRE)
292: PetscObjectQueryFunction((PetscObject)A,"MatHYPRESetPreallocation_C",&hyp);
293: #endif
294: if (!aij && !is && !hyp) {
295: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
296: }
297: if (aij || is || hyp) {
298: if (bs == cbs && bs == 1) {
299: MatSeqAIJSetPreallocation(A,0,dnnz);
300: MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
301: MatISSetPreallocation(A,0,dnnz,0,onnz);
302: #if defined(PETSC_HAVE_HYPRE)
303: MatHYPRESetPreallocation(A,0,dnnz,0,onnz);
304: #endif
305: } else { /* Convert block-row precallocation to scalar-row */
306: PetscInt i,m,*sdnnz,*sonnz;
307: MatGetLocalSize(A,&m,NULL);
308: PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
309: for (i=0; i<m; i++) {
310: if (dnnz) sdnnz[i] = dnnz[i/bs] * cbs;
311: if (onnz) sonnz[i] = onnz[i/bs] * cbs;
312: }
313: MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
314: MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
315: MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
316: #if defined(PETSC_HAVE_HYPRE)
317: MatHYPRESetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
318: #endif
319: PetscFree2(sdnnz,sonnz);
320: }
321: }
322: return(0);
323: }
325: /*
326: Merges some information from Cs header to A; the C object is then destroyed
328: This is somewhat different from MatHeaderReplace() it would be nice to merge the code
329: */
330: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
331: {
333: PetscInt refct;
334: PetscOps Abops;
335: struct _MatOps Aops;
336: char *mtype,*mname,*mprefix;
339: /* save the parts of A we need */
340: Abops = ((PetscObject)A)->bops[0];
341: Aops = A->ops[0];
342: refct = ((PetscObject)A)->refct;
343: mtype = ((PetscObject)A)->type_name;
344: mname = ((PetscObject)A)->name;
345: mprefix = ((PetscObject)A)->prefix;
347: /* zero these so the destroy below does not free them */
348: ((PetscObject)A)->type_name = 0;
349: ((PetscObject)A)->name = 0;
351: /* free all the interior data structures from mat */
352: (*A->ops->destroy)(A);
354: PetscFree(A->defaultvectype);
355: PetscLayoutDestroy(&A->rmap);
356: PetscLayoutDestroy(&A->cmap);
357: PetscFunctionListDestroy(&((PetscObject)A)->qlist);
358: PetscObjectListDestroy(&((PetscObject)A)->olist);
360: /* copy C over to A */
361: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
363: /* return the parts of A we saved */
364: ((PetscObject)A)->bops[0] = Abops;
365: A->ops[0] = Aops;
366: ((PetscObject)A)->refct = refct;
367: ((PetscObject)A)->type_name = mtype;
368: ((PetscObject)A)->name = mname;
369: ((PetscObject)A)->prefix = mprefix;
371: /* since these two are copied into A we do not want them destroyed in C */
372: ((PetscObject)*C)->qlist = 0;
373: ((PetscObject)*C)->olist = 0;
375: PetscHeaderDestroy(C);
376: return(0);
377: }
378: /*
379: Replace A's header with that of C; the C object is then destroyed
381: This is essentially code moved from MatDestroy()
383: This is somewhat different from MatHeaderMerge() it would be nice to merge the code
385: Used in DM hence is declared PETSC_EXTERN
386: */
387: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C)
388: {
389: PetscErrorCode ierr;
390: PetscInt refct;
391: PetscObjectState state;
392: struct _p_Mat buffer;
397: if (A == *C) return(0);
399: 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);
401: /* swap C and A */
402: refct = ((PetscObject)A)->refct;
403: state = ((PetscObject)A)->state;
404: PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));
405: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
406: PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));
407: ((PetscObject)A)->refct = refct;
408: ((PetscObject)A)->state = state + 1;
410: ((PetscObject)*C)->refct = 1;
411: MatDestroy(C);
412: return(0);
413: }
415: /*@
416: MatPinToCPU - marks a matrix to temporarily stay on the CPU and perform computations on the CPU
418: Input Parameters:
419: + A - the matrix
420: - flg - pin to the CPU if value of PETSC_TRUE
422: Level: intermediate
423: @*/
424: PetscErrorCode MatPinToCPU(Mat A,PetscBool flg)
425: {
426: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
430: if (A->pinnedtocpu == flg) return(0);
431: A->pinnedtocpu = flg;
432: if (A->ops->pintocpu) {
433: (*A->ops->pintocpu)(A,flg);
434: }
435: return(0);
436: #else
437: return 0;
438: #endif
439: }