Actual source code: gcreate.c
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: #if defined(PETSC_HAVE_DEVICE)
91: B->boundtocpu = PETSC_TRUE;
92: #endif
93: *A = B;
94: return(0);
95: }
97: /*@
98: MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected.
100: Logically Collective on Mat
102: Input Parameters:
103: + mat - matrix obtained from MatCreate()
104: - flg - PETSC_TRUE indicates you want the error generated
106: Level: advanced
108: .seealso: PCSetErrorIfFailure()
109: @*/
110: PetscErrorCode MatSetErrorIfFailure(Mat mat,PetscBool flg)
111: {
115: mat->erroriffailure = flg;
116: return(0);
117: }
119: /*@
120: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
122: Collective on Mat
124: Input Parameters:
125: + A - the matrix
126: . m - number of local rows (or PETSC_DECIDE)
127: . n - number of local columns (or PETSC_DECIDE)
128: . M - number of global rows (or PETSC_DETERMINE)
129: - N - number of global columns (or PETSC_DETERMINE)
131: Notes:
132: m (n) and M (N) cannot be both PETSC_DECIDE
133: If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.
135: If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
136: user must ensure that they are chosen to be compatible with the
137: vectors. To do this, one first considers the matrix-vector product
138: 'y = A x'. The 'm' that is used in the above routine must match the
139: local size used in the vector creation routine VecCreateMPI() for 'y'.
140: Likewise, the 'n' used must match that used as the local size in
141: VecCreateMPI() for 'x'.
143: You cannot change the sizes once they have been set.
145: The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.
147: Level: beginner
149: .seealso: MatGetSize(), PetscSplitOwnership()
150: @*/
151: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
152: {
157: 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);
158: 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);
159: 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);
160: 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);
161: A->rmap->n = m;
162: A->cmap->n = n;
163: A->rmap->N = M > -1 ? M : A->rmap->N;
164: A->cmap->N = N > -1 ? N : A->cmap->N;
165: return(0);
166: }
168: /*@
169: MatSetFromOptions - Creates a matrix where the type is determined
170: from the options database. Generates a parallel MPI matrix if the
171: communicator has more than one processor. The default matrix type is
172: AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
173: you do not select a type in the options database.
175: Collective on Mat
177: Input Parameter:
178: . A - the matrix
180: Options Database Keys:
181: + -mat_type seqaij - AIJ type, uses MatCreateSeqAIJ()
182: . -mat_type mpiaij - AIJ type, uses MatCreateAIJ()
183: . -mat_type seqdense - dense type, uses MatCreateSeqDense()
184: . -mat_type mpidense - dense type, uses MatCreateDense()
185: . -mat_type seqbaij - block AIJ type, uses MatCreateSeqBAIJ()
186: - -mat_type mpibaij - block AIJ type, uses MatCreateBAIJ()
188: Even More Options Database Keys:
189: See the manpages for particular formats (e.g., MatCreateSeqAIJ())
190: for additional format-specific options.
192: Level: beginner
194: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
195: MatCreateSeqDense(), MatCreateDense(),
196: MatCreateSeqBAIJ(), MatCreateBAIJ(),
197: MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
198: MatConvert()
199: @*/
200: PetscErrorCode MatSetFromOptions(Mat B)
201: {
203: const char *deft = MATAIJ;
204: char type[256];
205: PetscBool flg,set;
210: PetscObjectOptionsBegin((PetscObject)B);
212: if (B->rmap->bs < 0) {
213: PetscInt newbs = -1;
214: PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
215: if (flg) {
216: PetscLayoutSetBlockSize(B->rmap,newbs);
217: PetscLayoutSetBlockSize(B->cmap,newbs);
218: }
219: }
221: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
222: if (flg) {
223: MatSetType(B,type);
224: } else if (!((PetscObject)B)->type_name) {
225: MatSetType(B,deft);
226: }
228: PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
229: PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
230: PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);
231: PetscOptionsBool("-mat_error_if_failure","Generate an error if an error occurs when factoring the matrix","MatSetErrorIfFailure",B->erroriffailure,&B->erroriffailure,NULL);
233: if (B->ops->setfromoptions) {
234: (*B->ops->setfromoptions)(PetscOptionsObject,B);
235: }
237: flg = PETSC_FALSE;
238: 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);
239: if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
240: flg = PETSC_FALSE;
241: 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);
242: if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
243: flg = PETSC_FALSE;
244: PetscOptionsBool("-mat_ignore_zero_entries","For AIJ/IS matrices this will stop zero values from creating a zero location in the matrix","MatSetOption",flg,&flg,&set);
245: if (set) {MatSetOption(B,MAT_IGNORE_ZERO_ENTRIES,flg);}
247: flg = PETSC_FALSE;
248: PetscOptionsBool("-mat_form_explicit_transpose","Hint to form an explicit transpose for operations like MatMultTranspose","MatSetOption",flg,&flg,&set);
249: if (set) {MatSetOption(B,MAT_FORM_EXPLICIT_TRANSPOSE,flg);}
251: /* process any options handlers added with PetscObjectAddOptionsHandler() */
252: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);
253: PetscOptionsEnd();
254: return(0);
255: }
257: /*@C
258: MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions.
260: Collective on Mat
262: Input Parameters:
263: + A - matrix being preallocated
264: . bs - block size
265: . dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix
266: . onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
267: . dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
268: - onnzu - number of nonzero column 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: PetscInt cbs;
279: void (*aij)(void);
280: void (*is)(void);
281: void (*hyp)(void) = NULL;
284: if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
285: MatSetBlockSize(A,bs);
286: }
287: PetscLayoutSetUp(A->rmap);
288: PetscLayoutSetUp(A->cmap);
289: MatGetBlockSizes(A,&bs,&cbs);
290: /* these routines assumes bs == cbs, this should be checked somehow */
291: MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
292: MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
293: MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
294: MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
295: /*
296: 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
297: good before going on with it.
298: */
299: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
300: PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
301: #if defined(PETSC_HAVE_HYPRE)
302: PetscObjectQueryFunction((PetscObject)A,"MatHYPRESetPreallocation_C",&hyp);
303: #endif
304: if (!aij && !is && !hyp) {
305: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
306: }
307: if (aij || is || hyp) {
308: if (bs == cbs && bs == 1) {
309: MatSeqAIJSetPreallocation(A,0,dnnz);
310: MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
311: MatISSetPreallocation(A,0,dnnz,0,onnz);
312: #if defined(PETSC_HAVE_HYPRE)
313: MatHYPRESetPreallocation(A,0,dnnz,0,onnz);
314: #endif
315: } else { /* Convert block-row precallocation to scalar-row */
316: PetscInt i,m,*sdnnz,*sonnz;
317: MatGetLocalSize(A,&m,NULL);
318: PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
319: for (i=0; i<m; i++) {
320: if (dnnz) sdnnz[i] = dnnz[i/bs] * cbs;
321: if (onnz) sonnz[i] = onnz[i/bs] * cbs;
322: }
323: MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
324: MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
325: MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
326: #if defined(PETSC_HAVE_HYPRE)
327: MatHYPRESetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
328: #endif
329: PetscFree2(sdnnz,sonnz);
330: }
331: }
332: return(0);
333: }
335: /*
336: Merges some information from Cs header to A; the C object is then destroyed
338: This is somewhat different from MatHeaderReplace() it would be nice to merge the code
339: */
340: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
341: {
343: PetscInt refct;
344: PetscOps Abops;
345: struct _MatOps Aops;
346: char *mtype,*mname,*mprefix;
347: Mat_Product *product;
352: if (A == *C) return(0);
354: /* save the parts of A we need */
355: Abops = ((PetscObject)A)->bops[0];
356: Aops = A->ops[0];
357: refct = ((PetscObject)A)->refct;
358: mtype = ((PetscObject)A)->type_name;
359: mname = ((PetscObject)A)->name;
360: mprefix = ((PetscObject)A)->prefix;
361: product = A->product;
363: /* zero these so the destroy below does not free them */
364: ((PetscObject)A)->type_name = NULL;
365: ((PetscObject)A)->name = NULL;
367: /* free all the interior data structures from mat */
368: (*A->ops->destroy)(A);
370: PetscFree(A->defaultvectype);
371: PetscLayoutDestroy(&A->rmap);
372: PetscLayoutDestroy(&A->cmap);
373: PetscFunctionListDestroy(&((PetscObject)A)->qlist);
374: PetscObjectListDestroy(&((PetscObject)A)->olist);
376: /* copy C over to A */
377: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
379: /* return the parts of A we saved */
380: ((PetscObject)A)->bops[0] = Abops;
381: A->ops[0] = Aops;
382: ((PetscObject)A)->refct = refct;
383: ((PetscObject)A)->type_name = mtype;
384: ((PetscObject)A)->name = mname;
385: ((PetscObject)A)->prefix = mprefix;
386: A->product = product;
388: /* since these two are copied into A we do not want them destroyed in C */
389: ((PetscObject)*C)->qlist = NULL;
390: ((PetscObject)*C)->olist = NULL;
392: PetscHeaderDestroy(C);
393: return(0);
394: }
395: /*
396: Replace A's header with that of C; the C object is then destroyed
398: This is essentially code moved from MatDestroy()
400: This is somewhat different from MatHeaderMerge() it would be nice to merge the code
402: Used in DM hence is declared PETSC_EXTERN
403: */
404: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C)
405: {
406: PetscErrorCode ierr;
407: PetscInt refct;
408: PetscObjectState state;
409: struct _p_Mat buffer;
410: MatStencilInfo stencil;
415: if (A == *C) return(0);
417: 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);
419: /* swap C and A */
420: refct = ((PetscObject)A)->refct;
421: state = ((PetscObject)A)->state;
422: stencil = A->stencil;
423: PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));
424: PetscMemcpy(A,*C,sizeof(struct _p_Mat));
425: PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));
426: ((PetscObject)A)->refct = refct;
427: ((PetscObject)A)->state = state + 1;
428: A->stencil = stencil;
430: ((PetscObject)*C)->refct = 1;
431: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))NULL);
432: MatDestroy(C);
433: return(0);
434: }
436: /*@
437: MatBindToCPU - marks a matrix to temporarily stay on the CPU and perform computations on the CPU
439: Logically collective on Mat
441: Input Parameters:
442: + A - the matrix
443: - flg - bind to the CPU if value of PETSC_TRUE
445: Level: intermediate
447: .seealso: MatBoundToCPU()
448: @*/
449: PetscErrorCode MatBindToCPU(Mat A,PetscBool flg)
450: {
454: #if defined(PETSC_HAVE_DEVICE)
455: if (A->boundtocpu == flg) return(0);
456: A->boundtocpu = flg;
457: if (A->ops->bindtocpu) {
459: (*A->ops->bindtocpu)(A,flg);
460: }
461: #endif
462: return(0);
463: }
465: /*@
466: MatBoundToCPU - query if a matrix is bound to the CPU
468: Input Parameter:
469: . A - the matrix
471: Output Parameter:
472: . flg - the logical flag
474: Level: intermediate
476: .seealso: MatBindToCPU()
477: @*/
478: PetscErrorCode MatBoundToCPU(Mat A,PetscBool *flg)
479: {
483: #if defined(PETSC_HAVE_DEVICE)
484: *flg = A->boundtocpu;
485: #else
486: *flg = PETSC_TRUE;
487: #endif
488: return(0);
489: }
491: PetscErrorCode MatSetValuesCOO_Basic(Mat A,const PetscScalar coo_v[],InsertMode imode)
492: {
493: IS is_coo_i,is_coo_j;
494: const PetscInt *coo_i,*coo_j;
495: PetscInt n,n_i,n_j;
496: PetscScalar zero = 0.;
500: PetscObjectQuery((PetscObject)A,"__PETSc_coo_i",(PetscObject*)&is_coo_i);
501: PetscObjectQuery((PetscObject)A,"__PETSc_coo_j",(PetscObject*)&is_coo_j);
502: if (!is_coo_i) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_COR,"Missing coo_i IS");
503: if (!is_coo_j) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_COR,"Missing coo_j IS");
504: ISGetLocalSize(is_coo_i,&n_i);
505: ISGetLocalSize(is_coo_j,&n_j);
506: if (n_i != n_j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_COR,"Wrong local size %D != %D",n_i,n_j);
507: ISGetIndices(is_coo_i,&coo_i);
508: ISGetIndices(is_coo_j,&coo_j);
509: if (imode != ADD_VALUES) {
510: MatZeroEntries(A);
511: }
512: for (n = 0; n < n_i; n++) {
513: MatSetValue(A,coo_i[n],coo_j[n],coo_v ? coo_v[n] : zero,ADD_VALUES);
514: }
515: ISRestoreIndices(is_coo_i,&coo_i);
516: ISRestoreIndices(is_coo_j,&coo_j);
517: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
518: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
519: return(0);
520: }
522: PetscErrorCode MatSetPreallocationCOO_Basic(Mat A,PetscInt ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
523: {
524: Mat preallocator;
525: IS is_coo_i,is_coo_j;
526: PetscScalar zero = 0.0;
527: PetscInt n;
531: PetscLayoutSetUp(A->rmap);
532: PetscLayoutSetUp(A->cmap);
533: MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
534: MatSetType(preallocator,MATPREALLOCATOR);
535: MatSetSizes(preallocator,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
536: MatSetLayouts(preallocator,A->rmap,A->cmap);
537: MatSetUp(preallocator);
538: for (n = 0; n < ncoo; n++) {
539: MatSetValue(preallocator,coo_i[n],coo_j[n],zero,INSERT_VALUES);
540: }
541: MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
542: MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);
543: MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);
544: MatDestroy(&preallocator);
545: ISCreateGeneral(PETSC_COMM_SELF,ncoo,coo_i,PETSC_COPY_VALUES,&is_coo_i);
546: ISCreateGeneral(PETSC_COMM_SELF,ncoo,coo_j,PETSC_COPY_VALUES,&is_coo_j);
547: PetscObjectCompose((PetscObject)A,"__PETSc_coo_i",(PetscObject)is_coo_i);
548: PetscObjectCompose((PetscObject)A,"__PETSc_coo_j",(PetscObject)is_coo_j);
549: ISDestroy(&is_coo_i);
550: ISDestroy(&is_coo_j);
551: return(0);
552: }
554: /*@
555: MatSetPreallocationCOO - set preallocation for matrices using a coordinate format of the entries
557: Collective on Mat
559: Input Parameters:
560: + A - matrix being preallocated
561: . ncoo - number of entries in the locally owned part of the parallel matrix
562: . coo_i - row indices
563: - coo_j - column indices
565: Level: beginner
567: Notes: Entries can be repeated, see MatSetValuesCOO(). Currently optimized for cuSPARSE matrices only.
569: .seealso: MatSetValuesCOO(), MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation()
570: @*/
571: PetscErrorCode MatSetPreallocationCOO(Mat A,PetscInt ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
572: {
573: PetscErrorCode (*f)(Mat,PetscInt,const PetscInt[],const PetscInt[]) = NULL;
581: PetscLayoutSetUp(A->rmap);
582: PetscLayoutSetUp(A->cmap);
583: if (PetscDefined(USE_DEBUG)) {
584: PetscInt i;
585: for (i = 0; i < ncoo; i++) {
586: if (coo_i[i] < A->rmap->rstart || coo_i[i] >= A->rmap->rend) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid row index %D! Must be in [%D,%D)",coo_i[i],A->rmap->rstart,A->rmap->rend);
587: if (coo_j[i] < 0 || coo_j[i] >= A->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid col index %D! Must be in [0,%D)",coo_j[i],A->cmap->N);
588: }
589: }
590: PetscObjectQueryFunction((PetscObject)A,"MatSetPreallocationCOO_C",&f);
591: PetscLogEventBegin(MAT_PreallCOO,A,0,0,0);
592: if (f) {
593: (*f)(A,ncoo,coo_i,coo_j);
594: } else { /* allow fallback, very slow */
595: MatSetPreallocationCOO_Basic(A,ncoo,coo_i,coo_j);
596: }
597: PetscLogEventEnd(MAT_PreallCOO,A,0,0,0);
598: return(0);
599: }
601: /*@
602: MatSetValuesCOO - set values at once in a matrix preallocated using MatSetPreallocationCOO()
604: Collective on Mat
606: Input Parameters:
607: + A - matrix being preallocated
608: . coo_v - the matrix values (can be NULL)
609: - imode - the insert mode
611: Level: beginner
613: Notes: The values must follow the order of the indices prescribed with MatSetPreallocationCOO().
614: When repeated entries are specified in the COO indices the coo_v values are first properly summed.
615: The imode flag indicates if coo_v must be added to the current values of the matrix (ADD_VALUES) or overwritten (INSERT_VALUES).
616: Currently optimized for cuSPARSE matrices only.
617: Passing coo_v == NULL is equivalent to passing an array of zeros.
619: .seealso: MatSetPreallocationCOO(), InsertMode, INSERT_VALUES, ADD_VALUES
620: @*/
621: PetscErrorCode MatSetValuesCOO(Mat A, const PetscScalar coo_v[], InsertMode imode)
622: {
623: PetscErrorCode (*f)(Mat,const PetscScalar[],InsertMode) = NULL;
629: MatCheckPreallocated(A,1);
631: PetscObjectQueryFunction((PetscObject)A,"MatSetValuesCOO_C",&f);
632: PetscLogEventBegin(MAT_SetVCOO,A,0,0,0);
633: if (f) {
634: (*f)(A,coo_v,imode);
635: } else { /* allow fallback */
636: MatSetValuesCOO_Basic(A,coo_v,imode);
637: }
638: PetscLogEventEnd(MAT_SetVCOO,A,0,0,0);
639: PetscObjectStateIncrease((PetscObject)A);
640: return(0);
641: }