Actual source code: matrix.c
2: /*
3: This is where the abstract matrix operations are defined
4: */
6: #include <private/matimpl.h> /*I "petscmat.h" I*/
7: #include <private/vecimpl.h>
9: /* Logging support */
10: PetscClassId MAT_CLASSID;
11: PetscClassId MAT_FDCOLORING_CLASSID;
13: PetscLogEvent MAT_Mult, MAT_Mults, MAT_MultConstrained, MAT_MultAdd, MAT_MultTranspose;
14: PetscLogEvent MAT_MultTransposeConstrained, MAT_MultTransposeAdd, MAT_Solve, MAT_Solves, MAT_SolveAdd, MAT_SolveTranspose, MAT_MatSolve;
15: PetscLogEvent MAT_SolveTransposeAdd, MAT_SOR, MAT_ForwardSolve, MAT_BackwardSolve, MAT_LUFactor, MAT_LUFactorSymbolic;
16: PetscLogEvent MAT_LUFactorNumeric, MAT_CholeskyFactor, MAT_CholeskyFactorSymbolic, MAT_CholeskyFactorNumeric, MAT_ILUFactor;
17: PetscLogEvent MAT_ILUFactorSymbolic, MAT_ICCFactorSymbolic, MAT_Copy, MAT_Convert, MAT_Scale, MAT_AssemblyBegin;
18: PetscLogEvent MAT_AssemblyEnd, MAT_SetValues, MAT_GetValues, MAT_GetRow, MAT_GetRowIJ, MAT_GetSubMatrices, MAT_GetColoring, MAT_GetOrdering, MAT_GetRedundantMatrix, MAT_GetSeqNonzeroStructure;
19: PetscLogEvent MAT_IncreaseOverlap, MAT_Partitioning, MAT_ZeroEntries, MAT_Load, MAT_View, MAT_AXPY, MAT_FDColoringCreate;
20: PetscLogEvent MAT_FDColoringApply,MAT_Transpose,MAT_FDColoringFunction;
21: PetscLogEvent MAT_MatMult, MAT_MatMultSymbolic, MAT_MatMultNumeric;
22: PetscLogEvent MAT_PtAP, MAT_PtAPSymbolic, MAT_PtAPNumeric;
23: PetscLogEvent MAT_MatMultTranspose, MAT_MatMultTransposeSymbolic, MAT_MatMultTransposeNumeric;
24: PetscLogEvent MAT_MultHermitianTranspose,MAT_MultHermitianTransposeAdd;
25: PetscLogEvent MAT_Getsymtranspose, MAT_Getsymtransreduced, MAT_Transpose_SeqAIJ, MAT_GetBrowsOfAcols;
26: PetscLogEvent MAT_GetBrowsOfAocols, MAT_Getlocalmat, MAT_Getlocalmatcondensed, MAT_Seqstompi, MAT_Seqstompinum, MAT_Seqstompisym;
27: PetscLogEvent MAT_Applypapt, MAT_Applypapt_numeric, MAT_Applypapt_symbolic, MAT_GetSequentialNonzeroStructure;
28: PetscLogEvent MAT_GetMultiProcBlock;
29: PetscLogEvent MAT_CUSPCopyToGPU, MAT_SetValuesBatch, MAT_SetValuesBatchI, MAT_SetValuesBatchII, MAT_SetValuesBatchIII, MAT_SetValuesBatchIV;
31: /* nasty global values for MatSetValue() */
32: PetscInt MatSetValue_Row = 0;
33: PetscInt MatSetValue_Column = 0;
34: PetscScalar MatSetValue_Value = 0.0;
36: const char *const MatFactorTypes[] = {"NONE","LU","CHOLESKY","ILU","ICC","ILUDT","MatFactorType","MAT_FACTOR_",0};
40: /*@C
41: MatFindNonzeroRows - Locate all rows that are not completely zero in the matrix
43: Input Parameter:
44: . A - the matrix
46: Output Parameter:
47: . keptrows - the rows that are not completely zero
49: Level: intermediate
51: @*/
52: PetscErrorCode MatFindNonzeroRows(Mat mat,IS *keptrows)
53: {
54: PetscErrorCode ierr;
58: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
59: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
60: if (!mat->ops->findnonzerorows) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Not coded for this matrix type");
61: (*mat->ops->findnonzerorows)(mat,keptrows);
62: return(0);
63: }
67: /*@
68: MatGetDiagonalBlock - Returns the part of the matrix associated with the on-process coupling
70: Not Collective
72: Input Parameters:
73: . A - the matrix
75: Output Parameters:
76: . a - the diagonal part (which is a SEQUENTIAL matrix)
78: Notes: see the manual page for MatCreateMPIAIJ() for more information on the "diagonal part" of the matrix.
80: Level: advanced
82: @*/
83: PetscErrorCode MatGetDiagonalBlock(Mat A,Mat *a)
84: {
85: PetscErrorCode ierr,(*f)(Mat,Mat*);
86: PetscMPIInt size;
92: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
93: if (A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
94: MPI_Comm_size(((PetscObject)A)->comm,&size);
95: PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);
96: if (f) {
97: (*f)(A,a);
98: return(0);
99: } else if (size == 1) {
100: *a = A;
101: } else {
102: const MatType mattype;
103: MatGetType(A,&mattype);
104: SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Matrix type %s does not support getting diagonal block",mattype);
105: }
106: return(0);
107: }
111: /*@
112: MatGetTrace - Gets the trace of a matrix. The sum of the diagonal entries.
114: Collective on Mat
116: Input Parameters:
117: . mat - the matrix
119: Output Parameter:
120: . trace - the sum of the diagonal entries
122: Level: advanced
124: @*/
125: PetscErrorCode MatGetTrace(Mat mat,PetscScalar *trace)
126: {
128: Vec diag;
131: MatGetVecs(mat,&diag,PETSC_NULL);
132: MatGetDiagonal(mat,diag);
133: VecSum(diag,trace);
134: VecDestroy(&diag);
135: return(0);
136: }
140: /*@
141: MatRealPart - Zeros out the imaginary part of the matrix
143: Logically Collective on Mat
145: Input Parameters:
146: . mat - the matrix
148: Level: advanced
151: .seealso: MatImaginaryPart()
152: @*/
153: PetscErrorCode MatRealPart(Mat mat)
154: {
160: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
161: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
162: if (!mat->ops->realpart) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
163: MatPreallocated(mat);
164: (*mat->ops->realpart)(mat);
165: #if defined(PETSC_HAVE_CUSP)
166: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
167: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
168: }
169: #endif
170: return(0);
171: }
175: /*@C
176: MatGetGhosts - Get the global index of all ghost nodes defined by the sparse matrix
178: Collective on Mat
180: Input Parameter:
181: . mat - the matrix
183: Output Parameters:
184: + nghosts - number of ghosts (note for BAIJ matrices there is one ghost for each block)
185: - ghosts - the global indices of the ghost points
187: Notes: the nghosts and ghosts are suitable to pass into VecCreateGhost()
189: Level: advanced
191: @*/
192: PetscErrorCode MatGetGhosts(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
193: {
199: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
200: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
201: if (!mat->ops->getghosts) {
202: if (nghosts) *nghosts = 0;
203: if (ghosts) *ghosts = 0;
204: } else {
205: (*mat->ops->getghosts)(mat,nghosts,ghosts);
206: }
207: return(0);
208: }
213: /*@
214: MatImaginaryPart - Moves the imaginary part of the matrix to the real part and zeros the imaginary part
216: Logically Collective on Mat
218: Input Parameters:
219: . mat - the matrix
221: Level: advanced
224: .seealso: MatRealPart()
225: @*/
226: PetscErrorCode MatImaginaryPart(Mat mat)
227: {
233: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
234: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
235: if (!mat->ops->imaginarypart) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
236: MatPreallocated(mat);
237: (*mat->ops->imaginarypart)(mat);
238: #if defined(PETSC_HAVE_CUSP)
239: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
240: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
241: }
242: #endif
243: return(0);
244: }
248: /*@
249: MatMissingDiagonal - Determine if sparse matrix is missing a diagonal entry (or block entry for BAIJ matrices)
251: Collective on Mat
253: Input Parameter:
254: . mat - the matrix
256: Output Parameters:
257: + missing - is any diagonal missing
258: - dd - first diagonal entry that is missing (optional)
260: Level: advanced
263: .seealso: MatRealPart()
264: @*/
265: PetscErrorCode MatMissingDiagonal(Mat mat,PetscBool *missing,PetscInt *dd)
266: {
272: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
273: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
274: if (!mat->ops->missingdiagonal) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
275: (*mat->ops->missingdiagonal)(mat,missing,dd);
276: return(0);
277: }
281: /*@C
282: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
283: for each row that you get to ensure that your application does
284: not bleed memory.
286: Not Collective
288: Input Parameters:
289: + mat - the matrix
290: - row - the row to get
292: Output Parameters:
293: + ncols - if not NULL, the number of nonzeros in the row
294: . cols - if not NULL, the column numbers
295: - vals - if not NULL, the values
297: Notes:
298: This routine is provided for people who need to have direct access
299: to the structure of a matrix. We hope that we provide enough
300: high-level matrix routines that few users will need it.
302: MatGetRow() always returns 0-based column indices, regardless of
303: whether the internal representation is 0-based (default) or 1-based.
305: For better efficiency, set cols and/or vals to PETSC_NULL if you do
306: not wish to extract these quantities.
308: The user can only examine the values extracted with MatGetRow();
309: the values cannot be altered. To change the matrix entries, one
310: must use MatSetValues().
312: You can only have one call to MatGetRow() outstanding for a particular
313: matrix at a time, per processor. MatGetRow() can only obtain rows
314: associated with the given processor, it cannot get rows from the
315: other processors; for that we suggest using MatGetSubMatrices(), then
316: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
317: is in the global number of rows.
319: Fortran Notes:
320: The calling sequence from Fortran is
321: .vb
322: MatGetRow(matrix,row,ncols,cols,values,ierr)
323: Mat matrix (input)
324: integer row (input)
325: integer ncols (output)
326: integer cols(maxcols) (output)
327: double precision (or double complex) values(maxcols) output
328: .ve
329: where maxcols >= maximum nonzeros in any row of the matrix.
332: Caution:
333: Do not try to change the contents of the output arrays (cols and vals).
334: In some cases, this may corrupt the matrix.
336: Level: advanced
338: Concepts: matrices^row access
340: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubMatrices(), MatGetDiagonal()
341: @*/
342: PetscErrorCode MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
343: {
345: PetscInt incols;
350: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
351: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
352: if (!mat->ops->getrow) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
353: MatPreallocated(mat);
354: PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
355: (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
356: if (ncols) *ncols = incols;
357: PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
358: return(0);
359: }
363: /*@
364: MatConjugate - replaces the matrix values with their complex conjugates
366: Logically Collective on Mat
368: Input Parameters:
369: . mat - the matrix
371: Level: advanced
373: .seealso: VecConjugate()
374: @*/
375: PetscErrorCode MatConjugate(Mat mat)
376: {
381: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
382: if (!mat->ops->conjugate) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Not provided for this matrix format, send email to petsc-maint@mcs.anl.gov");
383: (*mat->ops->conjugate)(mat);
384: #if defined(PETSC_HAVE_CUSP)
385: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
386: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
387: }
388: #endif
389: return(0);
390: }
394: /*@C
395: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
397: Not Collective
399: Input Parameters:
400: + mat - the matrix
401: . row - the row to get
402: . ncols, cols - the number of nonzeros and their columns
403: - vals - if nonzero the column values
405: Notes:
406: This routine should be called after you have finished examining the entries.
408: Fortran Notes:
409: The calling sequence from Fortran is
410: .vb
411: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
412: Mat matrix (input)
413: integer row (input)
414: integer ncols (output)
415: integer cols(maxcols) (output)
416: double precision (or double complex) values(maxcols) output
417: .ve
418: Where maxcols >= maximum nonzeros in any row of the matrix.
420: In Fortran MatRestoreRow() MUST be called after MatGetRow()
421: before another call to MatGetRow() can be made.
423: Level: advanced
425: .seealso: MatGetRow()
426: @*/
427: PetscErrorCode MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
428: {
434: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
435: if (!mat->ops->restorerow) return(0);
436: (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
437: return(0);
438: }
442: /*@
443: MatGetRowUpperTriangular - Sets a flag to enable calls to MatGetRow() for matrix in MATSBAIJ format.
444: You should call MatRestoreRowUpperTriangular() after calling MatGetRow/MatRestoreRow() to disable the flag.
446: Not Collective
448: Input Parameters:
449: + mat - the matrix
451: Notes:
452: 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.
454: Level: advanced
456: Concepts: matrices^row access
458: .seealso: MatRestoreRowRowUpperTriangular()
459: @*/
460: PetscErrorCode MatGetRowUpperTriangular(Mat mat)
461: {
467: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
468: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
469: if (!mat->ops->getrowuppertriangular) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
470: MatPreallocated(mat);
471: (*mat->ops->getrowuppertriangular)(mat);
472: return(0);
473: }
477: /*@
478: MatRestoreRowUpperTriangular - Disable calls to MatGetRow() for matrix in MATSBAIJ format.
480: Not Collective
482: Input Parameters:
483: + mat - the matrix
485: Notes:
486: This routine should be called after you have finished MatGetRow/MatRestoreRow().
489: Level: advanced
491: .seealso: MatGetRowUpperTriangular()
492: @*/
493: PetscErrorCode MatRestoreRowUpperTriangular(Mat mat)
494: {
499: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
500: if (!mat->ops->restorerowuppertriangular) return(0);
501: (*mat->ops->restorerowuppertriangular)(mat);
502: return(0);
503: }
507: /*@C
508: MatSetOptionsPrefix - Sets the prefix used for searching for all
509: Mat options in the database.
511: Logically Collective on Mat
513: Input Parameter:
514: + A - the Mat context
515: - prefix - the prefix to prepend to all option names
517: Notes:
518: A hyphen (-) must NOT be given at the beginning of the prefix name.
519: The first character of all runtime options is AUTOMATICALLY the hyphen.
521: Level: advanced
523: .keywords: Mat, set, options, prefix, database
525: .seealso: MatSetFromOptions()
526: @*/
527: PetscErrorCode MatSetOptionsPrefix(Mat A,const char prefix[])
528: {
533: PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
534: return(0);
535: }
539: /*@C
540: MatAppendOptionsPrefix - Appends to the prefix used for searching for all
541: Mat options in the database.
543: Logically Collective on Mat
545: Input Parameters:
546: + A - the Mat context
547: - prefix - the prefix to prepend to all option names
549: Notes:
550: A hyphen (-) must NOT be given at the beginning of the prefix name.
551: The first character of all runtime options is AUTOMATICALLY the hyphen.
553: Level: advanced
555: .keywords: Mat, append, options, prefix, database
557: .seealso: MatGetOptionsPrefix()
558: @*/
559: PetscErrorCode MatAppendOptionsPrefix(Mat A,const char prefix[])
560: {
562:
565: PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
566: return(0);
567: }
571: /*@C
572: MatGetOptionsPrefix - Sets the prefix used for searching for all
573: Mat options in the database.
575: Not Collective
577: Input Parameter:
578: . A - the Mat context
580: Output Parameter:
581: . prefix - pointer to the prefix string used
583: Notes: On the fortran side, the user should pass in a string 'prefix' of
584: sufficient length to hold the prefix.
586: Level: advanced
588: .keywords: Mat, get, options, prefix, database
590: .seealso: MatAppendOptionsPrefix()
591: @*/
592: PetscErrorCode MatGetOptionsPrefix(Mat A,const char *prefix[])
593: {
598: PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
599: return(0);
600: }
604: /*@
605: MatSetUp - Sets up the internal matrix data structures for the later use.
607: Collective on Mat
609: Input Parameters:
610: . A - the Mat context
612: Notes:
613: For basic use of the Mat classes the user need not explicitly call
614: MatSetUp(), since these actions will happen automatically.
616: Level: advanced
618: .keywords: Mat, setup
620: .seealso: MatCreate(), MatDestroy()
621: @*/
622: PetscErrorCode MatSetUp(Mat A)
623: {
624: PetscMPIInt size;
629: if (!((PetscObject)A)->type_name) {
630: MPI_Comm_size(((PetscObject)A)->comm, &size);
631: if (size == 1) {
632: MatSetType(A, MATSEQAIJ);
633: } else {
634: MatSetType(A, MATMPIAIJ);
635: }
636: }
637: MatSetUpPreallocation(A);
638: return(0);
639: }
644: /*@C
645: MatView - Visualizes a matrix object.
647: Collective on Mat
649: Input Parameters:
650: + mat - the matrix
651: - viewer - visualization context
653: Notes:
654: The available visualization contexts include
655: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
656: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
657: output where only the first processor opens
658: the file. All other processors send their
659: data to the first processor to print.
660: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
662: The user can open alternative visualization contexts with
663: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
664: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
665: specified file; corresponding input uses MatLoad()
666: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
667: an X window display
668: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
669: Currently only the sequential dense and AIJ
670: matrix types support the Socket viewer.
672: The user can call PetscViewerSetFormat() to specify the output
673: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
674: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
675: + PETSC_VIEWER_DEFAULT - default, prints matrix contents
676: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
677: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
678: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
679: format common among all matrix types
680: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
681: format (which is in many cases the same as the default)
682: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
683: size and structure (not the matrix entries)
684: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
685: the matrix structure
687: Options Database Keys:
688: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
689: . -mat_view_info_detailed - Prints more detailed info
690: . -mat_view - Prints matrix in ASCII format
691: . -mat_view_matlab - Prints matrix in Matlab format
692: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
693: . -display <name> - Sets display name (default is host)
694: . -draw_pause <sec> - Sets number of seconds to pause after display
695: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see the <a href="../../docs/manual.pdf">users manual</a> for details).
696: . -viewer_socket_machine <machine>
697: . -viewer_socket_port <port>
698: . -mat_view_binary - save matrix to file in binary format
699: - -viewer_binary_filename <name>
700: Level: beginner
702: Notes: see the manual page for MatLoad() for the exact format of the binary file when the binary
703: viewer is used.
705: See bin/matlab/PetscBinaryRead.m for a Matlab code that can read in the binary file when the binary
706: viewer is used.
708: Concepts: matrices^viewing
709: Concepts: matrices^plotting
710: Concepts: matrices^printing
712: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
713: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
714: @*/
715: PetscErrorCode MatView(Mat mat,PetscViewer viewer)
716: {
717: PetscErrorCode ierr;
718: PetscInt rows,cols;
719: PetscBool iascii;
720: PetscViewerFormat format;
725: if (!viewer) {
726: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
727: }
730: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
731: MatPreallocated(mat);
733: PetscLogEventBegin(MAT_View,mat,viewer,0,0);
734: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
735: if (iascii) {
736: PetscViewerGetFormat(viewer,&format);
737: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
738: PetscObjectPrintClassNamePrefixType((PetscObject)mat,viewer,"Matrix Object");
739: PetscViewerASCIIPushTab(viewer);
740: MatGetSize(mat,&rows,&cols);
741: PetscViewerASCIIPrintf(viewer,"rows=%D, cols=%D\n",rows,cols);
742: if (mat->factortype) {
743: const MatSolverPackage solver;
744: MatFactorGetSolverPackage(mat,&solver);
745: PetscViewerASCIIPrintf(viewer,"package used to perform factorization: %s\n",solver);
746: }
747: if (mat->ops->getinfo) {
748: MatInfo info;
749: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
750: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",(PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
751: PetscViewerASCIIPrintf(viewer,"total number of mallocs used during MatSetValues calls =%D\n",(PetscInt)info.mallocs);
752: }
753: }
754: }
755: if (mat->ops->view) {
756: PetscViewerASCIIPushTab(viewer);
757: (*mat->ops->view)(mat,viewer);
758: PetscViewerASCIIPopTab(viewer);
759: } else if (!iascii) {
760: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
761: }
762: if (iascii) {
763: PetscViewerGetFormat(viewer,&format);
764: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
765: PetscViewerASCIIPopTab(viewer);
766: }
767: }
768: PetscLogEventEnd(MAT_View,mat,viewer,0,0);
769: return(0);
770: }
772: #if defined(PETSC_USE_DEBUG)
773: #include <../src/sys/totalview/tv_data_display.h>
774: PETSC_UNUSED static int TV_display_type(const struct _p_Mat *mat)
775: {
776: TV_add_row("Local rows", "int", &mat->rmap->n);
777: TV_add_row("Local columns", "int", &mat->cmap->n);
778: TV_add_row("Global rows", "int", &mat->rmap->N);
779: TV_add_row("Global columns", "int", &mat->cmap->N);
780: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)mat)->type_name);
781: return TV_format_OK;
782: }
783: #endif
787: /*@C
788: MatLoad - Loads a matrix that has been stored in binary format
789: with MatView(). The matrix format is determined from the options database.
790: Generates a parallel MPI matrix if the communicator has more than one
791: processor. The default matrix type is AIJ.
793: Collective on PetscViewer
795: Input Parameters:
796: + newmat - the newly loaded matrix, this needs to have been created with MatCreate()
797: or some related function before a call to MatLoad()
798: - viewer - binary file viewer, created with PetscViewerBinaryOpen()
800: Options Database Keys:
801: Used with block matrix formats (MATSEQBAIJ, ...) to specify
802: block size
803: . -matload_block_size <bs>
805: Level: beginner
807: Notes:
808: If the Mat type has not yet been given then MATAIJ is used, call MatSetFromOptions() on the
809: Mat before calling this routine if you wish to set it from the options database.
811: MatLoad() automatically loads into the options database any options
812: given in the file filename.info where filename is the name of the file
813: that was passed to the PetscViewerBinaryOpen(). The options in the info
814: file will be ignored if you use the -viewer_binary_skip_info option.
816: If the type or size of newmat is not set before a call to MatLoad, PETSc
817: sets the default matrix type AIJ and sets the local and global sizes.
818: If type and/or size is already set, then the same are used.
820: In parallel, each processor can load a subset of rows (or the
821: entire matrix). This routine is especially useful when a large
822: matrix is stored on disk and only part of it is desired on each
823: processor. For example, a parallel solver may access only some of
824: the rows from each processor. The algorithm used here reads
825: relatively small blocks of data rather than reading the entire
826: matrix and then subsetting it.
828: Notes for advanced users:
829: Most users should not need to know the details of the binary storage
830: format, since MatLoad() and MatView() completely hide these details.
831: But for anyone who's interested, the standard binary matrix storage
832: format is
834: $ int MAT_FILE_CLASSID
835: $ int number of rows
836: $ int number of columns
837: $ int total number of nonzeros
838: $ int *number nonzeros in each row
839: $ int *column indices of all nonzeros (starting index is zero)
840: $ PetscScalar *values of all nonzeros
842: PETSc automatically does the byte swapping for
843: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
844: linux, Windows and the paragon; thus if you write your own binary
845: read/write routines you have to swap the bytes; see PetscBinaryRead()
846: and PetscBinaryWrite() to see how this may be done.
848: .keywords: matrix, load, binary, input
850: .seealso: PetscViewerBinaryOpen(), MatView(), VecLoad()
852: @*/
853: PetscErrorCode MatLoad(Mat newmat,PetscViewer viewer)
854: {
856: PetscBool isbinary,flg;
861: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
862: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
864: if (!((PetscObject)newmat)->type_name) {
865: MatSetType(newmat,MATAIJ);
866: }
868: if (!newmat->ops->load) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatLoad is not supported for type");
869: PetscLogEventBegin(MAT_Load,viewer,0,0,0);
870: (*newmat->ops->load)(newmat,viewer);
871: PetscLogEventEnd(MAT_Load,viewer,0,0,0);
873: flg = PETSC_FALSE;
874: PetscOptionsGetBool(((PetscObject)newmat)->prefix,"-matload_symmetric",&flg,PETSC_NULL);
875: if (flg) {
876: MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);
877: MatSetOption(newmat,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
878: }
879: flg = PETSC_FALSE;
880: PetscOptionsGetBool(((PetscObject)newmat)->prefix,"-matload_spd",&flg,PETSC_NULL);
881: if (flg) {
882: MatSetOption(newmat,MAT_SPD,PETSC_TRUE);
883: }
884: return(0);
885: }
889: /*@
890: MatScaleSystem - Scale a vector solution and right hand side to
891: match the scaling of a scaled matrix.
892:
893: Collective on Mat
895: Input Parameter:
896: + mat - the matrix
897: . b - right hand side vector (or PETSC_NULL)
898: - x - solution vector (or PETSC_NULL)
901: Notes:
902: For AIJ, and BAIJ matrix formats, the matrices are not
903: internally scaled, so this does nothing.
905: The KSP methods automatically call this routine when required
906: (via PCPreSolve()) so it is rarely used directly.
908: Level: Developer
910: Concepts: matrices^scaling
912: .seealso: MatUseScaledForm(), MatUnScaleSystem()
913: @*/
914: PetscErrorCode MatScaleSystem(Mat mat,Vec b,Vec x)
915: {
921: MatPreallocated(mat);
925: if (mat->ops->scalesystem) {
926: (*mat->ops->scalesystem)(mat,b,x);
927: }
928: return(0);
929: }
933: /*@
934: MatUnScaleSystem - Unscales a vector solution and right hand side to
935: match the original scaling of a scaled matrix.
936:
937: Collective on Mat
939: Input Parameter:
940: + mat - the matrix
941: . b - right hand side vector (or PETSC_NULL)
942: - x - solution vector (or PETSC_NULL)
945: Notes:
946: For AIJ and BAIJ matrix formats, the matrices are not
947: internally scaled, so this does nothing.
949: The KSP methods automatically call this routine when required
950: (via PCPreSolve()) so it is rarely used directly.
952: Level: Developer
954: .seealso: MatUseScaledForm(), MatScaleSystem()
955: @*/
956: PetscErrorCode MatUnScaleSystem(Mat mat,Vec b,Vec x)
957: {
963: MatPreallocated(mat);
966: if (mat->ops->unscalesystem) {
967: (*mat->ops->unscalesystem)(mat,b,x);
968: }
969: return(0);
970: }
974: /*@
975: MatUseScaledForm - For matrix storage formats that scale the
976: matrix indicates matrix operations (MatMult() etc) are
977: applied using the scaled matrix.
978:
979: Logically Collective on Mat
981: Input Parameter:
982: + mat - the matrix
983: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
984: applying the original matrix
986: Notes:
987: For scaled matrix formats, applying the original, unscaled matrix
988: will be slightly more expensive
990: Level: Developer
992: .seealso: MatScaleSystem(), MatUnScaleSystem()
993: @*/
994: PetscErrorCode MatUseScaledForm(Mat mat,PetscBool scaled)
995: {
1002: MatPreallocated(mat);
1003: if (mat->ops->usescaledform) {
1004: (*mat->ops->usescaledform)(mat,scaled);
1005: }
1006: return(0);
1007: }
1011: /*@
1012: MatDestroy - Frees space taken by a matrix.
1014: Collective on Mat
1016: Input Parameter:
1017: . A - the matrix
1019: Level: beginner
1021: @*/
1022: PetscErrorCode MatDestroy(Mat *A)
1023: {
1027: if (!*A) return(0);
1029: if (--((PetscObject)(*A))->refct > 0) {*A = PETSC_NULL; return(0);}
1031: /* if memory was published with AMS then destroy it */
1032: PetscObjectDepublish(*A);
1034: if ((*A)->ops->destroy) {
1035: (*(*A)->ops->destroy)(*A);
1036: }
1038: MatNullSpaceDestroy(&(*A)->nullsp);
1039: PetscLayoutDestroy(&(*A)->rmap);
1040: PetscLayoutDestroy(&(*A)->cmap);
1041: PetscHeaderDestroy(A);
1042: return(0);
1043: }
1047: /*@
1048: MatSetValues - Inserts or adds a block of values into a matrix.
1049: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1050: MUST be called after all calls to MatSetValues() have been completed.
1052: Not Collective
1054: Input Parameters:
1055: + mat - the matrix
1056: . v - a logically two-dimensional array of values
1057: . m, idxm - the number of rows and their global indices
1058: . n, idxn - the number of columns and their global indices
1059: - addv - either ADD_VALUES or INSERT_VALUES, where
1060: ADD_VALUES adds values to any existing entries, and
1061: INSERT_VALUES replaces existing entries with new values
1063: Notes:
1064: By default the values, v, are row-oriented. See MatSetOption() for other options.
1066: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
1067: options cannot be mixed without intervening calls to the assembly
1068: routines.
1070: MatSetValues() uses 0-based row and column numbers in Fortran
1071: as well as in C.
1073: Negative indices may be passed in idxm and idxn, these rows and columns are
1074: simply ignored. This allows easily inserting element stiffness matrices
1075: with homogeneous Dirchlet boundary conditions that you don't want represented
1076: in the matrix.
1078: Efficiency Alert:
1079: The routine MatSetValuesBlocked() may offer much better efficiency
1080: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1082: Level: beginner
1084: Concepts: matrices^putting entries in
1086: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1087: InsertMode, INSERT_VALUES, ADD_VALUES
1088: @*/
1089: PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1090: {
1096: if (!m || !n) return(0); /* no values to insert */
1100: MatPreallocated(mat);
1101: if (mat->insertmode == NOT_SET_VALUES) {
1102: mat->insertmode = addv;
1103: }
1104: #if defined(PETSC_USE_DEBUG)
1105: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1106: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1107: if (!mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1108: #endif
1110: if (mat->assembled) {
1111: mat->was_assembled = PETSC_TRUE;
1112: mat->assembled = PETSC_FALSE;
1113: }
1114: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1115: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
1116: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1117: #if defined(PETSC_HAVE_CUSP)
1118: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1119: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1120: }
1121: #endif
1122: return(0);
1123: }
1128: /*@
1129: MatSetValuesRowLocal - Inserts a row (block row for BAIJ matrices) of nonzero
1130: values into a matrix
1132: Not Collective
1134: Input Parameters:
1135: + mat - the matrix
1136: . row - the (block) row to set
1137: - v - a logically two-dimensional array of values
1139: Notes:
1140: By the values, v, are column-oriented (for the block version) and sorted
1142: All the nonzeros in the row must be provided
1144: The matrix must have previously had its column indices set
1146: The row must belong to this process
1148: Level: intermediate
1150: Concepts: matrices^putting entries in
1152: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1153: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues(), MatSetValuesRow(), MatSetLocalToGlobalMapping()
1154: @*/
1155: PetscErrorCode MatSetValuesRowLocal(Mat mat,PetscInt row,const PetscScalar v[])
1156: {
1163: MatSetValuesRow(mat, mat->rmap->mapping->indices[row],v);
1164: #if defined(PETSC_HAVE_CUSP)
1165: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1166: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1167: }
1168: #endif
1169: return(0);
1170: }
1174: /*@
1175: MatSetValuesRow - Inserts a row (block row for BAIJ matrices) of nonzero
1176: values into a matrix
1178: Not Collective
1180: Input Parameters:
1181: + mat - the matrix
1182: . row - the (block) row to set
1183: - v - a logically two-dimensional array of values
1185: Notes:
1186: The values, v, are column-oriented for the block version.
1188: All the nonzeros in the row must be provided
1190: THE MATRIX MUSAT HAVE PREVIOUSLY HAD ITS COLUMN INDICES SET. IT IS RARE THAT THIS ROUTINE IS USED, usually MatSetValues() is used.
1192: The row must belong to this process
1194: Level: advanced
1196: Concepts: matrices^putting entries in
1198: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1199: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
1200: @*/
1201: PetscErrorCode MatSetValuesRow(Mat mat,PetscInt row,const PetscScalar v[])
1202: {
1209: #if defined(PETSC_USE_DEBUG)
1210: if (mat->insertmode == ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add and insert values");
1211: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1212: #endif
1213: mat->insertmode = INSERT_VALUES;
1215: if (mat->assembled) {
1216: mat->was_assembled = PETSC_TRUE;
1217: mat->assembled = PETSC_FALSE;
1218: }
1219: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1220: if (!mat->ops->setvaluesrow) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1221: (*mat->ops->setvaluesrow)(mat,row,v);
1222: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1223: #if defined(PETSC_HAVE_CUSP)
1224: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1225: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1226: }
1227: #endif
1228: return(0);
1229: }
1233: /*@
1234: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
1235: Using structured grid indexing
1237: Not Collective
1239: Input Parameters:
1240: + mat - the matrix
1241: . m - number of rows being entered
1242: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
1243: . n - number of columns being entered
1244: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
1245: . v - a logically two-dimensional array of values
1246: - addv - either ADD_VALUES or INSERT_VALUES, where
1247: ADD_VALUES adds values to any existing entries, and
1248: INSERT_VALUES replaces existing entries with new values
1250: Notes:
1251: By default the values, v, are row-oriented. See MatSetOption() for other options.
1253: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
1254: options cannot be mixed without intervening calls to the assembly
1255: routines.
1257: The grid coordinates are across the entire grid, not just the local portion
1259: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
1260: as well as in C.
1262: For setting/accessing vector values via array coordinates you can use the DMDAVecGetArray() routine
1264: In order to use this routine you must either obtain the matrix with DMGetMatrix()
1265: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1267: The columns and rows in the stencil passed in MUST be contained within the
1268: ghost region of the given process as set with DMDACreateXXX() or MatSetStencil(). For example,
1269: if you create a DMDA with an overlap of one grid level and on a particular process its first
1270: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1271: first i index you can use in your column and row indices in MatSetStencil() is 5.
1273: In Fortran idxm and idxn should be declared as
1274: $ MatStencil idxm(4,m),idxn(4,n)
1275: and the values inserted using
1276: $ idxm(MatStencil_i,1) = i
1277: $ idxm(MatStencil_j,1) = j
1278: $ idxm(MatStencil_k,1) = k
1279: $ idxm(MatStencil_c,1) = c
1280: etc
1281:
1282: For periodic boundary conditions use negative indices for values to the left (below 0; that are to be
1283: obtained by wrapping values from right edge). For values to the right of the last entry using that index plus one
1284: etc to obtain values that obtained by wrapping the values from the left edge. This does not work for anything but the
1285: DMDA_BOUNDARY_PERIODIC boundary type.
1287: 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
1288: a single value per point) you can skip filling those indices.
1290: Inspired by the structured grid interface to the HYPRE package
1291: (http://www.llnl.gov/CASC/hypre)
1293: Efficiency Alert:
1294: The routine MatSetValuesBlockedStencil() may offer much better efficiency
1295: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1297: Level: beginner
1299: Concepts: matrices^putting entries in
1301: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1302: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DMGetMatrix(), DMDAVecGetArray(), MatStencil
1303: @*/
1304: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1305: {
1307: PetscInt buf[8192],*bufm=0,*bufn=0,*jdxm,*jdxn;
1308: PetscInt j,i,dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1309: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1312: if (!m || !n) return(0); /* no values to insert */
1319: if ((m+n) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1320: jdxm = buf; jdxn = buf+m;
1321: } else {
1322: PetscMalloc2(m,PetscInt,&bufm,n,PetscInt,&bufn);
1323: jdxm = bufm; jdxn = bufn;
1324: }
1325: for (i=0; i<m; i++) {
1326: for (j=0; j<3-sdim; j++) dxm++;
1327: tmp = *dxm++ - starts[0];
1328: for (j=0; j<dim-1; j++) {
1329: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1330: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1331: }
1332: if (mat->stencil.noc) dxm++;
1333: jdxm[i] = tmp;
1334: }
1335: for (i=0; i<n; i++) {
1336: for (j=0; j<3-sdim; j++) dxn++;
1337: tmp = *dxn++ - starts[0];
1338: for (j=0; j<dim-1; j++) {
1339: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1340: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1341: }
1342: if (mat->stencil.noc) dxn++;
1343: jdxn[i] = tmp;
1344: }
1345: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1346: PetscFree2(bufm,bufn);
1347: return(0);
1348: }
1352: /*@C
1353: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1354: Using structured grid indexing
1356: Not Collective
1358: Input Parameters:
1359: + mat - the matrix
1360: . m - number of rows being entered
1361: . idxm - grid coordinates for matrix rows being entered
1362: . n - number of columns being entered
1363: . idxn - grid coordinates for matrix columns being entered
1364: . v - a logically two-dimensional array of values
1365: - addv - either ADD_VALUES or INSERT_VALUES, where
1366: ADD_VALUES adds values to any existing entries, and
1367: INSERT_VALUES replaces existing entries with new values
1369: Notes:
1370: By default the values, v, are row-oriented and unsorted.
1371: See MatSetOption() for other options.
1373: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1374: options cannot be mixed without intervening calls to the assembly
1375: routines.
1377: The grid coordinates are across the entire grid, not just the local portion
1379: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1380: as well as in C.
1382: For setting/accessing vector values via array coordinates you can use the DMDAVecGetArray() routine
1384: In order to use this routine you must either obtain the matrix with DMGetMatrix()
1385: or call MatSetBlockSize(), MatSetLocalToGlobalMapping() and MatSetStencil() first.
1387: The columns and rows in the stencil passed in MUST be contained within the
1388: ghost region of the given process as set with DMDACreateXXX() or MatSetStencil(). For example,
1389: if you create a DMDA with an overlap of one grid level and on a particular process its first
1390: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1391: first i index you can use in your column and row indices in MatSetStencil() is 5.
1393: In Fortran idxm and idxn should be declared as
1394: $ MatStencil idxm(4,m),idxn(4,n)
1395: and the values inserted using
1396: $ idxm(MatStencil_i,1) = i
1397: $ idxm(MatStencil_j,1) = j
1398: $ idxm(MatStencil_k,1) = k
1399: etc
1401: Negative indices may be passed in idxm and idxn, these rows and columns are
1402: simply ignored. This allows easily inserting element stiffness matrices
1403: with homogeneous Dirchlet boundary conditions that you don't want represented
1404: in the matrix.
1406: Inspired by the structured grid interface to the HYPRE package
1407: (http://www.llnl.gov/CASC/hypre)
1409: Level: beginner
1411: Concepts: matrices^putting entries in
1413: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1414: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DMGetMatrix(), DMDAVecGetArray(), MatStencil,
1415: MatSetBlockSize(), MatSetLocalToGlobalMapping()
1416: @*/
1417: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1418: {
1420: PetscInt buf[8192],*bufm=0,*bufn=0,*jdxm,*jdxn;
1421: PetscInt j,i,dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1422: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1425: if (!m || !n) return(0); /* no values to insert */
1432: if ((m+n) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1433: jdxm = buf; jdxn = buf+m;
1434: } else {
1435: PetscMalloc2(m,PetscInt,&bufm,n,PetscInt,&bufn);
1436: jdxm = bufm; jdxn = bufn;
1437: }
1438: for (i=0; i<m; i++) {
1439: for (j=0; j<3-sdim; j++) dxm++;
1440: tmp = *dxm++ - starts[0];
1441: for (j=0; j<sdim-1; j++) {
1442: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1443: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1444: }
1445: dxm++;
1446: jdxm[i] = tmp;
1447: }
1448: for (i=0; i<n; i++) {
1449: for (j=0; j<3-sdim; j++) dxn++;
1450: tmp = *dxn++ - starts[0];
1451: for (j=0; j<sdim-1; j++) {
1452: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1453: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1454: }
1455: dxn++;
1456: jdxn[i] = tmp;
1457: }
1458: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1459: PetscFree2(bufm,bufn);
1460: #if defined(PETSC_HAVE_CUSP)
1461: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1462: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1463: }
1464: #endif
1465: return(0);
1466: }
1470: /*@
1471: MatSetStencil - Sets the grid information for setting values into a matrix via
1472: MatSetValuesStencil()
1474: Not Collective
1476: Input Parameters:
1477: + mat - the matrix
1478: . dim - dimension of the grid 1, 2, or 3
1479: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1480: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1481: - dof - number of degrees of freedom per node
1484: Inspired by the structured grid interface to the HYPRE package
1485: (www.llnl.gov/CASC/hyper)
1487: For matrices generated with DMGetMatrix() this routine is automatically called and so not needed by the
1488: user.
1489:
1490: Level: beginner
1492: Concepts: matrices^putting entries in
1494: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1495: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1496: @*/
1497: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1498: {
1499: PetscInt i;
1506: mat->stencil.dim = dim + (dof > 1);
1507: for (i=0; i<dim; i++) {
1508: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1509: mat->stencil.starts[i] = starts[dim-i-1];
1510: }
1511: mat->stencil.dims[dim] = dof;
1512: mat->stencil.starts[dim] = 0;
1513: mat->stencil.noc = (PetscBool)(dof == 1);
1514: return(0);
1515: }
1519: /*@
1520: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1522: Not Collective
1524: Input Parameters:
1525: + mat - the matrix
1526: . v - a logically two-dimensional array of values
1527: . m, idxm - the number of block rows and their global block indices
1528: . n, idxn - the number of block columns and their global block indices
1529: - addv - either ADD_VALUES or INSERT_VALUES, where
1530: ADD_VALUES adds values to any existing entries, and
1531: INSERT_VALUES replaces existing entries with new values
1533: Notes:
1534: The m and n count the NUMBER of blocks in the row direction and column direction,
1535: NOT the total number of rows/columns; for example, if the block size is 2 and
1536: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1537: The values in idxm would be 1 2; that is the first index for each block divided by
1538: the block size.
1540: Note that you must call MatSetBlockSize() when constructing this matrix (after
1541: preallocating it).
1543: By default the values, v, are row-oriented, so the layout of
1544: v is the same as for MatSetValues(). See MatSetOption() for other options.
1546: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1547: options cannot be mixed without intervening calls to the assembly
1548: routines.
1550: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1551: as well as in C.
1553: Negative indices may be passed in idxm and idxn, these rows and columns are
1554: simply ignored. This allows easily inserting element stiffness matrices
1555: with homogeneous Dirchlet boundary conditions that you don't want represented
1556: in the matrix.
1558: Each time an entry is set within a sparse matrix via MatSetValues(),
1559: internal searching must be done to determine where to place the the
1560: data in the matrix storage space. By instead inserting blocks of
1561: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1562: reduced.
1564: Example:
1565: $ Suppose m=n=2 and block size(bs) = 2 The array is
1566: $
1567: $ 1 2 | 3 4
1568: $ 5 6 | 7 8
1569: $ - - - | - - -
1570: $ 9 10 | 11 12
1571: $ 13 14 | 15 16
1572: $
1573: $ v[] should be passed in like
1574: $ v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
1575: $
1576: $ If you are not using row oriented storage of v (that is you called MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE)) then
1577: $ v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
1579: Level: intermediate
1581: Concepts: matrices^putting entries in blocked
1583: .seealso: MatSetBlockSize(), MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1584: @*/
1585: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1586: {
1592: if (!m || !n) return(0); /* no values to insert */
1596: MatPreallocated(mat);
1597: if (mat->insertmode == NOT_SET_VALUES) {
1598: mat->insertmode = addv;
1599: }
1600: #if defined(PETSC_USE_DEBUG)
1601: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1602: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1603: if (!mat->ops->setvaluesblocked && !mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1604: #endif
1606: if (mat->assembled) {
1607: mat->was_assembled = PETSC_TRUE;
1608: mat->assembled = PETSC_FALSE;
1609: }
1610: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1611: if (mat->ops->setvaluesblocked) {
1612: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1613: } else {
1614: PetscInt buf[8192],*bufr=0,*bufc=0,*iidxm,*iidxn;
1615: PetscInt i,j,bs=mat->rmap->bs;
1616: if ((m+n)*bs <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1617: iidxm = buf; iidxn = buf + m*bs;
1618: } else {
1619: PetscMalloc2(m*bs,PetscInt,&bufr,n*bs,PetscInt,&bufc);
1620: iidxm = bufr; iidxn = bufc;
1621: }
1622: for (i=0; i<m; i++)
1623: for (j=0; j<bs; j++)
1624: iidxm[i*bs+j] = bs*idxm[i] + j;
1625: for (i=0; i<n; i++)
1626: for (j=0; j<bs; j++)
1627: iidxn[i*bs+j] = bs*idxn[i] + j;
1628: MatSetValues(mat,m*bs,iidxm,n*bs,iidxn,v,addv);
1629: PetscFree2(bufr,bufc);
1630: }
1631: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1632: #if defined(PETSC_HAVE_CUSP)
1633: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1634: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1635: }
1636: #endif
1637: return(0);
1638: }
1642: /*@
1643: MatGetValues - Gets a block of values from a matrix.
1645: Not Collective; currently only returns a local block
1647: Input Parameters:
1648: + mat - the matrix
1649: . v - a logically two-dimensional array for storing the values
1650: . m, idxm - the number of rows and their global indices
1651: - n, idxn - the number of columns and their global indices
1653: Notes:
1654: The user must allocate space (m*n PetscScalars) for the values, v.
1655: The values, v, are then returned in a row-oriented format,
1656: analogous to that used by default in MatSetValues().
1658: MatGetValues() uses 0-based row and column numbers in
1659: Fortran as well as in C.
1661: MatGetValues() requires that the matrix has been assembled
1662: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1663: MatSetValues() and MatGetValues() CANNOT be made in succession
1664: without intermediate matrix assembly.
1666: Negative row or column indices will be ignored and those locations in v[] will be
1667: left unchanged.
1669: Level: advanced
1671: Concepts: matrices^accessing values
1673: .seealso: MatGetRow(), MatGetSubMatrices(), MatSetValues()
1674: @*/
1675: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1676: {
1682: if (!m || !n) return(0);
1686: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1687: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1688: if (!mat->ops->getvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1689: MatPreallocated(mat);
1691: PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
1692: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1693: PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
1694: return(0);
1695: }
1699: /*@
1700: MatSetValuesBatch - Adds (ADD_VALUES) many blocks of values into a matrix at once. The blocks must all be square and
1701: the same size. Currently, this can only be called once and creates the given matrix.
1703: Not Collective
1705: Input Parameters:
1706: + mat - the matrix
1707: . nb - the number of blocks
1708: . bs - the number of rows (and columns) in each block
1709: . rows - a concatenation of the rows for each block
1710: - v - a concatenation of logically two-dimensional arrays of values
1712: Notes:
1713: In the future, we will extend this routine to handle rectangular blocks, and to allow multiple calls for a given matrix.
1715: Level: advanced
1717: Concepts: matrices^putting entries in
1719: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1720: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
1721: @*/
1722: PetscErrorCode MatSetValuesBatch(Mat mat, PetscInt nb, PetscInt bs, PetscInt rows[], const PetscScalar v[])
1723: {
1731: #if defined(PETSC_USE_DEBUG)
1732: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1733: #endif
1735: if (mat->ops->setvaluesbatch) {
1736: PetscLogEventBegin(MAT_SetValuesBatch,mat,0,0,0);
1737: (*mat->ops->setvaluesbatch)(mat,nb,bs,rows,v);
1738: PetscLogEventEnd(MAT_SetValuesBatch,mat,0,0,0);
1739: } else {
1740: PetscInt b;
1741: for(b = 0; b < nb; ++b) {
1742: MatSetValues(mat, bs, &rows[b*bs], bs, &rows[b*bs], &v[b*bs*bs], ADD_VALUES);
1743: }
1744: }
1745: return(0);
1746: }
1750: /*@
1751: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1752: the routine MatSetValuesLocal() to allow users to insert matrix entries
1753: using a local (per-processor) numbering.
1755: Not Collective
1757: Input Parameters:
1758: + x - the matrix
1759: . rmapping - row mapping created with ISLocalToGlobalMappingCreate()
1760: or ISLocalToGlobalMappingCreateIS()
1761: - cmapping - column mapping
1763: Level: intermediate
1765: Concepts: matrices^local to global mapping
1766: Concepts: local to global mapping^for matrices
1768: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1769: @*/
1770: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1771: {
1778: MatPreallocated(x);
1780: if (x->ops->setlocaltoglobalmapping) {
1781: (*x->ops->setlocaltoglobalmapping)(x,rmapping,cmapping);
1782: } else {
1783: PetscLayoutSetISLocalToGlobalMapping(x->rmap,rmapping);
1784: PetscLayoutSetISLocalToGlobalMapping(x->cmap,cmapping);
1785: }
1786: return(0);
1787: }
1791: /*@
1792: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
1793: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
1794: entries using a local (per-processor) numbering.
1796: Not Collective
1798: Input Parameters:
1799: + x - the matrix
1800: . rmapping - row mapping created with ISLocalToGlobalMappingCreate() or
1801: ISLocalToGlobalMappingCreateIS()
1802: - cmapping - column mapping
1804: Level: intermediate
1806: Concepts: matrices^local to global mapping blocked
1807: Concepts: local to global mapping^for matrices, blocked
1809: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
1810: MatSetValuesBlocked(), MatSetValuesLocal()
1811: @*/
1812: PetscErrorCode MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1813: {
1820: MatPreallocated(x);
1822: PetscLayoutSetISLocalToGlobalMappingBlock(x->rmap,rmapping);
1823: PetscLayoutSetISLocalToGlobalMappingBlock(x->cmap,cmapping);
1824: return(0);
1825: }
1829: /*@
1830: MatGetLocalToGlobalMapping - Gets the local-to-global numbering set by MatSetLocalToGlobalMapping()
1832: Not Collective
1834: Input Parameters:
1835: . A - the matrix
1837: Output Parameters:
1838: + rmapping - row mapping
1839: - cmapping - column mapping
1841: Level: advanced
1843: Concepts: matrices^local to global mapping
1844: Concepts: local to global mapping^for matrices
1846: .seealso: MatSetValuesLocal(), MatGetLocalToGlobalMappingBlock()
1847: @*/
1848: PetscErrorCode MatGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
1849: {
1855: if (rmapping) *rmapping = A->rmap->mapping;
1856: if (cmapping) *cmapping = A->cmap->mapping;
1857: return(0);
1858: }
1862: /*@
1863: MatGetLocalToGlobalMappingBlock - Gets the local-to-global numbering set by MatSetLocalToGlobalMappingBlock()
1865: Not Collective
1867: Input Parameters:
1868: . A - the matrix
1870: Output Parameters:
1871: + rmapping - row mapping
1872: - cmapping - column mapping
1874: Level: advanced
1876: Concepts: matrices^local to global mapping blocked
1877: Concepts: local to global mapping^for matrices, blocked
1879: .seealso: MatSetValuesBlockedLocal(), MatGetLocalToGlobalMapping()
1880: @*/
1881: PetscErrorCode MatGetLocalToGlobalMappingBlock(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
1882: {
1888: if (rmapping) *rmapping = A->rmap->bmapping;
1889: if (cmapping) *cmapping = A->cmap->bmapping;
1890: return(0);
1891: }
1895: /*@
1896: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
1897: using a local ordering of the nodes.
1899: Not Collective
1901: Input Parameters:
1902: + x - the matrix
1903: . nrow, irow - number of rows and their local indices
1904: . ncol, icol - number of columns and their local indices
1905: . y - a logically two-dimensional array of values
1906: - addv - either INSERT_VALUES or ADD_VALUES, where
1907: ADD_VALUES adds values to any existing entries, and
1908: INSERT_VALUES replaces existing entries with new values
1910: Notes:
1911: Before calling MatSetValuesLocal(), the user must first set the
1912: local-to-global mapping by calling MatSetLocalToGlobalMapping().
1914: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
1915: options cannot be mixed without intervening calls to the assembly
1916: routines.
1918: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1919: MUST be called after all calls to MatSetValuesLocal() have been completed.
1921: Level: intermediate
1923: Concepts: matrices^putting entries in with local numbering
1925: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
1926: MatSetValueLocal()
1927: @*/
1928: PetscErrorCode MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1929: {
1935: if (!nrow || !ncol) return(0); /* no values to insert */
1939: MatPreallocated(mat);
1940: if (mat->insertmode == NOT_SET_VALUES) {
1941: mat->insertmode = addv;
1942: }
1943: #if defined(PETSC_USE_DEBUG)
1944: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1945: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1946: if (!mat->ops->setvalueslocal && !mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1947: #endif
1949: if (mat->assembled) {
1950: mat->was_assembled = PETSC_TRUE;
1951: mat->assembled = PETSC_FALSE;
1952: }
1953: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1954: if (mat->ops->setvalueslocal) {
1955: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1956: } else {
1957: PetscInt buf[8192],*bufr=0,*bufc=0,*irowm,*icolm;
1958: if ((nrow+ncol) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1959: irowm = buf; icolm = buf+nrow;
1960: } else {
1961: PetscMalloc2(nrow,PetscInt,&bufr,ncol,PetscInt,&bufc);
1962: irowm = bufr; icolm = bufc;
1963: }
1964: ISLocalToGlobalMappingApply(mat->rmap->mapping,nrow,irow,irowm);
1965: ISLocalToGlobalMappingApply(mat->cmap->mapping,ncol,icol,icolm);
1966: MatSetValues(mat,nrow,irowm,ncol,icolm,y,addv);
1967: PetscFree2(bufr,bufc);
1968: }
1969: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1970: #if defined(PETSC_HAVE_CUSP)
1971: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1972: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1973: }
1974: #endif
1975: return(0);
1976: }
1980: /*@
1981: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1982: using a local ordering of the nodes a block at a time.
1984: Not Collective
1986: Input Parameters:
1987: + x - the matrix
1988: . nrow, irow - number of rows and their local indices
1989: . ncol, icol - number of columns and their local indices
1990: . y - a logically two-dimensional array of values
1991: - addv - either INSERT_VALUES or ADD_VALUES, where
1992: ADD_VALUES adds values to any existing entries, and
1993: INSERT_VALUES replaces existing entries with new values
1995: Notes:
1996: Before calling MatSetValuesBlockedLocal(), the user must first set the
1997: block size using MatSetBlockSize(), and the local-to-global mapping by
1998: calling MatSetLocalToGlobalMappingBlock(), where the mapping MUST be
1999: set for matrix blocks, not for matrix elements.
2001: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
2002: options cannot be mixed without intervening calls to the assembly
2003: routines.
2005: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
2006: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
2008: Level: intermediate
2010: Concepts: matrices^putting blocked values in with local numbering
2012: .seealso: MatSetBlockSize(), MatSetLocalToGlobalMappingBlock(), MatAssemblyBegin(), MatAssemblyEnd(),
2013: MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
2014: @*/
2015: PetscErrorCode MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
2016: {
2022: if (!nrow || !ncol) return(0); /* no values to insert */
2026: MatPreallocated(mat);
2027: if (mat->insertmode == NOT_SET_VALUES) {
2028: mat->insertmode = addv;
2029: }
2030: #if defined(PETSC_USE_DEBUG)
2031: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
2032: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2033: if (!mat->ops->setvaluesblockedlocal && !mat->ops->setvaluesblocked && !mat->ops->setvalueslocal && !mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2034: #endif
2036: if (mat->assembled) {
2037: mat->was_assembled = PETSC_TRUE;
2038: mat->assembled = PETSC_FALSE;
2039: }
2040: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
2041: if (mat->ops->setvaluesblockedlocal) {
2042: (*mat->ops->setvaluesblockedlocal)(mat,nrow,irow,ncol,icol,y,addv);
2043: } else {
2044: PetscInt buf[8192],*bufr=0,*bufc=0,*irowm,*icolm;
2045: if (mat->rmap->bmapping && mat->cmap->bmapping) {
2046: if ((nrow+ncol) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
2047: irowm = buf; icolm = buf + nrow;
2048: } else {
2049: PetscMalloc2(nrow,PetscInt,&bufr,ncol,PetscInt,&bufc);
2050: irowm = bufr; icolm = bufc;
2051: }
2052: ISLocalToGlobalMappingApply(mat->rmap->bmapping,nrow,irow,irowm);
2053: ISLocalToGlobalMappingApply(mat->cmap->bmapping,ncol,icol,icolm);
2054: MatSetValuesBlocked(mat,nrow,irowm,ncol,icolm,y,addv);
2055: PetscFree2(bufr,bufc);
2056: } else {
2057: PetscInt i,j,bs=mat->rmap->bs;
2058: if ((nrow+ncol)*bs <=(PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
2059: irowm = buf; icolm = buf + nrow;
2060: } else {
2061: PetscMalloc2(nrow*bs,PetscInt,&bufr,ncol*bs,PetscInt,&bufc);
2062: irowm = bufr; icolm = bufc;
2063: }
2064: for (i=0; i<nrow; i++)
2065: for (j=0; j<bs; j++)
2066: irowm[i*bs+j] = irow[i]*bs+j;
2067: for (i=0; i<ncol; i++)
2068: for (j=0; j<bs; j++)
2069: icolm[i*bs+j] = icol[i]*bs+j;
2070: MatSetValuesLocal(mat,nrow*bs,irowm,ncol*bs,icolm,y,addv);
2071: PetscFree2(bufr,bufc);
2072: }
2073: }
2074: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
2075: #if defined(PETSC_HAVE_CUSP)
2076: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
2077: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
2078: }
2079: #endif
2080: return(0);
2081: }
2085: /*@
2086: MatMultDiagonalBlock - Computes the matrix-vector product, y = Dx. Where D is defined by the inode or block structure of the diagonal
2088: Collective on Mat and Vec
2090: Input Parameters:
2091: + mat - the matrix
2092: - x - the vector to be multiplied
2094: Output Parameters:
2095: . y - the result
2097: Notes:
2098: The vectors x and y cannot be the same. I.e., one cannot
2099: call MatMult(A,y,y).
2101: Level: developer
2103: Concepts: matrix-vector product
2105: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2106: @*/
2107: PetscErrorCode MatMultDiagonalBlock(Mat mat,Vec x,Vec y)
2108: {
2117: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2118: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2119: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2120: MatPreallocated(mat);
2122: if (!mat->ops->multdiagonalblock) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
2123: (*mat->ops->multdiagonalblock)(mat,x,y);
2124: PetscObjectStateIncrease((PetscObject)y);
2125: return(0);
2126: }
2128: /* --------------------------------------------------------*/
2131: /*@
2132: MatMult - Computes the matrix-vector product, y = Ax.
2134: Neighbor-wise Collective on Mat and Vec
2136: Input Parameters:
2137: + mat - the matrix
2138: - x - the vector to be multiplied
2140: Output Parameters:
2141: . y - the result
2143: Notes:
2144: The vectors x and y cannot be the same. I.e., one cannot
2145: call MatMult(A,y,y).
2147: Level: beginner
2149: Concepts: matrix-vector product
2151: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2152: @*/
2153: PetscErrorCode MatMult(Mat mat,Vec x,Vec y)
2154: {
2162: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2163: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2164: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2165: #ifndef PETSC_HAVE_CONSTRAINTS
2166: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2167: if (mat->rmap->N != y->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
2168: if (mat->rmap->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->rmap->n,y->map->n);
2169: #endif
2170: MatPreallocated(mat);
2172: if (mat->nullsp) {
2173: MatNullSpaceRemove(mat->nullsp,x,&x);
2174: }
2176: if (!mat->ops->mult) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
2177: PetscLogEventBegin(MAT_Mult,mat,x,y,0);
2178: (*mat->ops->mult)(mat,x,y);
2179: PetscLogEventEnd(MAT_Mult,mat,x,y,0);
2181: if (mat->nullsp) {
2182: MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
2183: }
2184: PetscObjectStateIncrease((PetscObject)y);
2185: return(0);
2186: }
2190: /*@
2191: MatMultTranspose - Computes matrix transpose times a vector.
2193: Neighbor-wise Collective on Mat and Vec
2195: Input Parameters:
2196: + mat - the matrix
2197: - x - the vector to be multilplied
2199: Output Parameters:
2200: . y - the result
2202: Notes:
2203: The vectors x and y cannot be the same. I.e., one cannot
2204: call MatMultTranspose(A,y,y).
2206: For complex numbers this does NOT compute the Hermitian (complex conjugate) transpose multiple,
2207: use MatMultHermitianTranspose()
2209: Level: beginner
2211: Concepts: matrix vector product^transpose
2213: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd(), MatMultHermitianTranspose(), MatTranspose()
2214: @*/
2215: PetscErrorCode MatMultTranspose(Mat mat,Vec x,Vec y)
2216: {
2225: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2226: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2227: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2228: #ifndef PETSC_HAVE_CONSTRAINTS
2229: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
2230: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->cmap->N,y->map->N);
2231: #endif
2232: MatPreallocated(mat);
2234: if (!mat->ops->multtranspose) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
2235: PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
2236: (*mat->ops->multtranspose)(mat,x,y);
2237: PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
2238: PetscObjectStateIncrease((PetscObject)y);
2239: return(0);
2240: }
2244: /*@
2245: MatMultHermitianTranspose - Computes matrix Hermitian transpose times a vector.
2247: Neighbor-wise Collective on Mat and Vec
2249: Input Parameters:
2250: + mat - the matrix
2251: - x - the vector to be multilplied
2253: Output Parameters:
2254: . y - the result
2256: Notes:
2257: The vectors x and y cannot be the same. I.e., one cannot
2258: call MatMultHermitianTranspose(A,y,y).
2260: Also called the conjugate transpose, complex conjugate transpose, or adjoint.
2262: For real numbers MatMultTranspose() and MatMultHermitianTranspose() are identical.
2264: Level: beginner
2266: Concepts: matrix vector product^transpose
2268: .seealso: MatMult(), MatMultAdd(), MatMultHermitianTransposeAdd(), MatMultTranspose()
2269: @*/
2270: PetscErrorCode MatMultHermitianTranspose(Mat mat,Vec x,Vec y)
2271: {
2280: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2281: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2282: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2283: #ifndef PETSC_HAVE_CONSTRAINTS
2284: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
2285: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->cmap->N,y->map->N);
2286: #endif
2287: MatPreallocated(mat);
2289: if (!mat->ops->multhermitiantranspose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2290: PetscLogEventBegin(MAT_MultHermitianTranspose,mat,x,y,0);
2291: (*mat->ops->multhermitiantranspose)(mat,x,y);
2292: PetscLogEventEnd(MAT_MultHermitianTranspose,mat,x,y,0);
2293: PetscObjectStateIncrease((PetscObject)y);
2294: return(0);
2295: }
2299: /*@
2300: MatMultAdd - Computes v3 = v2 + A * v1.
2302: Neighbor-wise Collective on Mat and Vec
2304: Input Parameters:
2305: + mat - the matrix
2306: - v1, v2 - the vectors
2308: Output Parameters:
2309: . v3 - the result
2311: Notes:
2312: The vectors v1 and v3 cannot be the same. I.e., one cannot
2313: call MatMultAdd(A,v1,v2,v1).
2315: Level: beginner
2317: Concepts: matrix vector product^addition
2319: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
2320: @*/
2321: PetscErrorCode MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2322: {
2332: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2333: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2334: if (mat->cmap->N != v1->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->cmap->N,v1->map->N);
2335: /* if (mat->rmap->N != v2->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->rmap->N,v2->map->N);
2336: if (mat->rmap->N != v3->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->rmap->N,v3->map->N); */
2337: if (mat->rmap->n != v3->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %D %D",mat->rmap->n,v3->map->n);
2338: if (mat->rmap->n != v2->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %D %D",mat->rmap->n,v2->map->n);
2339: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2340: MatPreallocated(mat);
2342: if (!mat->ops->multadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"No MatMultAdd() for matrix type '%s'",((PetscObject)mat)->type_name);
2343: PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
2344: (*mat->ops->multadd)(mat,v1,v2,v3);
2345: PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
2346: PetscObjectStateIncrease((PetscObject)v3);
2347: return(0);
2348: }
2352: /*@
2353: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
2355: Neighbor-wise Collective on Mat and Vec
2357: Input Parameters:
2358: + mat - the matrix
2359: - v1, v2 - the vectors
2361: Output Parameters:
2362: . v3 - the result
2364: Notes:
2365: The vectors v1 and v3 cannot be the same. I.e., one cannot
2366: call MatMultTransposeAdd(A,v1,v2,v1).
2368: Level: beginner
2370: Concepts: matrix vector product^transpose and addition
2372: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
2373: @*/
2374: PetscErrorCode MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2375: {
2385: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2386: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2387: if (!mat->ops->multtransposeadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2388: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2389: if (mat->rmap->N != v1->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->rmap->N,v1->map->N);
2390: if (mat->cmap->N != v2->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->cmap->N,v2->map->N);
2391: if (mat->cmap->N != v3->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->cmap->N,v3->map->N);
2392: MatPreallocated(mat);
2394: PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
2395: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
2396: PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
2397: PetscObjectStateIncrease((PetscObject)v3);
2398: return(0);
2399: }
2403: /*@
2404: MatMultHermitianTransposeAdd - Computes v3 = v2 + A^H * v1.
2406: Neighbor-wise Collective on Mat and Vec
2408: Input Parameters:
2409: + mat - the matrix
2410: - v1, v2 - the vectors
2412: Output Parameters:
2413: . v3 - the result
2415: Notes:
2416: The vectors v1 and v3 cannot be the same. I.e., one cannot
2417: call MatMultHermitianTransposeAdd(A,v1,v2,v1).
2419: Level: beginner
2421: Concepts: matrix vector product^transpose and addition
2423: .seealso: MatMultHermitianTranspose(), MatMultTranspose(), MatMultAdd(), MatMult()
2424: @*/
2425: PetscErrorCode MatMultHermitianTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2426: {
2436: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2437: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2438: if (!mat->ops->multhermitiantransposeadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2439: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2440: if (mat->rmap->N != v1->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->rmap->N,v1->map->N);
2441: if (mat->cmap->N != v2->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->cmap->N,v2->map->N);
2442: if (mat->cmap->N != v3->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->cmap->N,v3->map->N);
2443: MatPreallocated(mat);
2445: PetscLogEventBegin(MAT_MultHermitianTransposeAdd,mat,v1,v2,v3);
2446: (*mat->ops->multhermitiantransposeadd)(mat,v1,v2,v3);
2447: PetscLogEventEnd(MAT_MultHermitianTransposeAdd,mat,v1,v2,v3);
2448: PetscObjectStateIncrease((PetscObject)v3);
2449: return(0);
2450: }
2454: /*@
2455: MatMultConstrained - The inner multiplication routine for a
2456: constrained matrix P^T A P.
2458: Neighbor-wise Collective on Mat and Vec
2460: Input Parameters:
2461: + mat - the matrix
2462: - x - the vector to be multilplied
2464: Output Parameters:
2465: . y - the result
2467: Notes:
2468: The vectors x and y cannot be the same. I.e., one cannot
2469: call MatMult(A,y,y).
2471: Level: beginner
2473: .keywords: matrix, multiply, matrix-vector product, constraint
2474: .seealso: MatMult(), MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2475: @*/
2476: PetscErrorCode MatMultConstrained(Mat mat,Vec x,Vec y)
2477: {
2484: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2485: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2486: if (x == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2487: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2488: if (mat->rmap->N != y->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
2489: if (mat->rmap->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->rmap->n,y->map->n);
2491: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2492: (*mat->ops->multconstrained)(mat,x,y);
2493: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2494: PetscObjectStateIncrease((PetscObject)y);
2496: return(0);
2497: }
2501: /*@
2502: MatMultTransposeConstrained - The inner multiplication routine for a
2503: constrained matrix P^T A^T P.
2505: Neighbor-wise Collective on Mat and Vec
2507: Input Parameters:
2508: + mat - the matrix
2509: - x - the vector to be multilplied
2511: Output Parameters:
2512: . y - the result
2514: Notes:
2515: The vectors x and y cannot be the same. I.e., one cannot
2516: call MatMult(A,y,y).
2518: Level: beginner
2520: .keywords: matrix, multiply, matrix-vector product, constraint
2521: .seealso: MatMult(), MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2522: @*/
2523: PetscErrorCode MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
2524: {
2531: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2532: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2533: if (x == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2534: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2535: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
2537: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2538: (*mat->ops->multtransposeconstrained)(mat,x,y);
2539: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2540: PetscObjectStateIncrease((PetscObject)y);
2542: return(0);
2543: }
2547: /*@C
2548: MatGetFactorType - gets the type of factorization it is
2550: Note Collective
2551: as the flag
2553: Input Parameters:
2554: . mat - the matrix
2556: Output Parameters:
2557: . t - the type, one of MAT_FACTOR_NONE, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ILU, MAT_FACTOR_ICC,MAT_FACTOR_ILUDT
2559: Level: intermediate
2561: .seealso: MatFactorType, MatGetFactor()
2562: @*/
2563: PetscErrorCode MatGetFactorType(Mat mat,MatFactorType *t)
2564: {
2568: *t = mat->factortype;
2569: return(0);
2570: }
2572: /* ------------------------------------------------------------*/
2575: /*@C
2576: MatGetInfo - Returns information about matrix storage (number of
2577: nonzeros, memory, etc.).
2579: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used as the flag
2581: Input Parameters:
2582: . mat - the matrix
2584: Output Parameters:
2585: + flag - flag indicating the type of parameters to be returned
2586: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
2587: MAT_GLOBAL_SUM - sum over all processors)
2588: - info - matrix information context
2590: Notes:
2591: The MatInfo context contains a variety of matrix data, including
2592: number of nonzeros allocated and used, number of mallocs during
2593: matrix assembly, etc. Additional information for factored matrices
2594: is provided (such as the fill ratio, number of mallocs during
2595: factorization, etc.). Much of this info is printed to PETSC_STDOUT
2596: when using the runtime options
2597: $ -info -mat_view_info
2599: Example for C/C++ Users:
2600: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
2601: data within the MatInfo context. For example,
2602: .vb
2603: MatInfo info;
2604: Mat A;
2605: double mal, nz_a, nz_u;
2607: MatGetInfo(A,MAT_LOCAL,&info);
2608: mal = info.mallocs;
2609: nz_a = info.nz_allocated;
2610: .ve
2612: Example for Fortran Users:
2613: Fortran users should declare info as a double precision
2614: array of dimension MAT_INFO_SIZE, and then extract the parameters
2615: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
2616: a complete list of parameter names.
2617: .vb
2618: double precision info(MAT_INFO_SIZE)
2619: double precision mal, nz_a
2620: Mat A
2621: integer ierr
2623: call MatGetInfo(A,MAT_LOCAL,info,ierr)
2624: mal = info(MAT_INFO_MALLOCS)
2625: nz_a = info(MAT_INFO_NZ_ALLOCATED)
2626: .ve
2628: Level: intermediate
2630: Concepts: matrices^getting information on
2631:
2632: Developer Note: fortran interface is not autogenerated as the f90
2633: interface defintion cannot be generated correctly [due to MatInfo]
2634:
2635: @*/
2636: PetscErrorCode MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
2637: {
2644: if (!mat->ops->getinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2645: MatPreallocated(mat);
2646: (*mat->ops->getinfo)(mat,flag,info);
2647: return(0);
2648: }
2650: /* ----------------------------------------------------------*/
2654: /*@C
2655: MatLUFactor - Performs in-place LU factorization of matrix.
2657: Collective on Mat
2659: Input Parameters:
2660: + mat - the matrix
2661: . row - row permutation
2662: . col - column permutation
2663: - info - options for factorization, includes
2664: $ fill - expected fill as ratio of original fill.
2665: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2666: $ Run with the option -info to determine an optimal value to use
2668: Notes:
2669: Most users should employ the simplified KSP interface for linear solvers
2670: instead of working directly with matrix algebra routines such as this.
2671: See, e.g., KSPCreate().
2673: This changes the state of the matrix to a factored matrix; it cannot be used
2674: for example with MatSetValues() unless one first calls MatSetUnfactored().
2676: Level: developer
2678: Concepts: matrices^LU factorization
2680: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
2681: MatGetOrdering(), MatSetUnfactored(), MatFactorInfo
2683: Developer Note: fortran interface is not autogenerated as the f90
2684: interface defintion cannot be generated correctly [due to MatFactorInfo]
2686: @*/
2687: PetscErrorCode MatLUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2688: {
2697: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2698: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2699: if (!mat->ops->lufactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2700: MatPreallocated(mat);
2702: PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
2703: (*mat->ops->lufactor)(mat,row,col,info);
2704: PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
2705: PetscObjectStateIncrease((PetscObject)mat);
2706: return(0);
2707: }
2711: /*@C
2712: MatILUFactor - Performs in-place ILU factorization of matrix.
2714: Collective on Mat
2716: Input Parameters:
2717: + mat - the matrix
2718: . row - row permutation
2719: . col - column permutation
2720: - info - structure containing
2721: $ levels - number of levels of fill.
2722: $ expected fill - as ratio of original fill.
2723: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2724: missing diagonal entries)
2726: Notes:
2727: Probably really in-place only when level of fill is zero, otherwise allocates
2728: new space to store factored matrix and deletes previous memory.
2730: Most users should employ the simplified KSP interface for linear solvers
2731: instead of working directly with matrix algebra routines such as this.
2732: See, e.g., KSPCreate().
2734: Level: developer
2736: Concepts: matrices^ILU factorization
2738: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2740: Developer Note: fortran interface is not autogenerated as the f90
2741: interface defintion cannot be generated correctly [due to MatFactorInfo]
2743: @*/
2744: PetscErrorCode MatILUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2745: {
2754: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"matrix must be square");
2755: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2756: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2757: if (!mat->ops->ilufactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2758: MatPreallocated(mat);
2760: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2761: (*mat->ops->ilufactor)(mat,row,col,info);
2762: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2763: PetscObjectStateIncrease((PetscObject)mat);
2764: return(0);
2765: }
2769: /*@C
2770: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
2771: Call this routine before calling MatLUFactorNumeric().
2773: Collective on Mat
2775: Input Parameters:
2776: + fact - the factor matrix obtained with MatGetFactor()
2777: . mat - the matrix
2778: . row, col - row and column permutations
2779: - info - options for factorization, includes
2780: $ fill - expected fill as ratio of original fill.
2781: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2782: $ Run with the option -info to determine an optimal value to use
2785: Notes:
2786: See the <a href="../../docs/manual.pdf">users manual</a> for additional information about
2787: choosing the fill factor for better efficiency.
2789: Most users should employ the simplified KSP interface for linear solvers
2790: instead of working directly with matrix algebra routines such as this.
2791: See, e.g., KSPCreate().
2793: Level: developer
2795: Concepts: matrices^LU symbolic factorization
2797: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2799: Developer Note: fortran interface is not autogenerated as the f90
2800: interface defintion cannot be generated correctly [due to MatFactorInfo]
2802: @*/
2803: PetscErrorCode MatLUFactorSymbolic(Mat fact,Mat mat,IS row,IS col,const MatFactorInfo *info)
2804: {
2814: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2815: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2816: if (!(fact)->ops->lufactorsymbolic) {
2817: const MatSolverPackage spackage;
2818: MatFactorGetSolverPackage(fact,&spackage);
2819: SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Matrix type %s symbolic LU using solver package %s",((PetscObject)mat)->type_name,spackage);
2820: }
2821: MatPreallocated(mat);
2823: PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
2824: (fact->ops->lufactorsymbolic)(fact,mat,row,col,info);
2825: PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
2826: PetscObjectStateIncrease((PetscObject)fact);
2827: return(0);
2828: }
2832: /*@C
2833: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
2834: Call this routine after first calling MatLUFactorSymbolic().
2836: Collective on Mat
2838: Input Parameters:
2839: + fact - the factor matrix obtained with MatGetFactor()
2840: . mat - the matrix
2841: - info - options for factorization
2843: Notes:
2844: See MatLUFactor() for in-place factorization. See
2845: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
2847: Most users should employ the simplified KSP interface for linear solvers
2848: instead of working directly with matrix algebra routines such as this.
2849: See, e.g., KSPCreate().
2851: Level: developer
2853: Concepts: matrices^LU numeric factorization
2855: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
2857: Developer Note: fortran interface is not autogenerated as the f90
2858: interface defintion cannot be generated correctly [due to MatFactorInfo]
2860: @*/
2861: PetscErrorCode MatLUFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2862: {
2870: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2871: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2872: SETERRQ4(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Mat fact: global dimensions are different %D should = %D %D should = %D",mat->rmap->N,(fact)->rmap->N,mat->cmap->N,(fact)->cmap->N);
2873: }
2874: if (!(fact)->ops->lufactornumeric) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s numeric LU",((PetscObject)mat)->type_name);
2875: MatPreallocated(mat);
2876: PetscLogEventBegin(MAT_LUFactorNumeric,mat,fact,0,0);
2877: (fact->ops->lufactornumeric)(fact,mat,info);
2878: PetscLogEventEnd(MAT_LUFactorNumeric,mat,fact,0,0);
2880: MatView_Private(fact);
2881: PetscObjectStateIncrease((PetscObject)fact);
2882: return(0);
2883: }
2887: /*@C
2888: MatCholeskyFactor - Performs in-place Cholesky factorization of a
2889: symmetric matrix.
2891: Collective on Mat
2893: Input Parameters:
2894: + mat - the matrix
2895: . perm - row and column permutations
2896: - f - expected fill as ratio of original fill
2898: Notes:
2899: See MatLUFactor() for the nonsymmetric case. See also
2900: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
2902: Most users should employ the simplified KSP interface for linear solvers
2903: instead of working directly with matrix algebra routines such as this.
2904: See, e.g., KSPCreate().
2906: Level: developer
2908: Concepts: matrices^Cholesky factorization
2910: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
2911: MatGetOrdering()
2913: Developer Note: fortran interface is not autogenerated as the f90
2914: interface defintion cannot be generated correctly [due to MatFactorInfo]
2916: @*/
2917: PetscErrorCode MatCholeskyFactor(Mat mat,IS perm,const MatFactorInfo *info)
2918: {
2926: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"Matrix must be square");
2927: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2928: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2929: if (!mat->ops->choleskyfactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2930: MatPreallocated(mat);
2932: PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
2933: (*mat->ops->choleskyfactor)(mat,perm,info);
2934: PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
2935: PetscObjectStateIncrease((PetscObject)mat);
2936: return(0);
2937: }
2941: /*@C
2942: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
2943: of a symmetric matrix.
2945: Collective on Mat
2947: Input Parameters:
2948: + fact - the factor matrix obtained with MatGetFactor()
2949: . mat - the matrix
2950: . perm - row and column permutations
2951: - info - options for factorization, includes
2952: $ fill - expected fill as ratio of original fill.
2953: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2954: $ Run with the option -info to determine an optimal value to use
2956: Notes:
2957: See MatLUFactorSymbolic() for the nonsymmetric case. See also
2958: MatCholeskyFactor() and MatCholeskyFactorNumeric().
2960: Most users should employ the simplified KSP interface for linear solvers
2961: instead of working directly with matrix algebra routines such as this.
2962: See, e.g., KSPCreate().
2964: Level: developer
2966: Concepts: matrices^Cholesky symbolic factorization
2968: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
2969: MatGetOrdering()
2971: Developer Note: fortran interface is not autogenerated as the f90
2972: interface defintion cannot be generated correctly [due to MatFactorInfo]
2974: @*/
2975: PetscErrorCode MatCholeskyFactorSymbolic(Mat fact,Mat mat,IS perm,const MatFactorInfo *info)
2976: {
2985: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"Matrix must be square");
2986: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2987: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2988: if (!(fact)->ops->choleskyfactorsymbolic) {
2989: const MatSolverPackage spackage;
2990: MatFactorGetSolverPackage(fact,&spackage);
2991: SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s symbolic factor Cholesky using solver package %s",((PetscObject)mat)->type_name,spackage);
2992: }
2993: MatPreallocated(mat);
2995: PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2996: (fact->ops->choleskyfactorsymbolic)(fact,mat,perm,info);
2997: PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2998: PetscObjectStateIncrease((PetscObject)fact);
2999: return(0);
3000: }
3004: /*@C
3005: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
3006: of a symmetric matrix. Call this routine after first calling
3007: MatCholeskyFactorSymbolic().
3009: Collective on Mat
3011: Input Parameters:
3012: + fact - the factor matrix obtained with MatGetFactor()
3013: . mat - the initial matrix
3014: . info - options for factorization
3015: - fact - the symbolic factor of mat
3018: Notes:
3019: Most users should employ the simplified KSP interface for linear solvers
3020: instead of working directly with matrix algebra routines such as this.
3021: See, e.g., KSPCreate().
3023: Level: developer
3025: Concepts: matrices^Cholesky numeric factorization
3027: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
3029: Developer Note: fortran interface is not autogenerated as the f90
3030: interface defintion cannot be generated correctly [due to MatFactorInfo]
3032: @*/
3033: PetscErrorCode MatCholeskyFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
3034: {
3042: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3043: if (!(fact)->ops->choleskyfactornumeric) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s numeric factor Cholesky",((PetscObject)mat)->type_name);
3044: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
3045: SETERRQ4(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Mat fact: global dim %D should = %D %D should = %D",mat->rmap->N,(fact)->rmap->N,mat->cmap->N,(fact)->cmap->N);
3046: }
3047: MatPreallocated(mat);
3049: PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,fact,0,0);
3050: (fact->ops->choleskyfactornumeric)(fact,mat,info);
3051: PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,fact,0,0);
3053: MatView_Private(fact);
3054: PetscObjectStateIncrease((PetscObject)fact);
3055: return(0);
3056: }
3058: /* ----------------------------------------------------------------*/
3061: /*@
3062: MatSolve - Solves A x = b, given a factored matrix.
3064: Neighbor-wise Collective on Mat and Vec
3066: Input Parameters:
3067: + mat - the factored matrix
3068: - b - the right-hand-side vector
3070: Output Parameter:
3071: . x - the result vector
3073: Notes:
3074: The vectors b and x cannot be the same. I.e., one cannot
3075: call MatSolve(A,x,x).
3077: Notes:
3078: Most users should employ the simplified KSP interface for linear solvers
3079: instead of working directly with matrix algebra routines such as this.
3080: See, e.g., KSPCreate().
3082: Level: developer
3084: Concepts: matrices^triangular solves
3086: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
3087: @*/
3088: PetscErrorCode MatSolve(Mat mat,Vec b,Vec x)
3089: {
3099: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3100: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3101: if (mat->cmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3102: if (mat->rmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3103: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3104: if (!mat->rmap->N && !mat->cmap->N) return(0);
3105: if (!mat->ops->solve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3106: MatPreallocated(mat);
3108: PetscLogEventBegin(MAT_Solve,mat,b,x,0);
3109: (*mat->ops->solve)(mat,b,x);
3110: PetscLogEventEnd(MAT_Solve,mat,b,x,0);
3111: PetscObjectStateIncrease((PetscObject)x);
3112: return(0);
3113: }
3117: PetscErrorCode MatMatSolve_Basic(Mat A,Mat B,Mat X)
3118: {
3120: Vec b,x;
3121: PetscInt m,N,i;
3122: PetscScalar *bb,*xx;
3123: PetscBool flg;
3126: PetscTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
3127: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
3128: PetscTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
3129: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
3131: MatGetArray(B,&bb);
3132: MatGetArray(X,&xx);
3133: MatGetLocalSize(B,&m,PETSC_NULL); /* number local rows */
3134: MatGetSize(B,PETSC_NULL,&N); /* total columns in dense matrix */
3135: MatGetVecs(A,&x,&b);
3136: for (i=0; i<N; i++) {
3137: VecPlaceArray(b,bb + i*m);
3138: VecPlaceArray(x,xx + i*m);
3139: MatSolve(A,b,x);
3140: VecResetArray(x);
3141: VecResetArray(b);
3142: }
3143: VecDestroy(&b);
3144: VecDestroy(&x);
3145: MatRestoreArray(B,&bb);
3146: MatRestoreArray(X,&xx);
3147: return(0);
3148: }
3152: /*@
3153: MatMatSolve - Solves A X = B, given a factored matrix.
3155: Neighbor-wise Collective on Mat
3157: Input Parameters:
3158: + mat - the factored matrix
3159: - B - the right-hand-side matrix (dense matrix)
3161: Output Parameter:
3162: . X - the result matrix (dense matrix)
3164: Notes:
3165: The matrices b and x cannot be the same. I.e., one cannot
3166: call MatMatSolve(A,x,x).
3168: Notes:
3169: Most users should usually employ the simplified KSP interface for linear solvers
3170: instead of working directly with matrix algebra routines such as this.
3171: See, e.g., KSPCreate(). However KSP can only solve for one vector (column of X)
3172: at a time.
3174: Level: developer
3176: Concepts: matrices^triangular solves
3178: .seealso: MatMatSolveAdd(), MatMatSolveTranspose(), MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()
3179: @*/
3180: PetscErrorCode MatMatSolve(Mat A,Mat B,Mat X)
3181: {
3191: if (X == B) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_IDN,"X and B must be different matrices");
3192: if (!A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3193: if (A->cmap->N != X->rmap->N) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Mat A,Mat X: global dim %D %D",A->cmap->N,X->rmap->N);
3194: if (A->rmap->N != B->rmap->N) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %D %D",A->rmap->N,B->rmap->N);
3195: if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D",A->rmap->n,B->rmap->n);
3196: if (!A->rmap->N && !A->cmap->N) return(0);
3197: MatPreallocated(A);
3199: PetscLogEventBegin(MAT_MatSolve,A,B,X,0);
3200: if (!A->ops->matsolve) {
3201: PetscInfo1(A,"Mat type %s using basic MatMatSolve",((PetscObject)A)->type_name);
3202: MatMatSolve_Basic(A,B,X);
3203: } else {
3204: (*A->ops->matsolve)(A,B,X);
3205: }
3206: PetscLogEventEnd(MAT_MatSolve,A,B,X,0);
3207: PetscObjectStateIncrease((PetscObject)X);
3208: return(0);
3209: }
3214: /*@
3215: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU, or
3216: U^T*D^(1/2) x = b, given a factored symmetric matrix, A = U^T*D*U,
3218: Neighbor-wise Collective on Mat and Vec
3220: Input Parameters:
3221: + mat - the factored matrix
3222: - b - the right-hand-side vector
3224: Output Parameter:
3225: . x - the result vector
3227: Notes:
3228: MatSolve() should be used for most applications, as it performs
3229: a forward solve followed by a backward solve.
3231: The vectors b and x cannot be the same, i.e., one cannot
3232: call MatForwardSolve(A,x,x).
3234: For matrix in seqsbaij format with block size larger than 1,
3235: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
3236: MatForwardSolve() solves U^T*D y = b, and
3237: MatBackwardSolve() solves U x = y.
3238: Thus they do not provide a symmetric preconditioner.
3240: Most users should employ the simplified KSP interface for linear solvers
3241: instead of working directly with matrix algebra routines such as this.
3242: See, e.g., KSPCreate().
3244: Level: developer
3246: Concepts: matrices^forward solves
3248: .seealso: MatSolve(), MatBackwardSolve()
3249: @*/
3250: PetscErrorCode MatForwardSolve(Mat mat,Vec b,Vec x)
3251: {
3261: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3262: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3263: if (!mat->ops->forwardsolve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3264: if (mat->cmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3265: if (mat->rmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3266: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3267: MatPreallocated(mat);
3268: PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
3269: (*mat->ops->forwardsolve)(mat,b,x);
3270: PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
3271: PetscObjectStateIncrease((PetscObject)x);
3272: return(0);
3273: }
3277: /*@
3278: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
3279: D^(1/2) U x = b, given a factored symmetric matrix, A = U^T*D*U,
3281: Neighbor-wise Collective on Mat and Vec
3283: Input Parameters:
3284: + mat - the factored matrix
3285: - b - the right-hand-side vector
3287: Output Parameter:
3288: . x - the result vector
3290: Notes:
3291: MatSolve() should be used for most applications, as it performs
3292: a forward solve followed by a backward solve.
3294: The vectors b and x cannot be the same. I.e., one cannot
3295: call MatBackwardSolve(A,x,x).
3297: For matrix in seqsbaij format with block size larger than 1,
3298: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
3299: MatForwardSolve() solves U^T*D y = b, and
3300: MatBackwardSolve() solves U x = y.
3301: Thus they do not provide a symmetric preconditioner.
3303: Most users should employ the simplified KSP interface for linear solvers
3304: instead of working directly with matrix algebra routines such as this.
3305: See, e.g., KSPCreate().
3307: Level: developer
3309: Concepts: matrices^backward solves
3311: .seealso: MatSolve(), MatForwardSolve()
3312: @*/
3313: PetscErrorCode MatBackwardSolve(Mat mat,Vec b,Vec x)
3314: {
3324: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3325: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3326: if (!mat->ops->backwardsolve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3327: if (mat->cmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3328: if (mat->rmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3329: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3330: MatPreallocated(mat);
3332: PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
3333: (*mat->ops->backwardsolve)(mat,b,x);
3334: PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
3335: PetscObjectStateIncrease((PetscObject)x);
3336: return(0);
3337: }
3341: /*@
3342: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
3344: Neighbor-wise Collective on Mat and Vec
3346: Input Parameters:
3347: + mat - the factored matrix
3348: . b - the right-hand-side vector
3349: - y - the vector to be added to
3351: Output Parameter:
3352: . x - the result vector
3354: Notes:
3355: The vectors b and x cannot be the same. I.e., one cannot
3356: call MatSolveAdd(A,x,y,x).
3358: Most users should employ the simplified KSP interface for linear solvers
3359: instead of working directly with matrix algebra routines such as this.
3360: See, e.g., KSPCreate().
3362: Level: developer
3364: Concepts: matrices^triangular solves
3366: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
3367: @*/
3368: PetscErrorCode MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
3369: {
3370: PetscScalar one = 1.0;
3371: Vec tmp;
3383: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3384: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3385: if (mat->cmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3386: if (mat->rmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3387: if (mat->rmap->N != y->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
3388: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3389: if (x->map->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->map->n,y->map->n);
3390: MatPreallocated(mat);
3392: PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
3393: if (mat->ops->solveadd) {
3394: (*mat->ops->solveadd)(mat,b,y,x);
3395: } else {
3396: /* do the solve then the add manually */
3397: if (x != y) {
3398: MatSolve(mat,b,x);
3399: VecAXPY(x,one,y);
3400: } else {
3401: VecDuplicate(x,&tmp);
3402: PetscLogObjectParent(mat,tmp);
3403: VecCopy(x,tmp);
3404: MatSolve(mat,b,x);
3405: VecAXPY(x,one,tmp);
3406: VecDestroy(&tmp);
3407: }
3408: }
3409: PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
3410: PetscObjectStateIncrease((PetscObject)x);
3411: return(0);
3412: }
3416: /*@
3417: MatSolveTranspose - Solves A' x = b, given a factored matrix.
3419: Neighbor-wise Collective on Mat and Vec
3421: Input Parameters:
3422: + mat - the factored matrix
3423: - b - the right-hand-side vector
3425: Output Parameter:
3426: . x - the result vector
3428: Notes:
3429: The vectors b and x cannot be the same. I.e., one cannot
3430: call MatSolveTranspose(A,x,x).
3432: Most users should employ the simplified KSP interface for linear solvers
3433: instead of working directly with matrix algebra routines such as this.
3434: See, e.g., KSPCreate().
3436: Level: developer
3438: Concepts: matrices^triangular solves
3440: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
3441: @*/
3442: PetscErrorCode MatSolveTranspose(Mat mat,Vec b,Vec x)
3443: {
3453: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3454: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3455: if (!mat->ops->solvetranspose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Matrix type %s",((PetscObject)mat)->type_name);
3456: if (mat->rmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
3457: if (mat->cmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->cmap->N,b->map->N);
3458: MatPreallocated(mat);
3459: PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
3460: (*mat->ops->solvetranspose)(mat,b,x);
3461: PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
3462: PetscObjectStateIncrease((PetscObject)x);
3463: return(0);
3464: }
3468: /*@
3469: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
3470: factored matrix.
3472: Neighbor-wise Collective on Mat and Vec
3474: Input Parameters:
3475: + mat - the factored matrix
3476: . b - the right-hand-side vector
3477: - y - the vector to be added to
3479: Output Parameter:
3480: . x - the result vector
3482: Notes:
3483: The vectors b and x cannot be the same. I.e., one cannot
3484: call MatSolveTransposeAdd(A,x,y,x).
3486: Most users should employ the simplified KSP interface for linear solvers
3487: instead of working directly with matrix algebra routines such as this.
3488: See, e.g., KSPCreate().
3490: Level: developer
3492: Concepts: matrices^triangular solves
3494: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
3495: @*/
3496: PetscErrorCode MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
3497: {
3498: PetscScalar one = 1.0;
3500: Vec tmp;
3511: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3512: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3513: if (mat->rmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
3514: if (mat->cmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->cmap->N,b->map->N);
3515: if (mat->cmap->N != y->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->cmap->N,y->map->N);
3516: if (x->map->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->map->n,y->map->n);
3517: MatPreallocated(mat);
3519: PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
3520: if (mat->ops->solvetransposeadd) {
3521: (*mat->ops->solvetransposeadd)(mat,b,y,x);
3522: } else {
3523: /* do the solve then the add manually */
3524: if (x != y) {
3525: MatSolveTranspose(mat,b,x);
3526: VecAXPY(x,one,y);
3527: } else {
3528: VecDuplicate(x,&tmp);
3529: PetscLogObjectParent(mat,tmp);
3530: VecCopy(x,tmp);
3531: MatSolveTranspose(mat,b,x);
3532: VecAXPY(x,one,tmp);
3533: VecDestroy(&tmp);
3534: }
3535: }
3536: PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
3537: PetscObjectStateIncrease((PetscObject)x);
3538: return(0);
3539: }
3540: /* ----------------------------------------------------------------*/
3544: /*@
3545: MatSOR - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3547: Neighbor-wise Collective on Mat and Vec
3549: Input Parameters:
3550: + mat - the matrix
3551: . b - the right hand side
3552: . omega - the relaxation factor
3553: . flag - flag indicating the type of SOR (see below)
3554: . shift - diagonal shift
3555: . its - the number of iterations
3556: - lits - the number of local iterations
3558: Output Parameters:
3559: . x - the solution (can contain an initial guess, use option SOR_ZERO_INITIAL_GUESS to indicate no guess)
3561: SOR Flags:
3562: . SOR_FORWARD_SWEEP - forward SOR
3563: . SOR_BACKWARD_SWEEP - backward SOR
3564: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
3565: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
3566: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
3567: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
3568: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
3569: upper/lower triangular part of matrix to
3570: vector (with omega)
3571: . SOR_ZERO_INITIAL_GUESS - zero initial guess
3573: Notes:
3574: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
3575: SOR_LOCAL_SYMMETRIC_SWEEP perform separate independent smoothings
3576: on each processor.
3578: Application programmers will not generally use MatSOR() directly,
3579: but instead will employ the KSP/PC interface.
3581: Notes: for BAIJ, SBAIJ, and AIJ matrices with Inodes this does a block SOR smoothing, otherwise it does a pointwise smoothing
3583: Notes for Advanced Users:
3584: The flags are implemented as bitwise inclusive or operations.
3585: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
3586: to specify a zero initial guess for SSOR.
3588: Most users should employ the simplified KSP interface for linear solvers
3589: instead of working directly with matrix algebra routines such as this.
3590: See, e.g., KSPCreate().
3592: Vectors x and b CANNOT be the same
3594: Developer Note: We should add block SOR support for AIJ matrices with block size set to great than one and no inodes
3596: Level: developer
3598: Concepts: matrices^relaxation
3599: Concepts: matrices^SOR
3600: Concepts: matrices^Gauss-Seidel
3602: @*/
3603: PetscErrorCode MatSOR(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3604: {
3614: if (!mat->ops->sor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3615: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3616: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3617: if (mat->cmap->N != x->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3618: if (mat->rmap->N != b->map->N) SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3619: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3620: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
3621: if (lits <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires local its %D positive",lits);
3622: if (b == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"b and x vector cannot be the same");
3624: MatPreallocated(mat);
3625: PetscLogEventBegin(MAT_SOR,mat,b,x,0);
3626: ierr =(*mat->ops->sor)(mat,b,omega,flag,shift,its,lits,x);
3627: PetscLogEventEnd(MAT_SOR,mat,b,x,0);
3628: PetscObjectStateIncrease((PetscObject)x);
3629: return(0);
3630: }
3634: /*
3635: Default matrix copy routine.
3636: */
3637: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
3638: {
3639: PetscErrorCode ierr;
3640: PetscInt i,rstart = 0,rend = 0,nz;
3641: const PetscInt *cwork;
3642: const PetscScalar *vwork;
3645: if (B->assembled) {
3646: MatZeroEntries(B);
3647: }
3648: MatGetOwnershipRange(A,&rstart,&rend);
3649: for (i=rstart; i<rend; i++) {
3650: MatGetRow(A,i,&nz,&cwork,&vwork);
3651: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
3652: MatRestoreRow(A,i,&nz,&cwork,&vwork);
3653: }
3654: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3655: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3656: PetscObjectStateIncrease((PetscObject)B);
3657: return(0);
3658: }
3662: /*@
3663: MatCopy - Copys a matrix to another matrix.
3665: Collective on Mat
3667: Input Parameters:
3668: + A - the matrix
3669: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
3671: Output Parameter:
3672: . B - where the copy is put
3674: Notes:
3675: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
3676: same nonzero pattern or the routine will crash.
3678: MatCopy() copies the matrix entries of a matrix to another existing
3679: matrix (after first zeroing the second matrix). A related routine is
3680: MatConvert(), which first creates a new matrix and then copies the data.
3682: Level: intermediate
3683:
3684: Concepts: matrices^copying
3686: .seealso: MatConvert(), MatDuplicate()
3688: @*/
3689: PetscErrorCode MatCopy(Mat A,Mat B,MatStructure str)
3690: {
3692: PetscInt i;
3700: MatPreallocated(B);
3701: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3702: if (A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3703: if (A->rmap->N != B->rmap->N || A->cmap->N != B->cmap->N) SETERRQ4(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim (%D,%D) (%D,%D)",A->rmap->N,B->rmap->N,A->cmap->N,B->cmap->N);
3704: MatPreallocated(A);
3706: PetscLogEventBegin(MAT_Copy,A,B,0,0);
3707: if (A->ops->copy) {
3708: (*A->ops->copy)(A,B,str);
3709: } else { /* generic conversion */
3710: MatCopy_Basic(A,B,str);
3711: }
3713: B->stencil.dim = A->stencil.dim;
3714: B->stencil.noc = A->stencil.noc;
3715: for (i=0; i<=A->stencil.dim; i++) {
3716: B->stencil.dims[i] = A->stencil.dims[i];
3717: B->stencil.starts[i] = A->stencil.starts[i];
3718: }
3720: PetscLogEventEnd(MAT_Copy,A,B,0,0);
3721: PetscObjectStateIncrease((PetscObject)B);
3722: return(0);
3723: }
3727: /*@C
3728: MatConvert - Converts a matrix to another matrix, either of the same
3729: or different type.
3731: Collective on Mat
3733: Input Parameters:
3734: + mat - the matrix
3735: . newtype - new matrix type. Use MATSAME to create a new matrix of the
3736: same type as the original matrix.
3737: - reuse - denotes if the destination matrix is to be created or reused. Currently
3738: MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
3739: MAT_INITIAL_MATRIX.
3741: Output Parameter:
3742: . M - pointer to place new matrix
3744: Notes:
3745: MatConvert() first creates a new matrix and then copies the data from
3746: the first matrix. A related routine is MatCopy(), which copies the matrix
3747: entries of one matrix to another already existing matrix context.
3749: Cannot be used to convert a sequential matrix to parallel or parallel to sequential,
3750: the MPI communicator of the generated matrix is always the same as the communicator
3751: of the input matrix.
3753: Level: intermediate
3755: Concepts: matrices^converting between storage formats
3757: .seealso: MatCopy(), MatDuplicate()
3758: @*/
3759: PetscErrorCode MatConvert(Mat mat, const MatType newtype,MatReuse reuse,Mat *M)
3760: {
3761: PetscErrorCode ierr;
3762: PetscBool sametype,issame,flg;
3763: char convname[256],mtype[256];
3764: Mat B;
3770: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3771: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3772: MatPreallocated(mat);
3774: PetscOptionsGetString(((PetscObject)mat)->prefix,"-matconvert_type",mtype,256,&flg);
3775: if (flg) {
3776: newtype = mtype;
3777: }
3778: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
3779: PetscStrcmp(newtype,"same",&issame);
3780: if ((reuse == MAT_REUSE_MATRIX) && (mat != *M)) {
3781: SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for in-place conversion currently");
3782: }
3784: if ((reuse == MAT_REUSE_MATRIX) && (issame || sametype)) return(0);
3785:
3786: if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
3787: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
3788: } else {
3789: PetscErrorCode (*conv)(Mat, const MatType,MatReuse,Mat*)=PETSC_NULL;
3790: const char *prefix[3] = {"seq","mpi",""};
3791: PetscInt i;
3792: /*
3793: Order of precedence:
3794: 1) See if a specialized converter is known to the current matrix.
3795: 2) See if a specialized converter is known to the desired matrix class.
3796: 3) See if a good general converter is registered for the desired class
3797: (as of 6/27/03 only MATMPIADJ falls into this category).
3798: 4) See if a good general converter is known for the current matrix.
3799: 5) Use a really basic converter.
3800: */
3801:
3802: /* 1) See if a specialized converter is known to the current matrix and the desired class */
3803: for (i=0; i<3; i++) {
3804: PetscStrcpy(convname,"MatConvert_");
3805: PetscStrcat(convname,((PetscObject)mat)->type_name);
3806: PetscStrcat(convname,"_");
3807: PetscStrcat(convname,prefix[i]);
3808: PetscStrcat(convname,issame?((PetscObject)mat)->type_name:newtype);
3809: PetscStrcat(convname,"_C");
3810: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3811: if (conv) goto foundconv;
3812: }
3814: /* 2) See if a specialized converter is known to the desired matrix class. */
3815: MatCreate(((PetscObject)mat)->comm,&B);
3816: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N);
3817: MatSetType(B,newtype);
3818: for (i=0; i<3; i++) {
3819: PetscStrcpy(convname,"MatConvert_");
3820: PetscStrcat(convname,((PetscObject)mat)->type_name);
3821: PetscStrcat(convname,"_");
3822: PetscStrcat(convname,prefix[i]);
3823: PetscStrcat(convname,newtype);
3824: PetscStrcat(convname,"_C");
3825: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
3826: if (conv) {
3827: MatDestroy(&B);
3828: goto foundconv;
3829: }
3830: }
3832: /* 3) See if a good general converter is registered for the desired class */
3833: conv = B->ops->convertfrom;
3834: MatDestroy(&B);
3835: if (conv) goto foundconv;
3837: /* 4) See if a good general converter is known for the current matrix */
3838: if (mat->ops->convert) {
3839: conv = mat->ops->convert;
3840: }
3841: if (conv) goto foundconv;
3843: /* 5) Use a really basic converter. */
3844: conv = MatConvert_Basic;
3846: foundconv:
3847: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3848: (*conv)(mat,newtype,reuse,M);
3849: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3850: }
3851: PetscObjectStateIncrease((PetscObject)*M);
3853: /* Copy Mat options */
3854: if (mat->symmetric){MatSetOption(*M,MAT_SYMMETRIC,PETSC_TRUE);}
3855: if (mat->hermitian){MatSetOption(*M,MAT_HERMITIAN,PETSC_TRUE);}
3856: return(0);
3857: }
3861: /*@C
3862: MatFactorGetSolverPackage - Returns name of the package providing the factorization routines
3864: Not Collective
3866: Input Parameter:
3867: . mat - the matrix, must be a factored matrix
3869: Output Parameter:
3870: . type - the string name of the package (do not free this string)
3872: Notes:
3873: In Fortran you pass in a empty string and the package name will be copied into it.
3874: (Make sure the string is long enough)
3876: Level: intermediate
3878: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable(), MatGetFactor()
3879: @*/
3880: PetscErrorCode MatFactorGetSolverPackage(Mat mat, const MatSolverPackage *type)
3881: {
3882: PetscErrorCode ierr;
3883: PetscErrorCode (*conv)(Mat,const MatSolverPackage*);
3888: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Only for factored matrix");
3889: PetscObjectQueryFunction((PetscObject)mat,"MatFactorGetSolverPackage_C",(void (**)(void))&conv);
3890: if (!conv) {
3891: *type = MATSOLVERPETSC;
3892: } else {
3893: (*conv)(mat,type);
3894: }
3895: return(0);
3896: }
3900: /*@C
3901: MatGetFactor - Returns a matrix suitable to calls to MatXXFactorSymbolic()
3903: Collective on Mat
3905: Input Parameters:
3906: + mat - the matrix
3907: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3908: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3910: Output Parameters:
3911: . f - the factor matrix used with MatXXFactorSymbolic() calls
3913: Notes:
3914: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3915: such as pastix, superlu, mumps, spooles etc.
3917: PETSc must have been ./configure to use the external solver, using the option --download-package
3919: Level: intermediate
3921: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable()
3922: @*/
3923: PetscErrorCode MatGetFactor(Mat mat, const MatSolverPackage type,MatFactorType ftype,Mat *f)
3924: {
3925: PetscErrorCode ierr,(*conv)(Mat,MatFactorType,Mat*);
3926: char convname[256];
3932: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3933: MatPreallocated(mat);
3935: PetscStrcpy(convname,"MatGetFactor_");
3936: PetscStrcat(convname,type);
3937: PetscStrcat(convname,"_C");
3938: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3939: if (!conv) {
3940: PetscBool flag;
3941: MPI_Comm comm;
3943: PetscObjectGetComm((PetscObject)mat,&comm);
3944: PetscStrcasecmp(MATSOLVERPETSC,type,&flag);
3945: if (flag) {
3946: SETERRQ2(comm,PETSC_ERR_SUP,"Matrix format %s does not have a built-in PETSc %s",((PetscObject)mat)->type_name,MatFactorTypes[ftype]);
3947: } else {
3948: SETERRQ4(comm,PETSC_ERR_SUP,"Matrix format %s does not have a solver package %s for %s. Perhaps you must ./configure with --download-%s",((PetscObject)mat)->type_name,type,MatFactorTypes[ftype],type);
3949: }
3950: }
3951: (*conv)(mat,ftype,f);
3952: return(0);
3953: }
3957: /*@C
3958: MatGetFactorAvailable - Returns a a flag if matrix supports particular package and factor type
3960: Not Collective
3962: Input Parameters:
3963: + mat - the matrix
3964: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3965: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3967: Output Parameter:
3968: . flg - PETSC_TRUE if the factorization is available
3970: Notes:
3971: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3972: such as pastix, superlu, mumps, spooles etc.
3974: PETSc must have been ./configure to use the external solver, using the option --download-package
3976: Level: intermediate
3978: .seealso: MatCopy(), MatDuplicate(), MatGetFactor()
3979: @*/
3980: PetscErrorCode MatGetFactorAvailable(Mat mat, const MatSolverPackage type,MatFactorType ftype,PetscBool *flg)
3981: {
3982: PetscErrorCode ierr;
3983: char convname[256];
3984: PetscErrorCode (*conv)(Mat,MatFactorType,PetscBool *);
3990: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3991: MatPreallocated(mat);
3993: PetscStrcpy(convname,"MatGetFactorAvailable_");
3994: PetscStrcat(convname,type);
3995: PetscStrcat(convname,"_C");
3996: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3997: if (!conv) {
3998: *flg = PETSC_FALSE;
3999: } else {
4000: (*conv)(mat,ftype,flg);
4001: }
4002: return(0);
4003: }
4008: /*@
4009: MatDuplicate - Duplicates a matrix including the non-zero structure.
4011: Collective on Mat
4013: Input Parameters:
4014: + mat - the matrix
4015: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy the numerical values in the matrix
4016: MAT_SHARE_NONZERO_PATTERN to share the nonzero patterns with the previous matrix and not copy them.
4018: Output Parameter:
4019: . M - pointer to place new matrix
4021: Level: intermediate
4023: Concepts: matrices^duplicating
4025: Notes: You cannot change the nonzero pattern for the parent or child matrix if you use MAT_SHARE_NONZERO_PATTERN.
4027: .seealso: MatCopy(), MatConvert()
4028: @*/
4029: PetscErrorCode MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
4030: {
4032: Mat B;
4033: PetscInt i;
4039: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4040: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4041: MatPreallocated(mat);
4043: *M = 0;
4044: if (!mat->ops->duplicate) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Not written for this matrix type");
4045: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
4046: (*mat->ops->duplicate)(mat,op,M);
4047: B = *M;
4048:
4049: B->stencil.dim = mat->stencil.dim;
4050: B->stencil.noc = mat->stencil.noc;
4051: for (i=0; i<=mat->stencil.dim; i++) {
4052: B->stencil.dims[i] = mat->stencil.dims[i];
4053: B->stencil.starts[i] = mat->stencil.starts[i];
4054: }
4056: B->nooffproczerorows = mat->nooffproczerorows;
4057: B->nooffprocentries = mat->nooffprocentries;
4058: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
4059: PetscObjectStateIncrease((PetscObject)B);
4060: return(0);
4061: }
4065: /*@
4066: MatGetDiagonal - Gets the diagonal of a matrix.
4068: Logically Collective on Mat and Vec
4070: Input Parameters:
4071: + mat - the matrix
4072: - v - the vector for storing the diagonal
4074: Output Parameter:
4075: . v - the diagonal of the matrix
4077: Level: intermediate
4079: Concepts: matrices^accessing diagonals
4081: .seealso: MatGetRow(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs()
4082: @*/
4083: PetscErrorCode MatGetDiagonal(Mat mat,Vec v)
4084: {
4091: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4092: if (!mat->ops->getdiagonal) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4093: MatPreallocated(mat);
4095: (*mat->ops->getdiagonal)(mat,v);
4096: PetscObjectStateIncrease((PetscObject)v);
4097: return(0);
4098: }
4102: /*@
4103: MatGetRowMin - Gets the minimum value (of the real part) of each
4104: row of the matrix
4106: Logically Collective on Mat and Vec
4108: Input Parameters:
4109: . mat - the matrix
4111: Output Parameter:
4112: + v - the vector for storing the maximums
4113: - idx - the indices of the column found for each row (optional)
4115: Level: intermediate
4117: Notes: The result of this call are the same as if one converted the matrix to dense format
4118: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
4120: This code is only implemented for a couple of matrix formats.
4122: Concepts: matrices^getting row maximums
4124: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(),
4125: MatGetRowMax()
4126: @*/
4127: PetscErrorCode MatGetRowMin(Mat mat,Vec v,PetscInt idx[])
4128: {
4135: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4136: if (!mat->ops->getrowmax) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4137: MatPreallocated(mat);
4139: (*mat->ops->getrowmin)(mat,v,idx);
4140: PetscObjectStateIncrease((PetscObject)v);
4141: return(0);
4142: }
4146: /*@
4147: MatGetRowMinAbs - Gets the minimum value (in absolute value) of each
4148: row of the matrix
4150: Logically Collective on Mat and Vec
4152: Input Parameters:
4153: . mat - the matrix
4155: Output Parameter:
4156: + v - the vector for storing the minimums
4157: - idx - the indices of the column found for each row (optional)
4159: Level: intermediate
4161: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
4162: row is 0 (the first column).
4164: This code is only implemented for a couple of matrix formats.
4166: Concepts: matrices^getting row maximums
4168: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMaxAbs(), MatGetRowMin()
4169: @*/
4170: PetscErrorCode MatGetRowMinAbs(Mat mat,Vec v,PetscInt idx[])
4171: {
4178: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4179: if (!mat->ops->getrowminabs) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4180: MatPreallocated(mat);
4181: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
4183: (*mat->ops->getrowminabs)(mat,v,idx);
4184: PetscObjectStateIncrease((PetscObject)v);
4185: return(0);
4186: }
4190: /*@
4191: MatGetRowMax - Gets the maximum value (of the real part) of each
4192: row of the matrix
4194: Logically Collective on Mat and Vec
4196: Input Parameters:
4197: . mat - the matrix
4199: Output Parameter:
4200: + v - the vector for storing the maximums
4201: - idx - the indices of the column found for each row (optional)
4203: Level: intermediate
4205: Notes: The result of this call are the same as if one converted the matrix to dense format
4206: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
4208: This code is only implemented for a couple of matrix formats.
4210: Concepts: matrices^getting row maximums
4212: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(), MatGetRowMin()
4213: @*/
4214: PetscErrorCode MatGetRowMax(Mat mat,Vec v,PetscInt idx[])
4215: {
4222: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4223: if (!mat->ops->getrowmax) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4224: MatPreallocated(mat);
4226: (*mat->ops->getrowmax)(mat,v,idx);
4227: PetscObjectStateIncrease((PetscObject)v);
4228: return(0);
4229: }
4233: /*@
4234: MatGetRowMaxAbs - Gets the maximum value (in absolute value) of each
4235: row of the matrix
4237: Logically Collective on Mat and Vec
4239: Input Parameters:
4240: . mat - the matrix
4242: Output Parameter:
4243: + v - the vector for storing the maximums
4244: - idx - the indices of the column found for each row (optional)
4246: Level: intermediate
4248: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
4249: row is 0 (the first column).
4251: This code is only implemented for a couple of matrix formats.
4253: Concepts: matrices^getting row maximums
4255: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
4256: @*/
4257: PetscErrorCode MatGetRowMaxAbs(Mat mat,Vec v,PetscInt idx[])
4258: {
4265: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4266: if (!mat->ops->getrowmaxabs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4267: MatPreallocated(mat);
4268: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
4270: (*mat->ops->getrowmaxabs)(mat,v,idx);
4271: PetscObjectStateIncrease((PetscObject)v);
4272: return(0);
4273: }
4277: /*@
4278: MatGetRowSum - Gets the sum of each row of the matrix
4280: Logically Collective on Mat and Vec
4282: Input Parameters:
4283: . mat - the matrix
4285: Output Parameter:
4286: . v - the vector for storing the sum of rows
4288: Level: intermediate
4290: Notes: This code is slow since it is not currently specialized for different formats
4292: Concepts: matrices^getting row sums
4294: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
4295: @*/
4296: PetscErrorCode MatGetRowSum(Mat mat, Vec v)
4297: {
4298: PetscInt start = 0, end = 0, row;
4299: PetscScalar *array;
4306: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4307: MatPreallocated(mat);
4308: MatGetOwnershipRange(mat, &start, &end);
4309: VecGetArray(v, &array);
4310: for(row = start; row < end; ++row) {
4311: PetscInt ncols, col;
4312: const PetscInt *cols;
4313: const PetscScalar *vals;
4315: array[row - start] = 0.0;
4316: MatGetRow(mat, row, &ncols, &cols, &vals);
4317: for(col = 0; col < ncols; col++) {
4318: array[row - start] += vals[col];
4319: }
4320: MatRestoreRow(mat, row, &ncols, &cols, &vals);
4321: }
4322: VecRestoreArray(v, &array);
4323: PetscObjectStateIncrease((PetscObject) v);
4324: return(0);
4325: }
4329: /*@
4330: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
4332: Collective on Mat
4334: Input Parameter:
4335: + mat - the matrix to transpose
4336: - reuse - store the transpose matrix in the provided B
4338: Output Parameters:
4339: . B - the transpose
4341: Notes:
4342: If you pass in &mat for B the transpose will be done in place, for example MatTranspose(mat,MAT_REUSE_MATRIX,&mat);
4344: Level: intermediate
4346: Concepts: matrices^transposing
4348: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose(), MatReuse
4349: @*/
4350: PetscErrorCode MatTranspose(Mat mat,MatReuse reuse,Mat *B)
4351: {
4357: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4358: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4359: if (!mat->ops->transpose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4360: MatPreallocated(mat);
4362: PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
4363: (*mat->ops->transpose)(mat,reuse,B);
4364: PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
4365: if (B) {PetscObjectStateIncrease((PetscObject)*B);}
4366: return(0);
4367: }
4371: /*@
4372: MatIsTranspose - Test whether a matrix is another one's transpose,
4373: or its own, in which case it tests symmetry.
4375: Collective on Mat
4377: Input Parameter:
4378: + A - the matrix to test
4379: - B - the matrix to test against, this can equal the first parameter
4381: Output Parameters:
4382: . flg - the result
4384: Notes:
4385: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
4386: has a running time of the order of the number of nonzeros; the parallel
4387: test involves parallel copies of the block-offdiagonal parts of the matrix.
4389: Level: intermediate
4391: Concepts: matrices^transposing, matrix^symmetry
4393: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
4394: @*/
4395: PetscErrorCode MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscBool *flg)
4396: {
4397: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscBool *),(*g)(Mat,Mat,PetscReal,PetscBool *);
4403: PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
4404: PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
4405: *flg = PETSC_FALSE;
4406: if (f && g) {
4407: if (f == g) {
4408: (*f)(A,B,tol,flg);
4409: } else {
4410: SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
4411: }
4412: } else {
4413: const MatType mattype;
4414: if (!f) {MatGetType(A,&mattype);}
4415: else {MatGetType(B,&mattype);}
4416: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix of type <%s> does not support checking for transpose",mattype);
4417: }
4418: return(0);
4419: }
4423: /*@
4424: MatHermitianTranspose - Computes an in-place or out-of-place transpose of a matrix in complex conjugate.
4426: Collective on Mat
4428: Input Parameter:
4429: + mat - the matrix to transpose and complex conjugate
4430: - reuse - store the transpose matrix in the provided B
4432: Output Parameters:
4433: . B - the Hermitian
4435: Notes:
4436: If you pass in &mat for B the Hermitian will be done in place
4438: Level: intermediate
4440: Concepts: matrices^transposing, complex conjugatex
4442: .seealso: MatTranspose(), MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose(), MatReuse
4443: @*/
4444: PetscErrorCode MatHermitianTranspose(Mat mat,MatReuse reuse,Mat *B)
4445: {
4449: MatTranspose(mat,reuse,B);
4450: #if defined(PETSC_USE_COMPLEX)
4451: MatConjugate(*B);
4452: #endif
4453: return(0);
4454: }
4458: /*@
4459: MatIsHermitianTranspose - Test whether a matrix is another one's Hermitian transpose,
4461: Collective on Mat
4463: Input Parameter:
4464: + A - the matrix to test
4465: - B - the matrix to test against, this can equal the first parameter
4467: Output Parameters:
4468: . flg - the result
4470: Notes:
4471: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
4472: has a running time of the order of the number of nonzeros; the parallel
4473: test involves parallel copies of the block-offdiagonal parts of the matrix.
4475: Level: intermediate
4477: Concepts: matrices^transposing, matrix^symmetry
4479: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian(), MatIsTranspose()
4480: @*/
4481: PetscErrorCode MatIsHermitianTranspose(Mat A,Mat B,PetscReal tol,PetscBool *flg)
4482: {
4483: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscBool *),(*g)(Mat,Mat,PetscReal,PetscBool *);
4489: PetscObjectQueryFunction((PetscObject)A,"MatIsHermitianTranspose_C",(void (**)(void))&f);
4490: PetscObjectQueryFunction((PetscObject)B,"MatIsHermitianTranspose_C",(void (**)(void))&g);
4491: if (f && g) {
4492: if (f==g) {
4493: (*f)(A,B,tol,flg);
4494: } else {
4495: SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for Hermitian test");
4496: }
4497: }
4498: return(0);
4499: }
4503: /*@
4504: MatPermute - Creates a new matrix with rows and columns permuted from the
4505: original.
4507: Collective on Mat
4509: Input Parameters:
4510: + mat - the matrix to permute
4511: . row - row permutation, each processor supplies only the permutation for its rows
4512: - col - column permutation, each processor needs the entire column permutation, that is
4513: this is the same size as the total number of columns in the matrix. It can often
4514: be obtained with ISAllGather() on the row permutation
4516: Output Parameters:
4517: . B - the permuted matrix
4519: Level: advanced
4521: Concepts: matrices^permuting
4523: .seealso: MatGetOrdering(), ISAllGather()
4525: @*/
4526: PetscErrorCode MatPermute(Mat mat,IS row,IS col,Mat *B)
4527: {
4536: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4537: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4538: if (!mat->ops->permute) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatPermute not available for Mat type %s",((PetscObject)mat)->type_name);
4539: MatPreallocated(mat);
4541: (*mat->ops->permute)(mat,row,col,B);
4542: PetscObjectStateIncrease((PetscObject)*B);
4543: return(0);
4544: }
4548: /*@
4549: MatEqual - Compares two matrices.
4551: Collective on Mat
4553: Input Parameters:
4554: + A - the first matrix
4555: - B - the second matrix
4557: Output Parameter:
4558: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
4560: Level: intermediate
4562: Concepts: matrices^equality between
4563: @*/
4564: PetscErrorCode MatEqual(Mat A,Mat B,PetscBool *flg)
4565: {
4575: MatPreallocated(B);
4576: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4577: if (!B->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4578: if (A->rmap->N != B->rmap->N || A->cmap->N != B->cmap->N) SETERRQ4(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %D %D %D %D",A->rmap->N,B->rmap->N,A->cmap->N,B->cmap->N);
4579: if (!A->ops->equal) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)A)->type_name);
4580: if (!B->ops->equal) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)B)->type_name);
4581: if (A->ops->equal != B->ops->equal) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_INCOMP,"A is type: %s\nB is type: %s",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
4582: MatPreallocated(A);
4584: (*A->ops->equal)(A,B,flg);
4585: return(0);
4586: }
4590: /*@
4591: MatDiagonalScale - Scales a matrix on the left and right by diagonal
4592: matrices that are stored as vectors. Either of the two scaling
4593: matrices can be PETSC_NULL.
4595: Collective on Mat
4597: Input Parameters:
4598: + mat - the matrix to be scaled
4599: . l - the left scaling vector (or PETSC_NULL)
4600: - r - the right scaling vector (or PETSC_NULL)
4602: Notes:
4603: MatDiagonalScale() computes A = LAR, where
4604: L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)
4605: The L scales the rows of the matrix, the R scales the columns of the matrix.
4607: Level: intermediate
4609: Concepts: matrices^diagonal scaling
4610: Concepts: diagonal scaling of matrices
4612: .seealso: MatScale()
4613: @*/
4614: PetscErrorCode MatDiagonalScale(Mat mat,Vec l,Vec r)
4615: {
4621: if (!mat->ops->diagonalscale) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4624: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4625: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4626: MatPreallocated(mat);
4628: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4629: (*mat->ops->diagonalscale)(mat,l,r);
4630: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4631: PetscObjectStateIncrease((PetscObject)mat);
4632: #if defined(PETSC_HAVE_CUSP)
4633: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
4634: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
4635: }
4636: #endif
4637: return(0);
4638: }
4642: /*@
4643: MatScale - Scales all elements of a matrix by a given number.
4645: Logically Collective on Mat
4647: Input Parameters:
4648: + mat - the matrix to be scaled
4649: - a - the scaling value
4651: Output Parameter:
4652: . mat - the scaled matrix
4654: Level: intermediate
4656: Concepts: matrices^scaling all entries
4658: .seealso: MatDiagonalScale()
4659: @*/
4660: PetscErrorCode MatScale(Mat mat,PetscScalar a)
4661: {
4667: if (a != (PetscScalar)1.0 && !mat->ops->scale) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4668: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4669: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4671: MatPreallocated(mat);
4673: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4674: if (a != (PetscScalar)1.0) {
4675: (*mat->ops->scale)(mat,a);
4676: PetscObjectStateIncrease((PetscObject)mat);
4677: }
4678: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4679: #if defined(PETSC_HAVE_CUSP)
4680: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
4681: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
4682: }
4683: #endif
4684: return(0);
4685: }
4689: /*@
4690: MatNorm - Calculates various norms of a matrix.
4692: Collective on Mat
4694: Input Parameters:
4695: + mat - the matrix
4696: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
4698: Output Parameters:
4699: . nrm - the resulting norm
4701: Level: intermediate
4703: Concepts: matrices^norm
4704: Concepts: norm^of matrix
4705: @*/
4706: PetscErrorCode MatNorm(Mat mat,NormType type,PetscReal *nrm)
4707: {
4715: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4716: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4717: if (!mat->ops->norm) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4718: MatPreallocated(mat);
4720: (*mat->ops->norm)(mat,type,nrm);
4721: return(0);
4722: }
4724: /*
4725: This variable is used to prevent counting of MatAssemblyBegin() that
4726: are called from within a MatAssemblyEnd().
4727: */
4728: static PetscInt MatAssemblyEnd_InUse = 0;
4731: /*@
4732: MatAssemblyBegin - Begins assembling the matrix. This routine should
4733: be called after completing all calls to MatSetValues().
4735: Collective on Mat
4737: Input Parameters:
4738: + mat - the matrix
4739: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4740:
4741: Notes:
4742: MatSetValues() generally caches the values. The matrix is ready to
4743: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4744: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4745: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4746: using the matrix.
4748: Level: beginner
4750: Concepts: matrices^assembling
4752: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
4753: @*/
4754: PetscErrorCode MatAssemblyBegin(Mat mat,MatAssemblyType type)
4755: {
4761: MatPreallocated(mat);
4762: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
4763: if (mat->assembled) {
4764: mat->was_assembled = PETSC_TRUE;
4765: mat->assembled = PETSC_FALSE;
4766: }
4767: if (!MatAssemblyEnd_InUse) {
4768: PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
4769: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4770: PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
4771: } else {
4772: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4773: }
4774: return(0);
4775: }
4779: /*@
4780: MatAssembled - Indicates if a matrix has been assembled and is ready for
4781: use; for example, in matrix-vector product.
4783: Not Collective
4785: Input Parameter:
4786: . mat - the matrix
4788: Output Parameter:
4789: . assembled - PETSC_TRUE or PETSC_FALSE
4791: Level: advanced
4793: Concepts: matrices^assembled?
4795: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
4796: @*/
4797: PetscErrorCode MatAssembled(Mat mat,PetscBool *assembled)
4798: {
4803: *assembled = mat->assembled;
4804: return(0);
4805: }
4809: /*
4810: Processes command line options to determine if/how a matrix
4811: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
4812: */
4813: PetscErrorCode MatView_Private(Mat mat)
4814: {
4815: PetscErrorCode ierr;
4816: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,flg3 = PETSC_FALSE,flg4 = PETSC_FALSE,flg6 = PETSC_FALSE,flg7 = PETSC_FALSE,flg8 = PETSC_FALSE;
4817: static PetscBool incall = PETSC_FALSE;
4818: #if defined(PETSC_USE_SOCKET_VIEWER)
4819: PetscBool flg5 = PETSC_FALSE;
4820: #endif
4823: if (incall) return(0);
4824: incall = PETSC_TRUE;
4825: PetscObjectOptionsBegin((PetscObject)mat);
4826: PetscOptionsBool("-mat_view_info","Information on matrix size","MatView",flg1,&flg1,PETSC_NULL);
4827: PetscOptionsBool("-mat_view_info_detailed","Nonzeros in the matrix","MatView",flg2,&flg2,PETSC_NULL);
4828: PetscOptionsBool("-mat_view","Print matrix to stdout","MatView",flg3,&flg3,PETSC_NULL);
4829: PetscOptionsBool("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",flg4,&flg4,PETSC_NULL);
4830: #if defined(PETSC_USE_SOCKET_VIEWER)
4831: PetscOptionsBool("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",flg5,&flg5,PETSC_NULL);
4832: #endif
4833: PetscOptionsBool("-mat_view_binary","Save matrix to file in binary format","MatView",flg6,&flg6,PETSC_NULL);
4834: PetscOptionsBool("-mat_view_draw","Draw the matrix nonzero structure","MatView",flg7,&flg7,PETSC_NULL);
4835: PetscOptionsEnd();
4837: if (flg1) {
4838: PetscViewer viewer;
4840: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4841: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
4842: MatView(mat,viewer);
4843: PetscViewerPopFormat(viewer);
4844: }
4845: if (flg2) {
4846: PetscViewer viewer;
4848: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4849: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
4850: MatView(mat,viewer);
4851: PetscViewerPopFormat(viewer);
4852: }
4853: if (flg3) {
4854: PetscViewer viewer;
4856: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4857: MatView(mat,viewer);
4858: }
4859: if (flg4) {
4860: PetscViewer viewer;
4862: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4863: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
4864: MatView(mat,viewer);
4865: PetscViewerPopFormat(viewer);
4866: }
4867: #if defined(PETSC_USE_SOCKET_VIEWER)
4868: if (flg5) {
4869: MatView(mat,PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4870: PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4871: }
4872: #endif
4873: if (flg6) {
4874: MatView(mat,PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4875: PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4876: }
4877: if (flg7) {
4878: PetscOptionsGetBool(((PetscObject)mat)->prefix,"-mat_view_contour",&flg8,PETSC_NULL);
4879: if (flg8) {
4880: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm),PETSC_VIEWER_DRAW_CONTOUR);
4881: }
4882: MatView(mat,PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4883: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4884: if (flg8) {
4885: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4886: }
4887: }
4888: incall = PETSC_FALSE;
4889: return(0);
4890: }
4894: /*@
4895: MatAssemblyEnd - Completes assembling the matrix. This routine should
4896: be called after MatAssemblyBegin().
4898: Collective on Mat
4900: Input Parameters:
4901: + mat - the matrix
4902: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4904: Options Database Keys:
4905: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
4906: . -mat_view_info_detailed - Prints more detailed info
4907: . -mat_view - Prints matrix in ASCII format
4908: . -mat_view_matlab - Prints matrix in Matlab format
4909: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
4910: . -display <name> - Sets display name (default is host)
4911: . -draw_pause <sec> - Sets number of seconds to pause after display
4912: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (See the <a href="../../docs/manual.pdf">users manual</a>)
4913: . -viewer_socket_machine <machine>
4914: . -viewer_socket_port <port>
4915: . -mat_view_binary - save matrix to file in binary format
4916: - -viewer_binary_filename <name>
4918: Notes:
4919: MatSetValues() generally caches the values. The matrix is ready to
4920: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4921: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4922: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4923: using the matrix.
4925: Level: beginner
4927: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
4928: @*/
4929: PetscErrorCode MatAssemblyEnd(Mat mat,MatAssemblyType type)
4930: {
4931: PetscErrorCode ierr;
4932: static PetscInt inassm = 0;
4933: PetscBool flg = PETSC_FALSE;
4939: inassm++;
4940: MatAssemblyEnd_InUse++;
4941: if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
4942: PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
4943: if (mat->ops->assemblyend) {
4944: (*mat->ops->assemblyend)(mat,type);
4945: }
4946: PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
4947: } else {
4948: if (mat->ops->assemblyend) {
4949: (*mat->ops->assemblyend)(mat,type);
4950: }
4951: }
4953: /* Flush assembly is not a true assembly */
4954: if (type != MAT_FLUSH_ASSEMBLY) {
4955: mat->assembled = PETSC_TRUE; mat->num_ass++;
4956: }
4957: mat->insertmode = NOT_SET_VALUES;
4958: MatAssemblyEnd_InUse--;
4959: PetscObjectStateIncrease((PetscObject)mat);
4960: if (!mat->symmetric_eternal) {
4961: mat->symmetric_set =