Actual source code: matrix.c
1: #define PETSCMAT_DLL
3: /*
4: This is where the abstract matrix operations are defined
5: */
7: #include include/private/matimpl.h
8: #include private/vecimpl.h
10: /* Logging support */
11: PetscCookie MAT_COOKIE = 0;
12: PetscEvent MAT_Mult = 0, MAT_Mults = 0, MAT_MultConstrained = 0, MAT_MultAdd = 0, MAT_MultTranspose = 0;
13: PetscEvent MAT_MultTransposeConstrained = 0, MAT_MultTransposeAdd = 0, MAT_Solve = 0, MAT_Solves = 0, MAT_SolveAdd = 0, MAT_SolveTranspose = 0, MAT_MatSolve = 0;
14: PetscEvent MAT_SolveTransposeAdd = 0, MAT_Relax = 0, MAT_ForwardSolve = 0, MAT_BackwardSolve = 0, MAT_LUFactor = 0, MAT_LUFactorSymbolic = 0;
15: PetscEvent MAT_LUFactorNumeric = 0, MAT_CholeskyFactor = 0, MAT_CholeskyFactorSymbolic = 0, MAT_CholeskyFactorNumeric = 0, MAT_ILUFactor = 0;
16: PetscEvent MAT_ILUFactorSymbolic = 0, MAT_ICCFactorSymbolic = 0, MAT_Copy = 0, MAT_Convert = 0, MAT_Scale = 0, MAT_AssemblyBegin = 0;
17: PetscEvent MAT_AssemblyEnd = 0, MAT_SetValues = 0, MAT_GetValues = 0, MAT_GetRow = 0, MAT_GetRowIJ = 0, MAT_GetSubMatrices = 0, MAT_GetColoring = 0, MAT_GetOrdering = 0, MAT_GetRedundantMatrix = 0;
18: PetscEvent MAT_IncreaseOverlap = 0, MAT_Partitioning = 0, MAT_ZeroEntries = 0, MAT_Load = 0, MAT_View = 0, MAT_AXPY = 0, MAT_FDColoringCreate = 0;
19: PetscEvent MAT_FDColoringApply = 0,MAT_Transpose = 0,MAT_FDColoringFunction = 0;
20: PetscEvent MAT_MatMult = 0, MAT_MatMultSymbolic = 0, MAT_MatMultNumeric = 0;
21: PetscEvent MAT_PtAP = 0, MAT_PtAPSymbolic = 0, MAT_PtAPNumeric = 0;
22: PetscEvent MAT_MatMultTranspose = 0, MAT_MatMultTransposeSymbolic = 0, MAT_MatMultTransposeNumeric = 0;
24: /* nasty global values for MatSetValue() */
25: PetscInt MatSetValue_Row = 0;
26: PetscInt MatSetValue_Column = 0;
27: PetscScalar MatSetValue_Value = 0.0;
31: /*@
32: MatRealPart - Zeros out the imaginary part of the matrix
34: Collective on Mat
36: Input Parameters:
37: . mat - the matrix
39: Level: advanced
42: .seealso: MatImaginaryPart()
43: @*/
45: PetscErrorCode MatRealPart(Mat mat)
46: {
52: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
53: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
54: if (!mat->ops->realpart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
55: MatPreallocated(mat);
56: (*mat->ops->realpart)(mat);
57: return(0);
58: }
62: /*@
63: MatImaginaryPart - Moves the imaginary part of the matrix to the real part and zeros the imaginary part
65: Collective on Mat
67: Input Parameters:
68: . mat - the matrix
70: Level: advanced
73: .seealso: MatRealPart()
74: @*/
76: PetscErrorCode MatImaginaryPart(Mat mat)
77: {
83: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
84: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
85: if (!mat->ops->imaginarypart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
86: MatPreallocated(mat);
87: (*mat->ops->imaginarypart)(mat);
88: return(0);
89: }
93: /*@C
94: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
95: for each row that you get to ensure that your application does
96: not bleed memory.
98: Not Collective
100: Input Parameters:
101: + mat - the matrix
102: - row - the row to get
104: Output Parameters:
105: + ncols - if not NULL, the number of nonzeros in the row
106: . cols - if not NULL, the column numbers
107: - vals - if not NULL, the values
109: Notes:
110: This routine is provided for people who need to have direct access
111: to the structure of a matrix. We hope that we provide enough
112: high-level matrix routines that few users will need it.
114: MatGetRow() always returns 0-based column indices, regardless of
115: whether the internal representation is 0-based (default) or 1-based.
117: For better efficiency, set cols and/or vals to PETSC_NULL if you do
118: not wish to extract these quantities.
120: The user can only examine the values extracted with MatGetRow();
121: the values cannot be altered. To change the matrix entries, one
122: must use MatSetValues().
124: You can only have one call to MatGetRow() outstanding for a particular
125: matrix at a time, per processor. MatGetRow() can only obtain rows
126: associated with the given processor, it cannot get rows from the
127: other processors; for that we suggest using MatGetSubMatrices(), then
128: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
129: is in the global number of rows.
131: Fortran Notes:
132: The calling sequence from Fortran is
133: .vb
134: MatGetRow(matrix,row,ncols,cols,values,ierr)
135: Mat matrix (input)
136: integer row (input)
137: integer ncols (output)
138: integer cols(maxcols) (output)
139: double precision (or double complex) values(maxcols) output
140: .ve
141: where maxcols >= maximum nonzeros in any row of the matrix.
144: Caution:
145: Do not try to change the contents of the output arrays (cols and vals).
146: In some cases, this may corrupt the matrix.
148: Level: advanced
150: Concepts: matrices^row access
152: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubMatrices(), MatGetDiagonal()
153: @*/
155: PetscErrorCode MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
156: {
158: PetscInt incols;
163: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
164: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
165: if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
166: MatPreallocated(mat);
168: (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
169: if (ncols) *ncols = incols;
171: return(0);
172: }
176: /*@
177: MatConjugate - replaces the matrix values with their complex conjugates
179: Collective on Mat
181: Input Parameters:
182: . mat - the matrix
184: Level: advanced
186: .seealso: VecConjugate()
187: @*/
188: PetscErrorCode MatConjugate(Mat mat)
189: {
194: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
195: if (!mat->ops->conjugate) SETERRQ(PETSC_ERR_SUP,"Not provided for this matrix format, send email to petsc-maint@mcs.anl.gov");
196: (*mat->ops->conjugate)(mat);
197: return(0);
198: }
202: /*@C
203: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
205: Not Collective
207: Input Parameters:
208: + mat - the matrix
209: . row - the row to get
210: . ncols, cols - the number of nonzeros and their columns
211: - vals - if nonzero the column values
213: Notes:
214: This routine should be called after you have finished examining the entries.
216: Fortran Notes:
217: The calling sequence from Fortran is
218: .vb
219: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
220: Mat matrix (input)
221: integer row (input)
222: integer ncols (output)
223: integer cols(maxcols) (output)
224: double precision (or double complex) values(maxcols) output
225: .ve
226: Where maxcols >= maximum nonzeros in any row of the matrix.
228: In Fortran MatRestoreRow() MUST be called after MatGetRow()
229: before another call to MatGetRow() can be made.
231: Level: advanced
233: .seealso: MatGetRow()
234: @*/
235: PetscErrorCode MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
236: {
242: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
243: if (!mat->ops->restorerow) return(0);
244: (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
245: return(0);
246: }
250: /*@C
251: MatGetRowUpperTriangular - Sets a flag to enable calls to MatGetRow() for matrix in MATSBAIJ format.
252: You should call MatRestoreRowUpperTriangular() after calling MatGetRow/MatRestoreRow() to disable the flag.
254: Not Collective
256: Input Parameters:
257: + mat - the matrix
259: Notes:
260: The flag is to ensure that users are aware of MatGetRow() only provides the upper trianglular part of the row for the matrices in MATSBAIJ format.
262: Level: advanced
264: Concepts: matrices^row access
266: .seealso: MatRestoreRowRowUpperTriangular()
267: @*/
269: PetscErrorCode MatGetRowUpperTriangular(Mat mat)
270: {
276: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
277: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
278: if (!mat->ops->getrowuppertriangular) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
279: MatPreallocated(mat);
280: (*mat->ops->getrowuppertriangular)(mat);
281: return(0);
282: }
286: /*@C
287: MatRestoreRowUpperTriangular - Disable calls to MatGetRow() for matrix in MATSBAIJ format.
289: Not Collective
291: Input Parameters:
292: + mat - the matrix
294: Notes:
295: This routine should be called after you have finished MatGetRow/MatRestoreRow().
298: Level: advanced
300: .seealso: MatGetRowUpperTriangular()
301: @*/
302: PetscErrorCode MatRestoreRowUpperTriangular(Mat mat)
303: {
308: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
309: if (!mat->ops->restorerowuppertriangular) return(0);
310: (*mat->ops->restorerowuppertriangular)(mat);
311: return(0);
312: }
316: /*@C
317: MatSetOptionsPrefix - Sets the prefix used for searching for all
318: Mat options in the database.
320: Collective on Mat
322: Input Parameter:
323: + A - the Mat context
324: - prefix - the prefix to prepend to all option names
326: Notes:
327: A hyphen (-) must NOT be given at the beginning of the prefix name.
328: The first character of all runtime options is AUTOMATICALLY the hyphen.
330: Level: advanced
332: .keywords: Mat, set, options, prefix, database
334: .seealso: MatSetFromOptions()
335: @*/
336: PetscErrorCode MatSetOptionsPrefix(Mat A,const char prefix[])
337: {
342: PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
343: return(0);
344: }
348: /*@C
349: MatAppendOptionsPrefix - Appends to the prefix used for searching for all
350: Mat options in the database.
352: Collective on Mat
354: Input Parameters:
355: + A - the Mat context
356: - prefix - the prefix to prepend to all option names
358: Notes:
359: A hyphen (-) must NOT be given at the beginning of the prefix name.
360: The first character of all runtime options is AUTOMATICALLY the hyphen.
362: Level: advanced
364: .keywords: Mat, append, options, prefix, database
366: .seealso: MatGetOptionsPrefix()
367: @*/
368: PetscErrorCode MatAppendOptionsPrefix(Mat A,const char prefix[])
369: {
371:
374: PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
375: return(0);
376: }
380: /*@C
381: MatGetOptionsPrefix - Sets the prefix used for searching for all
382: Mat options in the database.
384: Not Collective
386: Input Parameter:
387: . A - the Mat context
389: Output Parameter:
390: . prefix - pointer to the prefix string used
392: Notes: On the fortran side, the user should pass in a string 'prefix' of
393: sufficient length to hold the prefix.
395: Level: advanced
397: .keywords: Mat, get, options, prefix, database
399: .seealso: MatAppendOptionsPrefix()
400: @*/
401: PetscErrorCode MatGetOptionsPrefix(Mat A,const char *prefix[])
402: {
407: PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
408: return(0);
409: }
413: /*@
414: MatSetUp - Sets up the internal matrix data structures for the later use.
416: Collective on Mat
418: Input Parameters:
419: . A - the Mat context
421: Notes:
422: For basic use of the Mat classes the user need not explicitly call
423: MatSetUp(), since these actions will happen automatically.
425: Level: advanced
427: .keywords: Mat, setup
429: .seealso: MatCreate(), MatDestroy()
430: @*/
431: PetscErrorCode MatSetUp(Mat A)
432: {
437: MatSetUpPreallocation(A);
438: MatSetFromOptions(A);
439: return(0);
440: }
444: /*@C
445: MatView - Visualizes a matrix object.
447: Collective on Mat
449: Input Parameters:
450: + mat - the matrix
451: - viewer - visualization context
453: Notes:
454: The available visualization contexts include
455: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
456: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
457: output where only the first processor opens
458: the file. All other processors send their
459: data to the first processor to print.
460: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
462: The user can open alternative visualization contexts with
463: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
464: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
465: specified file; corresponding input uses MatLoad()
466: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
467: an X window display
468: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
469: Currently only the sequential dense and AIJ
470: matrix types support the Socket viewer.
472: The user can call PetscViewerSetFormat() to specify the output
473: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
474: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
475: + PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents
476: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
477: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
478: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
479: format common among all matrix types
480: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
481: format (which is in many cases the same as the default)
482: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
483: size and structure (not the matrix entries)
484: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
485: the matrix structure
487: Options Database Keys:
488: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
489: . -mat_view_info_detailed - Prints more detailed info
490: . -mat_view - Prints matrix in ASCII format
491: . -mat_view_matlab - Prints matrix in Matlab format
492: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
493: . -display <name> - Sets display name (default is host)
494: . -draw_pause <sec> - Sets number of seconds to pause after display
495: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
496: . -viewer_socket_machine <machine>
497: . -viewer_socket_port <port>
498: . -mat_view_binary - save matrix to file in binary format
499: - -viewer_binary_filename <name>
500: Level: beginner
502: Concepts: matrices^viewing
503: Concepts: matrices^plotting
504: Concepts: matrices^printing
506: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
507: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
508: @*/
509: PetscErrorCode MatView(Mat mat,PetscViewer viewer)
510: {
511: PetscErrorCode ierr;
512: PetscInt rows,cols;
513: PetscTruth iascii;
514: const char *cstr;
515: PetscViewerFormat format;
520: if (!viewer) {
521: PetscViewerASCIIGetStdout(mat->comm,&viewer);
522: }
525: if (!mat->assembled) SETERRQ(PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
526: MatPreallocated(mat);
528: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
529: if (iascii) {
530: PetscViewerGetFormat(viewer,&format);
531: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
532: if (mat->prefix) {
533: PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)\n",mat->prefix);
534: } else {
535: PetscViewerASCIIPrintf(viewer,"Matrix Object:\n");
536: }
537: PetscViewerASCIIPushTab(viewer);
538: MatGetType(mat,&cstr);
539: MatGetSize(mat,&rows,&cols);
540: PetscViewerASCIIPrintf(viewer,"type=%s, rows=%D, cols=%D\n",cstr,rows,cols);
541: if (mat->ops->getinfo) {
542: MatInfo info;
543: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
544: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",
545: (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
546: }
547: }
548: }
549: if (mat->ops->view) {
550: PetscViewerASCIIPushTab(viewer);
551: (*mat->ops->view)(mat,viewer);
552: PetscViewerASCIIPopTab(viewer);
553: } else if (!iascii) {
554: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
555: }
556: if (iascii) {
557: PetscViewerGetFormat(viewer,&format);
558: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
559: PetscViewerASCIIPopTab(viewer);
560: }
561: }
562: return(0);
563: }
567: /*@C
568: MatScaleSystem - Scale a vector solution and right hand side to
569: match the scaling of a scaled matrix.
570:
571: Collective on Mat
573: Input Parameter:
574: + mat - the matrix
575: . b - right hand side vector (or PETSC_NULL)
576: - x - solution vector (or PETSC_NULL)
579: Notes:
580: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
581: internally scaled, so this does nothing. For MPIROWBS it
582: permutes and diagonally scales.
584: The KSP methods automatically call this routine when required
585: (via PCPreSolve()) so it is rarely used directly.
587: Level: Developer
589: Concepts: matrices^scaling
591: .seealso: MatUseScaledForm(), MatUnScaleSystem()
592: @*/
593: PetscErrorCode MatScaleSystem(Mat mat,Vec b,Vec x)
594: {
600: MatPreallocated(mat);
604: if (mat->ops->scalesystem) {
605: (*mat->ops->scalesystem)(mat,b,x);
606: }
607: return(0);
608: }
612: /*@C
613: MatUnScaleSystem - Unscales a vector solution and right hand side to
614: match the original scaling of a scaled matrix.
615:
616: Collective on Mat
618: Input Parameter:
619: + mat - the matrix
620: . b - right hand side vector (or PETSC_NULL)
621: - x - solution vector (or PETSC_NULL)
624: Notes:
625: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
626: internally scaled, so this does nothing. For MPIROWBS it
627: permutes and diagonally scales.
629: The KSP methods automatically call this routine when required
630: (via PCPreSolve()) so it is rarely used directly.
632: Level: Developer
634: .seealso: MatUseScaledForm(), MatScaleSystem()
635: @*/
636: PetscErrorCode MatUnScaleSystem(Mat mat,Vec b,Vec x)
637: {
643: MatPreallocated(mat);
646: if (mat->ops->unscalesystem) {
647: (*mat->ops->unscalesystem)(mat,b,x);
648: }
649: return(0);
650: }
654: /*@C
655: MatUseScaledForm - For matrix storage formats that scale the
656: matrix (for example MPIRowBS matrices are diagonally scaled on
657: assembly) indicates matrix operations (MatMult() etc) are
658: applied using the scaled matrix.
659:
660: Collective on Mat
662: Input Parameter:
663: + mat - the matrix
664: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
665: applying the original matrix
667: Notes:
668: For scaled matrix formats, applying the original, unscaled matrix
669: will be slightly more expensive
671: Level: Developer
673: .seealso: MatScaleSystem(), MatUnScaleSystem()
674: @*/
675: PetscErrorCode MatUseScaledForm(Mat mat,PetscTruth scaled)
676: {
682: MatPreallocated(mat);
683: if (mat->ops->usescaledform) {
684: (*mat->ops->usescaledform)(mat,scaled);
685: }
686: return(0);
687: }
691: /*@
692: MatDestroy - Frees space taken by a matrix.
693:
694: Collective on Mat
696: Input Parameter:
697: . A - the matrix
699: Level: beginner
701: @*/
702: PetscErrorCode MatDestroy(Mat A)
703: {
707: if (--A->refct > 0) return(0);
708: MatPreallocated(A);
709: /* if memory was published with AMS then destroy it */
710: PetscObjectDepublish(A);
711: if (A->ops->destroy) {
712: (*A->ops->destroy)(A);
713: }
714: if (A->mapping) {
715: ISLocalToGlobalMappingDestroy(A->mapping);
716: }
717: if (A->bmapping) {
718: ISLocalToGlobalMappingDestroy(A->bmapping);
719: }
720: PetscFree(A->rmap.range);
721: PetscFree(A->cmap.range);
722: if (A->spptr){PetscFree(A->spptr);}
723: PetscHeaderDestroy(A);
724: return(0);
725: }
729: /*@
730: MatValid - Checks whether a matrix object is valid.
732: Collective on Mat
734: Input Parameter:
735: . m - the matrix to check
737: Output Parameter:
738: flg - flag indicating matrix status, either
739: PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
741: Level: developer
743: Concepts: matrices^validity
744: @*/
745: PetscErrorCode MatValid(Mat m,PetscTruth *flg)
746: {
749: if (!m) *flg = PETSC_FALSE;
750: else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
751: else *flg = PETSC_TRUE;
752: return(0);
753: }
757: /*@
758: MatSetValues - Inserts or adds a block of values into a matrix.
759: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
760: MUST be called after all calls to MatSetValues() have been completed.
762: Not Collective
764: Input Parameters:
765: + mat - the matrix
766: . v - a logically two-dimensional array of values
767: . m, idxm - the number of rows and their global indices
768: . n, idxn - the number of columns and their global indices
769: - addv - either ADD_VALUES or INSERT_VALUES, where
770: ADD_VALUES adds values to any existing entries, and
771: INSERT_VALUES replaces existing entries with new values
773: Notes:
774: By default the values, v, are row-oriented and unsorted.
775: See MatSetOption() for other options.
777: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
778: options cannot be mixed without intervening calls to the assembly
779: routines.
781: MatSetValues() uses 0-based row and column numbers in Fortran
782: as well as in C.
784: Negative indices may be passed in idxm and idxn, these rows and columns are
785: simply ignored. This allows easily inserting element stiffness matrices
786: with homogeneous Dirchlet boundary conditions that you don't want represented
787: in the matrix.
789: Efficiency Alert:
790: The routine MatSetValuesBlocked() may offer much better efficiency
791: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
793: Level: beginner
795: Concepts: matrices^putting entries in
797: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
798: InsertMode, INSERT_VALUES, ADD_VALUES
799: @*/
800: PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
801: {
805: if (!m || !n) return(0); /* no values to insert */
811: MatPreallocated(mat);
812: if (mat->insertmode == NOT_SET_VALUES) {
813: mat->insertmode = addv;
814: }
815: #if defined(PETSC_USE_DEBUG)
816: else if (mat->insertmode != addv) {
817: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
818: }
819: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
820: #endif
822: if (mat->assembled) {
823: mat->was_assembled = PETSC_TRUE;
824: mat->assembled = PETSC_FALSE;
825: }
827: if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
828: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
830: return(0);
831: }
836: /*@
837: MatSetValuesRowLocal - Inserts a row (block row for BAIJ matrices) of nonzero
838: values into a matrix
840: Not Collective
842: Input Parameters:
843: + mat - the matrix
844: . row - the (block) row to set
845: - v - a logically two-dimensional array of values
847: Notes:
848: By the values, v, are column-oriented (for the block version) and sorted
850: All the nonzeros in the row must be provided
852: The matrix must have previously had its column indices set
854: The row must belong to this process
856: Level: intermediate
858: Concepts: matrices^putting entries in
860: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
861: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues(), MatSetValuesRow(), MatSetLocalToGlobalMapping()
862: @*/
863: PetscErrorCode MatSetValuesRowLocal(Mat mat,PetscInt row,const PetscScalar v[])
864: {
871: MatSetValuesRow(mat, mat->mapping->indices[row],v);
872: return(0);
873: }
877: /*@
878: MatSetValuesRow - Inserts a row (block row for BAIJ matrices) of nonzero
879: values into a matrix
881: Not Collective
883: Input Parameters:
884: + mat - the matrix
885: . row - the (block) row to set
886: - v - a logically two-dimensional array of values
888: Notes:
889: By the values, v, are column-oriented (for the block version) and sorted
891: All the nonzeros in the row must be provided
893: The matrix must have previously had its column indices set
895: The row must belong to this process
897: Level: intermediate
899: Concepts: matrices^putting entries in
901: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
902: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
903: @*/
904: PetscErrorCode MatSetValuesRow(Mat mat,PetscInt row,const PetscScalar v[])
905: {
912: #if defined(PETSC_USE_DEBUG)
913: if (mat->insertmode == ADD_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add and insert values");
914: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
915: #endif
916: mat->insertmode = INSERT_VALUES;
918: if (mat->assembled) {
919: mat->was_assembled = PETSC_TRUE;
920: mat->assembled = PETSC_FALSE;
921: }
923: if (!mat->ops->setvaluesrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
924: (*mat->ops->setvaluesrow)(mat,row,v);
926: return(0);
927: }
931: /*@
932: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
933: Using structured grid indexing
935: Not Collective
937: Input Parameters:
938: + mat - the matrix
939: . v - a logically two-dimensional array of values
940: . m - number of rows being entered
941: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
942: . n - number of columns being entered
943: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
944: - addv - either ADD_VALUES or INSERT_VALUES, where
945: ADD_VALUES adds values to any existing entries, and
946: INSERT_VALUES replaces existing entries with new values
948: Notes:
949: By default the values, v, are row-oriented and unsorted.
950: See MatSetOption() for other options.
952: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
953: options cannot be mixed without intervening calls to the assembly
954: routines.
956: The grid coordinates are across the entire grid, not just the local portion
958: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
959: as well as in C.
961: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
963: In order to use this routine you must either obtain the matrix with DAGetMatrix()
964: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
966: The columns and rows in the stencil passed in MUST be contained within the
967: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
968: if you create a DA with an overlap of one grid level and on a particular process its first
969: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
970: first i index you can use in your column and row indices in MatSetStencil() is 5.
972: In Fortran idxm and idxn should be declared as
973: $ MatStencil idxm(4,m),idxn(4,n)
974: and the values inserted using
975: $ idxm(MatStencil_i,1) = i
976: $ idxm(MatStencil_j,1) = j
977: $ idxm(MatStencil_k,1) = k
978: $ idxm(MatStencil_c,1) = c
979: etc
981: For periodic boundary conditions use negative indices for values to the left (below 0; that are to be
982: obtained by wrapping values from right edge). For values to the right of the last entry using that index plus one
983: etc to obtain values that obtained by wrapping the values from the left edge.
985: For indices that don't mean anything for your case (like the k index when working in 2d) or the c index when you have
986: a single value per point) you can skip filling those indices.
988: Inspired by the structured grid interface to the HYPRE package
989: (http://www.llnl.gov/CASC/hypre)
991: Efficiency Alert:
992: The routine MatSetValuesBlockedStencil() may offer much better efficiency
993: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
995: Level: beginner
997: Concepts: matrices^putting entries in
999: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1000: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1001: @*/
1002: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1003: {
1005: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1006: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1009: if (!m || !n) return(0); /* no values to insert */
1016: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1017: if (n > 256) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1019: for (i=0; i<m; i++) {
1020: for (j=0; j<3-sdim; j++) dxm++;
1021: tmp = *dxm++ - starts[0];
1022: for (j=0; j<dim-1; j++) {
1023: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1024: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1025: }
1026: if (mat->stencil.noc) dxm++;
1027: jdxm[i] = tmp;
1028: }
1029: for (i=0; i<n; i++) {
1030: for (j=0; j<3-sdim; j++) dxn++;
1031: tmp = *dxn++ - starts[0];
1032: for (j=0; j<dim-1; j++) {
1033: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1034: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1035: }
1036: if (mat->stencil.noc) dxn++;
1037: jdxn[i] = tmp;
1038: }
1039: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1040: return(0);
1041: }
1045: /*@C
1046: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1047: Using structured grid indexing
1049: Not Collective
1051: Input Parameters:
1052: + mat - the matrix
1053: . v - a logically two-dimensional array of values
1054: . m - number of rows being entered
1055: . idxm - grid coordinates for matrix rows being entered
1056: . n - number of columns being entered
1057: . idxn - grid coordinates for matrix columns being entered
1058: - addv - either ADD_VALUES or INSERT_VALUES, where
1059: ADD_VALUES adds values to any existing entries, and
1060: INSERT_VALUES replaces existing entries with new values
1062: Notes:
1063: By default the values, v, are row-oriented and unsorted.
1064: See MatSetOption() for other options.
1066: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1067: options cannot be mixed without intervening calls to the assembly
1068: routines.
1070: The grid coordinates are across the entire grid, not just the local portion
1072: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1073: as well as in C.
1075: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
1077: In order to use this routine you must either obtain the matrix with DAGetMatrix()
1078: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1080: The columns and rows in the stencil passed in MUST be contained within the
1081: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
1082: if you create a DA with an overlap of one grid level and on a particular process its first
1083: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1084: first i index you can use in your column and row indices in MatSetStencil() is 5.
1086: In Fortran idxm and idxn should be declared as
1087: $ MatStencil idxm(4,m),idxn(4,n)
1088: and the values inserted using
1089: $ idxm(MatStencil_i,1) = i
1090: $ idxm(MatStencil_j,1) = j
1091: $ idxm(MatStencil_k,1) = k
1092: etc
1094: Negative indices may be passed in idxm and idxn, these rows and columns are
1095: simply ignored. This allows easily inserting element stiffness matrices
1096: with homogeneous Dirchlet boundary conditions that you don't want represented
1097: in the matrix.
1099: Inspired by the structured grid interface to the HYPRE package
1100: (http://www.llnl.gov/CASC/hypre)
1102: Level: beginner
1104: Concepts: matrices^putting entries in
1106: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1107: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1108: @*/
1109: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1110: {
1112: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1113: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1116: if (!m || !n) return(0); /* no values to insert */
1123: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1124: if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1126: for (i=0; i<m; i++) {
1127: for (j=0; j<3-sdim; j++) dxm++;
1128: tmp = *dxm++ - starts[0];
1129: for (j=0; j<sdim-1; j++) {
1130: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1131: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1132: }
1133: dxm++;
1134: jdxm[i] = tmp;
1135: }
1136: for (i=0; i<n; i++) {
1137: for (j=0; j<3-sdim; j++) dxn++;
1138: tmp = *dxn++ - starts[0];
1139: for (j=0; j<sdim-1; j++) {
1140: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1141: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1142: }
1143: dxn++;
1144: jdxn[i] = tmp;
1145: }
1146: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1147: return(0);
1148: }
1152: /*@
1153: MatSetStencil - Sets the grid information for setting values into a matrix via
1154: MatSetValuesStencil()
1156: Not Collective
1158: Input Parameters:
1159: + mat - the matrix
1160: . dim - dimension of the grid 1, 2, or 3
1161: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1162: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1163: - dof - number of degrees of freedom per node
1166: Inspired by the structured grid interface to the HYPRE package
1167: (www.llnl.gov/CASC/hyper)
1169: For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the
1170: user.
1171:
1172: Level: beginner
1174: Concepts: matrices^putting entries in
1176: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1177: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1178: @*/
1179: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1180: {
1181: PetscInt i;
1188: mat->stencil.dim = dim + (dof > 1);
1189: for (i=0; i<dim; i++) {
1190: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1191: mat->stencil.starts[i] = starts[dim-i-1];
1192: }
1193: mat->stencil.dims[dim] = dof;
1194: mat->stencil.starts[dim] = 0;
1195: mat->stencil.noc = (PetscTruth)(dof == 1);
1196: return(0);
1197: }
1201: /*@
1202: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1204: Not Collective
1206: Input Parameters:
1207: + mat - the matrix
1208: . v - a logically two-dimensional array of values
1209: . m, idxm - the number of block rows and their global block indices
1210: . n, idxn - the number of block columns and their global block indices
1211: - addv - either ADD_VALUES or INSERT_VALUES, where
1212: ADD_VALUES adds values to any existing entries, and
1213: INSERT_VALUES replaces existing entries with new values
1215: Notes:
1216: The m and n count the NUMBER of blocks in the row direction and column direction,
1217: NOT the total number of rows/columns; for example, if the block size is 2 and
1218: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1220: By default the values, v, are row-oriented and unsorted. So the layout of
1221: v is the same as for MatSetValues(). See MatSetOption() for other options.
1223: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1224: options cannot be mixed without intervening calls to the assembly
1225: routines.
1227: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1228: as well as in C.
1230: Negative indices may be passed in idxm and idxn, these rows and columns are
1231: simply ignored. This allows easily inserting element stiffness matrices
1232: with homogeneous Dirchlet boundary conditions that you don't want represented
1233: in the matrix.
1235: Each time an entry is set within a sparse matrix via MatSetValues(),
1236: internal searching must be done to determine where to place the the
1237: data in the matrix storage space. By instead inserting blocks of
1238: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1239: reduced.
1241: Example:
1242: $ Suppose m=n=2 and block size(bs) = 2 The matrix is
1243: $
1244: $ 1 2 | 3 4
1245: $ 5 6 | 7 8
1246: $ - - - | - - -
1247: $ 9 10 | 11 12
1248: $ 13 14 | 15 16
1249: $
1250: $ v[] should be passed in like
1251: $ v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
1253: Restrictions:
1254: MatSetValuesBlocked() is currently supported only for the BAIJ and SBAIJ formats
1256: Level: intermediate
1258: Concepts: matrices^putting entries in blocked
1260: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1261: @*/
1262: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1263: {
1267: if (!m || !n) return(0); /* no values to insert */
1273: MatPreallocated(mat);
1274: if (mat->insertmode == NOT_SET_VALUES) {
1275: mat->insertmode = addv;
1276: }
1277: #if defined(PETSC_USE_DEBUG)
1278: else if (mat->insertmode != addv) {
1279: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1280: }
1281: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1282: #endif
1284: if (mat->assembled) {
1285: mat->was_assembled = PETSC_TRUE;
1286: mat->assembled = PETSC_FALSE;
1287: }
1289: if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1290: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1292: return(0);
1293: }
1297: /*@
1298: MatGetValues - Gets a block of values from a matrix.
1300: Not Collective; currently only returns a local block
1302: Input Parameters:
1303: + mat - the matrix
1304: . v - a logically two-dimensional array for storing the values
1305: . m, idxm - the number of rows and their global indices
1306: - n, idxn - the number of columns and their global indices
1308: Notes:
1309: The user must allocate space (m*n PetscScalars) for the values, v.
1310: The values, v, are then returned in a row-oriented format,
1311: analogous to that used by default in MatSetValues().
1313: MatGetValues() uses 0-based row and column numbers in
1314: Fortran as well as in C.
1316: MatGetValues() requires that the matrix has been assembled
1317: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1318: MatSetValues() and MatGetValues() CANNOT be made in succession
1319: without intermediate matrix assembly.
1321: Level: advanced
1323: Concepts: matrices^accessing values
1325: .seealso: MatGetRow(), MatGetSubMatrices(), MatSetValues()
1326: @*/
1327: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1328: {
1337: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1338: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1339: if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1340: MatPreallocated(mat);
1343: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1345: return(0);
1346: }
1350: /*@
1351: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1352: the routine MatSetValuesLocal() to allow users to insert matrix entries
1353: using a local (per-processor) numbering.
1355: Not Collective
1357: Input Parameters:
1358: + x - the matrix
1359: - mapping - mapping created with ISLocalToGlobalMappingCreate()
1360: or ISLocalToGlobalMappingCreateIS()
1362: Level: intermediate
1364: Concepts: matrices^local to global mapping
1365: Concepts: local to global mapping^for matrices
1367: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1368: @*/
1369: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
1370: {
1376: if (x->mapping) {
1377: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1378: }
1379: MatPreallocated(x);
1381: if (x->ops->setlocaltoglobalmapping) {
1382: (*x->ops->setlocaltoglobalmapping)(x,mapping);
1383: } else {
1384: