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;
12: PetscClassId MAT_TRANSPOSECOLORING_CLASSID;
14: PetscLogEvent MAT_Mult, MAT_Mults, MAT_MultConstrained, MAT_MultAdd, MAT_MultTranspose;
15: PetscLogEvent MAT_MultTransposeConstrained, MAT_MultTransposeAdd, MAT_Solve, MAT_Solves, MAT_SolveAdd, MAT_SolveTranspose, MAT_MatSolve;
16: PetscLogEvent MAT_SolveTransposeAdd, MAT_SOR, MAT_ForwardSolve, MAT_BackwardSolve, MAT_LUFactor, MAT_LUFactorSymbolic;
17: PetscLogEvent MAT_LUFactorNumeric, MAT_CholeskyFactor, MAT_CholeskyFactorSymbolic, MAT_CholeskyFactorNumeric, MAT_ILUFactor;
18: PetscLogEvent MAT_ILUFactorSymbolic, MAT_ICCFactorSymbolic, MAT_Copy, MAT_Convert, MAT_Scale, MAT_AssemblyBegin;
19: PetscLogEvent MAT_AssemblyEnd, MAT_SetValues, MAT_GetValues, MAT_GetRow, MAT_GetRowIJ, MAT_GetSubMatrices, MAT_GetColoring, MAT_GetOrdering, MAT_GetRedundantMatrix, MAT_GetSeqNonzeroStructure;
20: PetscLogEvent MAT_IncreaseOverlap, MAT_Partitioning, MAT_ZeroEntries, MAT_Load, MAT_View, MAT_AXPY, MAT_FDColoringCreate;
21: PetscLogEvent MAT_FDColoringApply,MAT_Transpose,MAT_FDColoringFunction;
22: PetscLogEvent MAT_TransposeColoringCreate;
23: PetscLogEvent MAT_MatMult, MAT_MatMultSymbolic, MAT_MatMultNumeric;
24: PetscLogEvent MAT_PtAP, MAT_PtAPSymbolic, MAT_PtAPNumeric,MAT_RARt, MAT_RARtSymbolic, MAT_RARtNumeric;
25: PetscLogEvent MAT_MatTransposeMult, MAT_MatTransposeMultSymbolic, MAT_MatTransposeMultNumeric;
26: PetscLogEvent MAT_TransposeMatMult, MAT_TransposeMatMultSymbolic, MAT_TransposeMatMultNumeric;
27: PetscLogEvent MAT_MultHermitianTranspose,MAT_MultHermitianTransposeAdd;
28: PetscLogEvent MAT_Getsymtranspose, MAT_Getsymtransreduced, MAT_Transpose_SeqAIJ, MAT_GetBrowsOfAcols;
29: PetscLogEvent MAT_GetBrowsOfAocols, MAT_Getlocalmat, MAT_Getlocalmatcondensed, MAT_Seqstompi, MAT_Seqstompinum, MAT_Seqstompisym;
30: PetscLogEvent MAT_Applypapt, MAT_Applypapt_numeric, MAT_Applypapt_symbolic, MAT_GetSequentialNonzeroStructure;
31: PetscLogEvent MAT_GetMultiProcBlock;
32: PetscLogEvent MAT_CUSPCopyToGPU, MAT_SetValuesBatch, MAT_SetValuesBatchI, MAT_SetValuesBatchII, MAT_SetValuesBatchIII, MAT_SetValuesBatchIV;
33: PetscLogEvent MAT_Merge;
35: /* nasty global values for MatSetValue() */
36: PetscInt MatSetValue_Row = 0;
37: PetscInt MatSetValue_Column = 0;
38: PetscScalar MatSetValue_Value = 0.0;
40: const char *const MatFactorTypes[] = {"NONE","LU","CHOLESKY","ILU","ICC","ILUDT","MatFactorType","MAT_FACTOR_",0};
44: /*@C
45: MatFindNonzeroRows - Locate all rows that are not completely zero in the matrix
47: Input Parameter:
48: . A - the matrix
50: Output Parameter:
51: . keptrows - the rows that are not completely zero
53: Level: intermediate
55: @*/
56: PetscErrorCode MatFindNonzeroRows(Mat mat,IS *keptrows)
57: {
58: PetscErrorCode ierr;
62: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
63: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
64: if (!mat->ops->findnonzerorows) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Not coded for this matrix type");
65: (*mat->ops->findnonzerorows)(mat,keptrows);
66: return(0);
67: }
71: /*@
72: MatGetDiagonalBlock - Returns the part of the matrix associated with the on-process coupling
74: Not Collective
76: Input Parameters:
77: . A - the matrix
79: Output Parameters:
80: . a - the diagonal part (which is a SEQUENTIAL matrix)
82: Notes: see the manual page for MatCreateMPIAIJ() for more information on the "diagonal part" of the matrix.
84: Level: advanced
86: @*/
87: PetscErrorCode MatGetDiagonalBlock(Mat A,Mat *a)
88: {
89: PetscErrorCode ierr,(*f)(Mat,Mat*);
90: PetscMPIInt size;
96: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
97: if (A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
98: MPI_Comm_size(((PetscObject)A)->comm,&size);
99: PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);
100: if (f) {
101: (*f)(A,a);
102: return(0);
103: } else if (size == 1) {
104: *a = A;
105: } else {
106: const MatType mattype;
107: MatGetType(A,&mattype);
108: SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Matrix type %s does not support getting diagonal block",mattype);
109: }
110: return(0);
111: }
115: /*@
116: MatGetTrace - Gets the trace of a matrix. The sum of the diagonal entries.
118: Collective on Mat
120: Input Parameters:
121: . mat - the matrix
123: Output Parameter:
124: . trace - the sum of the diagonal entries
126: Level: advanced
128: @*/
129: PetscErrorCode MatGetTrace(Mat mat,PetscScalar *trace)
130: {
132: Vec diag;
135: MatGetVecs(mat,&diag,PETSC_NULL);
136: MatGetDiagonal(mat,diag);
137: VecSum(diag,trace);
138: VecDestroy(&diag);
139: return(0);
140: }
144: /*@
145: MatRealPart - Zeros out the imaginary part of the matrix
147: Logically Collective on Mat
149: Input Parameters:
150: . mat - the matrix
152: Level: advanced
155: .seealso: MatImaginaryPart()
156: @*/
157: PetscErrorCode MatRealPart(Mat mat)
158: {
164: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
165: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
166: if (!mat->ops->realpart) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
167: MatCheckPreallocated(mat,1);
168: (*mat->ops->realpart)(mat);
169: #if defined(PETSC_HAVE_CUSP)
170: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
171: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
172: }
173: #endif
174: return(0);
175: }
179: /*@C
180: MatGetGhosts - Get the global index of all ghost nodes defined by the sparse matrix
182: Collective on Mat
184: Input Parameter:
185: . mat - the matrix
187: Output Parameters:
188: + nghosts - number of ghosts (note for BAIJ matrices there is one ghost for each block)
189: - ghosts - the global indices of the ghost points
191: Notes: the nghosts and ghosts are suitable to pass into VecCreateGhost()
193: Level: advanced
195: @*/
196: PetscErrorCode MatGetGhosts(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[])
197: {
203: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
204: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
205: if (!mat->ops->getghosts) {
206: if (nghosts) *nghosts = 0;
207: if (ghosts) *ghosts = 0;
208: } else {
209: (*mat->ops->getghosts)(mat,nghosts,ghosts);
210: }
211: return(0);
212: }
217: /*@
218: MatImaginaryPart - Moves the imaginary part of the matrix to the real part and zeros the imaginary part
220: Logically Collective on Mat
222: Input Parameters:
223: . mat - the matrix
225: Level: advanced
228: .seealso: MatRealPart()
229: @*/
230: PetscErrorCode MatImaginaryPart(Mat mat)
231: {
237: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
238: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
239: if (!mat->ops->imaginarypart) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
240: MatCheckPreallocated(mat,1);
241: (*mat->ops->imaginarypart)(mat);
242: #if defined(PETSC_HAVE_CUSP)
243: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
244: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
245: }
246: #endif
247: return(0);
248: }
252: /*@
253: MatMissingDiagonal - Determine if sparse matrix is missing a diagonal entry (or block entry for BAIJ matrices)
255: Collective on Mat
257: Input Parameter:
258: . mat - the matrix
260: Output Parameters:
261: + missing - is any diagonal missing
262: - dd - first diagonal entry that is missing (optional)
264: Level: advanced
267: .seealso: MatRealPart()
268: @*/
269: PetscErrorCode MatMissingDiagonal(Mat mat,PetscBool *missing,PetscInt *dd)
270: {
276: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
277: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
278: if (!mat->ops->missingdiagonal) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
279: (*mat->ops->missingdiagonal)(mat,missing,dd);
280: return(0);
281: }
285: /*@C
286: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
287: for each row that you get to ensure that your application does
288: not bleed memory.
290: Not Collective
292: Input Parameters:
293: + mat - the matrix
294: - row - the row to get
296: Output Parameters:
297: + ncols - if not NULL, the number of nonzeros in the row
298: . cols - if not NULL, the column numbers
299: - vals - if not NULL, the values
301: Notes:
302: This routine is provided for people who need to have direct access
303: to the structure of a matrix. We hope that we provide enough
304: high-level matrix routines that few users will need it.
306: MatGetRow() always returns 0-based column indices, regardless of
307: whether the internal representation is 0-based (default) or 1-based.
309: For better efficiency, set cols and/or vals to PETSC_NULL if you do
310: not wish to extract these quantities.
312: The user can only examine the values extracted with MatGetRow();
313: the values cannot be altered. To change the matrix entries, one
314: must use MatSetValues().
316: You can only have one call to MatGetRow() outstanding for a particular
317: matrix at a time, per processor. MatGetRow() can only obtain rows
318: associated with the given processor, it cannot get rows from the
319: other processors; for that we suggest using MatGetSubMatrices(), then
320: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
321: is in the global number of rows.
323: Fortran Notes:
324: The calling sequence from Fortran is
325: .vb
326: MatGetRow(matrix,row,ncols,cols,values,ierr)
327: Mat matrix (input)
328: integer row (input)
329: integer ncols (output)
330: integer cols(maxcols) (output)
331: double precision (or double complex) values(maxcols) output
332: .ve
333: where maxcols >= maximum nonzeros in any row of the matrix.
336: Caution:
337: Do not try to change the contents of the output arrays (cols and vals).
338: In some cases, this may corrupt the matrix.
340: Level: advanced
342: Concepts: matrices^row access
344: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubMatrices(), MatGetDiagonal()
345: @*/
346: PetscErrorCode MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
347: {
349: PetscInt incols;
354: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
355: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
356: if (!mat->ops->getrow) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
357: MatCheckPreallocated(mat,1);
358: PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
359: (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
360: if (ncols) *ncols = incols;
361: PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
362: return(0);
363: }
367: /*@
368: MatConjugate - replaces the matrix values with their complex conjugates
370: Logically Collective on Mat
372: Input Parameters:
373: . mat - the matrix
375: Level: advanced
377: .seealso: VecConjugate()
378: @*/
379: PetscErrorCode MatConjugate(Mat mat)
380: {
385: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
386: 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");
387: (*mat->ops->conjugate)(mat);
388: #if defined(PETSC_HAVE_CUSP)
389: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
390: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
391: }
392: #endif
393: return(0);
394: }
398: /*@C
399: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
401: Not Collective
403: Input Parameters:
404: + mat - the matrix
405: . row - the row to get
406: . ncols, cols - the number of nonzeros and their columns
407: - vals - if nonzero the column values
409: Notes:
410: This routine should be called after you have finished examining the entries.
412: Fortran Notes:
413: The calling sequence from Fortran is
414: .vb
415: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
416: Mat matrix (input)
417: integer row (input)
418: integer ncols (output)
419: integer cols(maxcols) (output)
420: double precision (or double complex) values(maxcols) output
421: .ve
422: Where maxcols >= maximum nonzeros in any row of the matrix.
424: In Fortran MatRestoreRow() MUST be called after MatGetRow()
425: before another call to MatGetRow() can be made.
427: Level: advanced
429: .seealso: MatGetRow()
430: @*/
431: PetscErrorCode MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
432: {
438: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
439: if (!mat->ops->restorerow) return(0);
440: (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
441: return(0);
442: }
446: /*@
447: MatGetRowUpperTriangular - Sets a flag to enable calls to MatGetRow() for matrix in MATSBAIJ format.
448: You should call MatRestoreRowUpperTriangular() after calling MatGetRow/MatRestoreRow() to disable the flag.
450: Not Collective
452: Input Parameters:
453: + mat - the matrix
455: Notes:
456: 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.
458: Level: advanced
460: Concepts: matrices^row access
462: .seealso: MatRestoreRowRowUpperTriangular()
463: @*/
464: PetscErrorCode MatGetRowUpperTriangular(Mat mat)
465: {
471: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
472: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
473: if (!mat->ops->getrowuppertriangular) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
474: MatCheckPreallocated(mat,1);
475: (*mat->ops->getrowuppertriangular)(mat);
476: return(0);
477: }
481: /*@
482: MatRestoreRowUpperTriangular - Disable calls to MatGetRow() for matrix in MATSBAIJ format.
484: Not Collective
486: Input Parameters:
487: + mat - the matrix
489: Notes:
490: This routine should be called after you have finished MatGetRow/MatRestoreRow().
493: Level: advanced
495: .seealso: MatGetRowUpperTriangular()
496: @*/
497: PetscErrorCode MatRestoreRowUpperTriangular(Mat mat)
498: {
503: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
504: if (!mat->ops->restorerowuppertriangular) return(0);
505: (*mat->ops->restorerowuppertriangular)(mat);
506: return(0);
507: }
511: /*@C
512: MatSetOptionsPrefix - Sets the prefix used for searching for all
513: Mat options in the database.
515: Logically Collective on Mat
517: Input Parameter:
518: + A - the Mat context
519: - prefix - the prefix to prepend to all option names
521: Notes:
522: A hyphen (-) must NOT be given at the beginning of the prefix name.
523: The first character of all runtime options is AUTOMATICALLY the hyphen.
525: Level: advanced
527: .keywords: Mat, set, options, prefix, database
529: .seealso: MatSetFromOptions()
530: @*/
531: PetscErrorCode MatSetOptionsPrefix(Mat A,const char prefix[])
532: {
537: PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
538: return(0);
539: }
543: /*@C
544: MatAppendOptionsPrefix - Appends to the prefix used for searching for all
545: Mat options in the database.
547: Logically Collective on Mat
549: Input Parameters:
550: + A - the Mat context
551: - prefix - the prefix to prepend to all option names
553: Notes:
554: A hyphen (-) must NOT be given at the beginning of the prefix name.
555: The first character of all runtime options is AUTOMATICALLY the hyphen.
557: Level: advanced
559: .keywords: Mat, append, options, prefix, database
561: .seealso: MatGetOptionsPrefix()
562: @*/
563: PetscErrorCode MatAppendOptionsPrefix(Mat A,const char prefix[])
564: {
566:
569: PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
570: return(0);
571: }
575: /*@C
576: MatGetOptionsPrefix - Sets the prefix used for searching for all
577: Mat options in the database.
579: Not Collective
581: Input Parameter:
582: . A - the Mat context
584: Output Parameter:
585: . prefix - pointer to the prefix string used
587: Notes: On the fortran side, the user should pass in a string 'prefix' of
588: sufficient length to hold the prefix.
590: Level: advanced
592: .keywords: Mat, get, options, prefix, database
594: .seealso: MatAppendOptionsPrefix()
595: @*/
596: PetscErrorCode MatGetOptionsPrefix(Mat A,const char *prefix[])
597: {
602: PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
603: return(0);
604: }
608: /*@
609: MatSetUp - Sets up the internal matrix data structures for the later use.
611: Collective on Mat
613: Input Parameters:
614: . A - the Mat context
616: Notes:
617: If the user has not set preallocation for this matrix then a default preallocation that is likely to be inefficient is used.
619: If a suitable preallocation routine is used, this function does not need to be called.
621: See the Performance chapter of the PETSc users manual for how to preallocate matrices
623: Level: beginner
625: .keywords: Mat, setup
627: .seealso: MatCreate(), MatDestroy()
628: @*/
629: PetscErrorCode MatSetUp(Mat A)
630: {
631: PetscMPIInt size;
636: if (!((PetscObject)A)->type_name) {
637: MPI_Comm_size(((PetscObject)A)->comm, &size);
638: if (size == 1) {
639: MatSetType(A, MATSEQAIJ);
640: } else {
641: MatSetType(A, MATMPIAIJ);
642: }
643: }
644: if (!A->preallocated && A->ops->setup) {
645: PetscInfo(A,"Warning not preallocating matrix storage\n");
646: (*A->ops->setup)(A);
647: }
648: A->preallocated = PETSC_TRUE;
649: return(0);
650: }
655: /*@C
656: MatView - Visualizes a matrix object.
658: Collective on Mat
660: Input Parameters:
661: + mat - the matrix
662: - viewer - visualization context
664: Notes:
665: The available visualization contexts include
666: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
667: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
668: output where only the first processor opens
669: the file. All other processors send their
670: data to the first processor to print.
671: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
673: The user can open alternative visualization contexts with
674: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
675: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
676: specified file; corresponding input uses MatLoad()
677: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
678: an X window display
679: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
680: Currently only the sequential dense and AIJ
681: matrix types support the Socket viewer.
683: The user can call PetscViewerSetFormat() to specify the output
684: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
685: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
686: + PETSC_VIEWER_DEFAULT - default, prints matrix contents
687: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
688: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
689: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
690: format common among all matrix types
691: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
692: format (which is in many cases the same as the default)
693: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
694: size and structure (not the matrix entries)
695: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
696: the matrix structure
698: Options Database Keys:
699: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
700: . -mat_view_info_detailed - Prints more detailed info
701: . -mat_view - Prints matrix in ASCII format
702: . -mat_view_matlab - Prints matrix in Matlab format
703: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
704: . -display <name> - Sets display name (default is host)
705: . -draw_pause <sec> - Sets number of seconds to pause after display
706: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see the <a href="../../docs/manual.pdf">users manual</a> for details).
707: . -viewer_socket_machine <machine>
708: . -viewer_socket_port <port>
709: . -mat_view_binary - save matrix to file in binary format
710: - -viewer_binary_filename <name>
711: Level: beginner
713: Notes: see the manual page for MatLoad() for the exact format of the binary file when the binary
714: viewer is used.
716: See bin/matlab/PetscBinaryRead.m for a Matlab code that can read in the binary file when the binary
717: viewer is used.
719: Concepts: matrices^viewing
720: Concepts: matrices^plotting
721: Concepts: matrices^printing
723: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
724: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
725: @*/
726: PetscErrorCode MatView(Mat mat,PetscViewer viewer)
727: {
728: PetscErrorCode ierr;
729: PetscInt rows,cols;
730: PetscBool iascii;
731: PetscViewerFormat format;
736: if (!viewer) {
737: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
738: }
741: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
742: MatCheckPreallocated(mat,1);
744: PetscLogEventBegin(MAT_View,mat,viewer,0,0);
745: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
746: if (iascii) {
747: PetscViewerGetFormat(viewer,&format);
748: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
749: PetscObjectPrintClassNamePrefixType((PetscObject)mat,viewer,"Matrix Object");
750: PetscViewerASCIIPushTab(viewer);
751: MatGetSize(mat,&rows,&cols);
752: PetscViewerASCIIPrintf(viewer,"rows=%D, cols=%D\n",rows,cols);
753: if (mat->factortype) {
754: const MatSolverPackage solver;
755: MatFactorGetSolverPackage(mat,&solver);
756: PetscViewerASCIIPrintf(viewer,"package used to perform factorization: %s\n",solver);
757: }
758: if (mat->ops->getinfo) {
759: MatInfo info;
760: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
761: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",(PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
762: PetscViewerASCIIPrintf(viewer,"total number of mallocs used during MatSetValues calls =%D\n",(PetscInt)info.mallocs);
763: }
764: }
765: }
766: if (mat->ops->view) {
767: PetscViewerASCIIPushTab(viewer);
768: (*mat->ops->view)(mat,viewer);
769: PetscViewerASCIIPopTab(viewer);
770: } else if (!iascii) {
771: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
772: }
773: if (iascii) {
774: PetscViewerGetFormat(viewer,&format);
775: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
776: PetscViewerASCIIPopTab(viewer);
777: }
778: }
779: PetscLogEventEnd(MAT_View,mat,viewer,0,0);
780: return(0);
781: }
783: #if defined(PETSC_USE_DEBUG)
784: #include <../src/sys/totalview/tv_data_display.h>
785: PETSC_UNUSED static int TV_display_type(const struct _p_Mat *mat)
786: {
787: TV_add_row("Local rows", "int", &mat->rmap->n);
788: TV_add_row("Local columns", "int", &mat->cmap->n);
789: TV_add_row("Global rows", "int", &mat->rmap->N);
790: TV_add_row("Global columns", "int", &mat->cmap->N);
791: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)mat)->type_name);
792: return TV_format_OK;
793: }
794: #endif
798: /*@C
799: MatLoad - Loads a matrix that has been stored in binary format
800: with MatView(). The matrix format is determined from the options database.
801: Generates a parallel MPI matrix if the communicator has more than one
802: processor. The default matrix type is AIJ.
804: Collective on PetscViewer
806: Input Parameters:
807: + newmat - the newly loaded matrix, this needs to have been created with MatCreate()
808: or some related function before a call to MatLoad()
809: - viewer - binary file viewer, created with PetscViewerBinaryOpen()
811: Options Database Keys:
812: Used with block matrix formats (MATSEQBAIJ, ...) to specify
813: block size
814: . -matload_block_size <bs>
816: Level: beginner
818: Notes:
819: If the Mat type has not yet been given then MATAIJ is used, call MatSetFromOptions() on the
820: Mat before calling this routine if you wish to set it from the options database.
822: MatLoad() automatically loads into the options database any options
823: given in the file filename.info where filename is the name of the file
824: that was passed to the PetscViewerBinaryOpen(). The options in the info
825: file will be ignored if you use the -viewer_binary_skip_info option.
827: If the type or size of newmat is not set before a call to MatLoad, PETSc
828: sets the default matrix type AIJ and sets the local and global sizes.
829: If type and/or size is already set, then the same are used.
831: In parallel, each processor can load a subset of rows (or the
832: entire matrix). This routine is especially useful when a large
833: matrix is stored on disk and only part of it is desired on each
834: processor. For example, a parallel solver may access only some of
835: the rows from each processor. The algorithm used here reads
836: relatively small blocks of data rather than reading the entire
837: matrix and then subsetting it.
839: Notes for advanced users:
840: Most users should not need to know the details of the binary storage
841: format, since MatLoad() and MatView() completely hide these details.
842: But for anyone who's interested, the standard binary matrix storage
843: format is
845: $ int MAT_FILE_CLASSID
846: $ int number of rows
847: $ int number of columns
848: $ int total number of nonzeros
849: $ int *number nonzeros in each row
850: $ int *column indices of all nonzeros (starting index is zero)
851: $ PetscScalar *values of all nonzeros
853: PETSc automatically does the byte swapping for
854: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
855: linux, Windows and the paragon; thus if you write your own binary
856: read/write routines you have to swap the bytes; see PetscBinaryRead()
857: and PetscBinaryWrite() to see how this may be done.
859: .keywords: matrix, load, binary, input
861: .seealso: PetscViewerBinaryOpen(), MatView(), VecLoad()
863: @*/
864: PetscErrorCode MatLoad(Mat newmat,PetscViewer viewer)
865: {
867: PetscBool isbinary,flg;
872: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
873: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
875: if (!((PetscObject)newmat)->type_name) {
876: MatSetType(newmat,MATAIJ);
877: }
879: if (!newmat->ops->load) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatLoad is not supported for type");
880: PetscLogEventBegin(MAT_Load,viewer,0,0,0);
881: (*newmat->ops->load)(newmat,viewer);
882: PetscLogEventEnd(MAT_Load,viewer,0,0,0);
884: flg = PETSC_FALSE;
885: PetscOptionsGetBool(((PetscObject)newmat)->prefix,"-matload_symmetric",&flg,PETSC_NULL);
886: if (flg) {
887: MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);
888: MatSetOption(newmat,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
889: }
890: flg = PETSC_FALSE;
891: PetscOptionsGetBool(((PetscObject)newmat)->prefix,"-matload_spd",&flg,PETSC_NULL);
892: if (flg) {
893: MatSetOption(newmat,MAT_SPD,PETSC_TRUE);
894: }
895: return(0);
896: }
900: /*@
901: MatScaleSystem - Scale a vector solution and right hand side to
902: match the scaling of a scaled matrix.
903:
904: Collective on Mat
906: Input Parameter:
907: + mat - the matrix
908: . b - right hand side vector (or PETSC_NULL)
909: - x - solution vector (or PETSC_NULL)
912: Notes:
913: For AIJ, and BAIJ matrix formats, the matrices are not
914: internally scaled, so this does nothing.
916: The KSP methods automatically call this routine when required
917: (via PCPreSolve()) so it is rarely used directly.
919: Level: Developer
921: Concepts: matrices^scaling
923: .seealso: MatUseScaledForm(), MatUnScaleSystem()
924: @*/
925: PetscErrorCode MatScaleSystem(Mat mat,Vec b,Vec x)
926: {
932: MatCheckPreallocated(mat,1);
936: if (mat->ops->scalesystem) {
937: (*mat->ops->scalesystem)(mat,b,x);
938: }
939: return(0);
940: }
944: /*@
945: MatUnScaleSystem - Unscales a vector solution and right hand side to
946: match the original scaling of a scaled matrix.
947:
948: Collective on Mat
950: Input Parameter:
951: + mat - the matrix
952: . b - right hand side vector (or PETSC_NULL)
953: - x - solution vector (or PETSC_NULL)
956: Notes:
957: For AIJ and BAIJ matrix formats, the matrices are not
958: internally scaled, so this does nothing.
960: The KSP methods automatically call this routine when required
961: (via PCPreSolve()) so it is rarely used directly.
963: Level: Developer
965: .seealso: MatUseScaledForm(), MatScaleSystem()
966: @*/
967: PetscErrorCode MatUnScaleSystem(Mat mat,Vec b,Vec x)
968: {
974: MatCheckPreallocated(mat,1);
977: if (mat->ops->unscalesystem) {
978: (*mat->ops->unscalesystem)(mat,b,x);
979: }
980: return(0);
981: }
985: /*@
986: MatUseScaledForm - For matrix storage formats that scale the
987: matrix indicates matrix operations (MatMult() etc) are
988: applied using the scaled matrix.
989:
990: Logically Collective on Mat
992: Input Parameter:
993: + mat - the matrix
994: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
995: applying the original matrix
997: Notes:
998: For scaled matrix formats, applying the original, unscaled matrix
999: will be slightly more expensive
1001: Level: Developer
1003: .seealso: MatScaleSystem(), MatUnScaleSystem()
1004: @*/
1005: PetscErrorCode MatUseScaledForm(Mat mat,PetscBool scaled)
1006: {
1013: MatCheckPreallocated(mat,1);
1014: if (mat->ops->usescaledform) {
1015: (*mat->ops->usescaledform)(mat,scaled);
1016: }
1017: return(0);
1018: }
1022: /*@
1023: MatDestroy - Frees space taken by a matrix.
1025: Collective on Mat
1027: Input Parameter:
1028: . A - the matrix
1030: Level: beginner
1032: @*/
1033: PetscErrorCode MatDestroy(Mat *A)
1034: {
1038: if (!*A) return(0);
1040: if (--((PetscObject)(*A))->refct > 0) {*A = PETSC_NULL; return(0);}
1042: /* if memory was published with AMS then destroy it */
1043: PetscObjectDepublish(*A);
1045: if ((*A)->ops->destroy) {
1046: (*(*A)->ops->destroy)(*A);
1047: }
1049: MatNullSpaceDestroy(&(*A)->nullsp);
1050: PetscLayoutDestroy(&(*A)->rmap);
1051: PetscLayoutDestroy(&(*A)->cmap);
1052: PetscHeaderDestroy(A);
1053: return(0);
1054: }
1058: /*@
1059: MatSetValues - Inserts or adds a block of values into a matrix.
1060: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1061: MUST be called after all calls to MatSetValues() have been completed.
1063: Not Collective
1065: Input Parameters:
1066: + mat - the matrix
1067: . v - a logically two-dimensional array of values
1068: . m, idxm - the number of rows and their global indices
1069: . n, idxn - the number of columns and their global indices
1070: - addv - either ADD_VALUES or INSERT_VALUES, where
1071: ADD_VALUES adds values to any existing entries, and
1072: INSERT_VALUES replaces existing entries with new values
1074: Notes:
1075: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call MatXXXXSetPreallocation() or
1076: MatSetUp() before using this routine
1078: By default the values, v, are row-oriented. See MatSetOption() for other options.
1080: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
1081: options cannot be mixed without intervening calls to the assembly
1082: routines.
1084: MatSetValues() uses 0-based row and column numbers in Fortran
1085: as well as in C.
1087: Negative indices may be passed in idxm and idxn, these rows and columns are
1088: simply ignored. This allows easily inserting element stiffness matrices
1089: with homogeneous Dirchlet boundary conditions that you don't want represented
1090: in the matrix.
1092: Efficiency Alert:
1093: The routine MatSetValuesBlocked() may offer much better efficiency
1094: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1096: Level: beginner
1098: Concepts: matrices^putting entries in
1100: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1101: InsertMode, INSERT_VALUES, ADD_VALUES
1102: @*/
1103: PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1104: {
1110: if (!m || !n) return(0); /* no values to insert */
1114: MatCheckPreallocated(mat,1);
1115: if (mat->insertmode == NOT_SET_VALUES) {
1116: mat->insertmode = addv;
1117: }
1118: #if defined(PETSC_USE_DEBUG)
1119: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1120: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1121: if (!mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1122: #endif
1124: if (mat->assembled) {
1125: mat->was_assembled = PETSC_TRUE;
1126: }
1127: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1128: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
1129: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1130: #if defined(PETSC_HAVE_CUSP)
1131: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1132: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1133: }
1134: #endif
1135: return(0);
1136: }
1141: /*@
1142: MatSetValuesRowLocal - Inserts a row (block row for BAIJ matrices) of nonzero
1143: values into a matrix
1145: Not Collective
1147: Input Parameters:
1148: + mat - the matrix
1149: . row - the (block) row to set
1150: - v - a logically two-dimensional array of values
1152: Notes:
1153: By the values, v, are column-oriented (for the block version) and sorted
1155: All the nonzeros in the row must be provided
1157: The matrix must have previously had its column indices set
1159: The row must belong to this process
1161: Level: intermediate
1163: Concepts: matrices^putting entries in
1165: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1166: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues(), MatSetValuesRow(), MatSetLocalToGlobalMapping()
1167: @*/
1168: PetscErrorCode MatSetValuesRowLocal(Mat mat,PetscInt row,const PetscScalar v[])
1169: {
1176: MatSetValuesRow(mat, mat->rmap->mapping->indices[row],v);
1177: #if defined(PETSC_HAVE_CUSP)
1178: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1179: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1180: }
1181: #endif
1182: return(0);
1183: }
1187: /*@
1188: MatSetValuesRow - Inserts a row (block row for BAIJ matrices) of nonzero
1189: values into a matrix
1191: Not Collective
1193: Input Parameters:
1194: + mat - the matrix
1195: . row - the (block) row to set
1196: - v - a logically two-dimensional array of values
1198: Notes:
1199: The values, v, are column-oriented for the block version.
1201: All the nonzeros in the row must be provided
1203: THE MATRIX MUSAT HAVE PREVIOUSLY HAD ITS COLUMN INDICES SET. IT IS RARE THAT THIS ROUTINE IS USED, usually MatSetValues() is used.
1205: The row must belong to this process
1207: Level: advanced
1209: Concepts: matrices^putting entries in
1211: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1212: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
1213: @*/
1214: PetscErrorCode MatSetValuesRow(Mat mat,PetscInt row,const PetscScalar v[])
1215: {
1221: MatCheckPreallocated(mat,1);
1223: #if defined(PETSC_USE_DEBUG)
1224: if (mat->insertmode == ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add and insert values");
1225: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1226: #endif
1227: mat->insertmode = INSERT_VALUES;
1229: if (mat->assembled) {
1230: mat->was_assembled = PETSC_TRUE;
1231: mat->assembled = PETSC_FALSE;
1232: }
1233: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1234: if (!mat->ops->setvaluesrow) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1235: (*mat->ops->setvaluesrow)(mat,row,v);
1236: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1237: #if defined(PETSC_HAVE_CUSP)
1238: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1239: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1240: }
1241: #endif
1242: return(0);
1243: }
1247: /*@
1248: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
1249: Using structured grid indexing
1251: Not Collective
1253: Input Parameters:
1254: + mat - the matrix
1255: . m - number of rows being entered
1256: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
1257: . n - number of columns being entered
1258: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
1259: . v - a logically two-dimensional array of values
1260: - addv - either ADD_VALUES or INSERT_VALUES, where
1261: ADD_VALUES adds values to any existing entries, and
1262: INSERT_VALUES replaces existing entries with new values
1264: Notes:
1265: By default the values, v, are row-oriented. See MatSetOption() for other options.
1267: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
1268: options cannot be mixed without intervening calls to the assembly
1269: routines.
1271: The grid coordinates are across the entire grid, not just the local portion
1273: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
1274: as well as in C.
1276: For setting/accessing vector values via array coordinates you can use the DMDAVecGetArray() routine
1278: In order to use this routine you must either obtain the matrix with DMCreateMatrix()
1279: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1281: The columns and rows in the stencil passed in MUST be contained within the
1282: ghost region of the given process as set with DMDACreateXXX() or MatSetStencil(). For example,
1283: if you create a DMDA with an overlap of one grid level and on a particular process its first
1284: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1285: first i index you can use in your column and row indices in MatSetStencil() is 5.
1287: In Fortran idxm and idxn should be declared as
1288: $ MatStencil idxm(4,m),idxn(4,n)
1289: and the values inserted using
1290: $ idxm(MatStencil_i,1) = i
1291: $ idxm(MatStencil_j,1) = j
1292: $ idxm(MatStencil_k,1) = k
1293: $ idxm(MatStencil_c,1) = c
1294: etc
1295:
1296: For periodic boundary conditions use negative indices for values to the left (below 0; that are to be
1297: obtained by wrapping values from right edge). For values to the right of the last entry using that index plus one
1298: etc to obtain values that obtained by wrapping the values from the left edge. This does not work for anything but the
1299: DMDA_BOUNDARY_PERIODIC boundary type.
1301: 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
1302: a single value per point) you can skip filling those indices.
1304: Inspired by the structured grid interface to the HYPRE package
1305: (http://www.llnl.gov/CASC/hypre)
1307: Efficiency Alert:
1308: The routine MatSetValuesBlockedStencil() may offer much better efficiency
1309: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1311: Level: beginner
1313: Concepts: matrices^putting entries in
1315: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1316: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DMCreateMatrix(), DMDAVecGetArray(), MatStencil
1317: @*/
1318: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1319: {
1321: PetscInt buf[8192],*bufm=0,*bufn=0,*jdxm,*jdxn;
1322: PetscInt j,i,dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1323: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1326: if (!m || !n) return(0); /* no values to insert */
1333: if ((m+n) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1334: jdxm = buf; jdxn = buf+m;
1335: } else {
1336: PetscMalloc2(m,PetscInt,&bufm,n,PetscInt,&bufn);
1337: jdxm = bufm; jdxn = bufn;
1338: }
1339: for (i=0; i<m; i++) {
1340: for (j=0; j<3-sdim; j++) dxm++;
1341: tmp = *dxm++ - starts[0];
1342: for (j=0; j<dim-1; j++) {
1343: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1344: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1345: }
1346: if (mat->stencil.noc) dxm++;
1347: jdxm[i] = tmp;
1348: }
1349: for (i=0; i<n; i++) {
1350: for (j=0; j<3-sdim; j++) dxn++;
1351: tmp = *dxn++ - starts[0];
1352: for (j=0; j<dim-1; j++) {
1353: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1354: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1355: }
1356: if (mat->stencil.noc) dxn++;
1357: jdxn[i] = tmp;
1358: }
1359: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1360: PetscFree2(bufm,bufn);
1361: return(0);
1362: }
1366: /*@C
1367: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1368: Using structured grid indexing
1370: Not Collective
1372: Input Parameters:
1373: + mat - the matrix
1374: . m - number of rows being entered
1375: . idxm - grid coordinates for matrix rows being entered
1376: . n - number of columns being entered
1377: . idxn - grid coordinates for matrix columns being entered
1378: . v - a logically two-dimensional array of values
1379: - addv - either ADD_VALUES or INSERT_VALUES, where
1380: ADD_VALUES adds values to any existing entries, and
1381: INSERT_VALUES replaces existing entries with new values
1383: Notes:
1384: By default the values, v, are row-oriented and unsorted.
1385: See MatSetOption() for other options.
1387: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1388: options cannot be mixed without intervening calls to the assembly
1389: routines.
1391: The grid coordinates are across the entire grid, not just the local portion
1393: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1394: as well as in C.
1396: For setting/accessing vector values via array coordinates you can use the DMDAVecGetArray() routine
1398: In order to use this routine you must either obtain the matrix with DMCreateMatrix()
1399: or call MatSetBlockSize(), MatSetLocalToGlobalMapping() and MatSetStencil() first.
1401: The columns and rows in the stencil passed in MUST be contained within the
1402: ghost region of the given process as set with DMDACreateXXX() or MatSetStencil(). For example,
1403: if you create a DMDA with an overlap of one grid level and on a particular process its first
1404: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1405: first i index you can use in your column and row indices in MatSetStencil() is 5.
1407: In Fortran idxm and idxn should be declared as
1408: $ MatStencil idxm(4,m),idxn(4,n)
1409: and the values inserted using
1410: $ idxm(MatStencil_i,1) = i
1411: $ idxm(MatStencil_j,1) = j
1412: $ idxm(MatStencil_k,1) = k
1413: etc
1415: Negative indices may be passed in idxm and idxn, these rows and columns are
1416: simply ignored. This allows easily inserting element stiffness matrices
1417: with homogeneous Dirchlet boundary conditions that you don't want represented
1418: in the matrix.
1420: Inspired by the structured grid interface to the HYPRE package
1421: (http://www.llnl.gov/CASC/hypre)
1423: Level: beginner
1425: Concepts: matrices^putting entries in
1427: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1428: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DMCreateMatrix(), DMDAVecGetArray(), MatStencil,
1429: MatSetBlockSize(), MatSetLocalToGlobalMapping()
1430: @*/
1431: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1432: {
1434: PetscInt buf[8192],*bufm=0,*bufn=0,*jdxm,*jdxn;
1435: PetscInt j,i,dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1436: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1439: if (!m || !n) return(0); /* no values to insert */
1446: if ((m+n) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1447: jdxm = buf; jdxn = buf+m;
1448: } else {
1449: PetscMalloc2(m,PetscInt,&bufm,n,PetscInt,&bufn);
1450: jdxm = bufm; jdxn = bufn;
1451: }
1452: for (i=0; i<m; i++) {
1453: for (j=0; j<3-sdim; j++) dxm++;
1454: tmp = *dxm++ - starts[0];
1455: for (j=0; j<sdim-1; j++) {
1456: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1457: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1458: }
1459: dxm++;
1460: jdxm[i] = tmp;
1461: }
1462: for (i=0; i<n; i++) {
1463: for (j=0; j<3-sdim; j++) dxn++;
1464: tmp = *dxn++ - starts[0];
1465: for (j=0; j<sdim-1; j++) {
1466: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1467: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1468: }
1469: dxn++;
1470: jdxn[i] = tmp;
1471: }
1472: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1473: PetscFree2(bufm,bufn);
1474: #if defined(PETSC_HAVE_CUSP)
1475: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1476: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1477: }
1478: #endif
1479: return(0);
1480: }
1484: /*@
1485: MatSetStencil - Sets the grid information for setting values into a matrix via
1486: MatSetValuesStencil()
1488: Not Collective
1490: Input Parameters:
1491: + mat - the matrix
1492: . dim - dimension of the grid 1, 2, or 3
1493: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1494: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1495: - dof - number of degrees of freedom per node
1498: Inspired by the structured grid interface to the HYPRE package
1499: (www.llnl.gov/CASC/hyper)
1501: For matrices generated with DMCreateMatrix() this routine is automatically called and so not needed by the
1502: user.
1503:
1504: Level: beginner
1506: Concepts: matrices^putting entries in
1508: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1509: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1510: @*/
1511: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1512: {
1513: PetscInt i;
1520: mat->stencil.dim = dim + (dof > 1);
1521: for (i=0; i<dim; i++) {
1522: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1523: mat->stencil.starts[i] = starts[dim-i-1];
1524: }
1525: mat->stencil.dims[dim] = dof;
1526: mat->stencil.starts[dim] = 0;
1527: mat->stencil.noc = (PetscBool)(dof == 1);
1528: return(0);
1529: }
1533: /*@
1534: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1536: Not Collective
1538: Input Parameters:
1539: + mat - the matrix
1540: . v - a logically two-dimensional array of values
1541: . m, idxm - the number of block rows and their global block indices
1542: . n, idxn - the number of block columns and their global block indices
1543: - addv - either ADD_VALUES or INSERT_VALUES, where
1544: ADD_VALUES adds values to any existing entries, and
1545: INSERT_VALUES replaces existing entries with new values
1547: Notes:
1548: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call
1549: MatXXXXSetPreallocation() or MatSetUp() before using this routine.
1551: The m and n count the NUMBER of blocks in the row direction and column direction,
1552: NOT the total number of rows/columns; for example, if the block size is 2 and
1553: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1554: The values in idxm would be 1 2; that is the first index for each block divided by
1555: the block size.
1557: Note that you must call MatSetBlockSize() when constructing this matrix (after
1558: preallocating it).
1560: By default the values, v, are row-oriented, so the layout of
1561: v is the same as for MatSetValues(). See MatSetOption() for other options.
1563: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1564: options cannot be mixed without intervening calls to the assembly
1565: routines.
1567: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1568: as well as in C.
1570: Negative indices may be passed in idxm and idxn, these rows and columns are
1571: simply ignored. This allows easily inserting element stiffness matrices
1572: with homogeneous Dirchlet boundary conditions that you don't want represented
1573: in the matrix.
1575: Each time an entry is set within a sparse matrix via MatSetValues(),
1576: internal searching must be done to determine where to place the the
1577: data in the matrix storage space. By instead inserting blocks of
1578: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1579: reduced.
1581: Example:
1582: $ Suppose m=n=2 and block size(bs) = 2 The array is
1583: $
1584: $ 1 2 | 3 4
1585: $ 5 6 | 7 8
1586: $ - - - | - - -
1587: $ 9 10 | 11 12
1588: $ 13 14 | 15 16
1589: $
1590: $ v[] should be passed in like
1591: $ v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
1592: $
1593: $ If you are not using row oriented storage of v (that is you called MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE)) then
1594: $ v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
1596: Level: intermediate
1598: Concepts: matrices^putting entries in blocked
1600: .seealso: MatSetBlockSize(), MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1601: @*/
1602: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1603: {
1609: if (!m || !n) return(0); /* no values to insert */
1613: MatCheckPreallocated(mat,1);
1614: if (mat->insertmode == NOT_SET_VALUES) {
1615: mat->insertmode = addv;
1616: }
1617: #if defined(PETSC_USE_DEBUG)
1618: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1619: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1620: if (!mat->ops->setvaluesblocked && !mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1621: #endif
1623: if (mat->assembled) {
1624: mat->was_assembled = PETSC_TRUE;
1625: mat->assembled = PETSC_FALSE;
1626: }
1627: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1628: if (mat->ops->setvaluesblocked) {
1629: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1630: } else {
1631: PetscInt buf[8192],*bufr=0,*bufc=0,*iidxm,*iidxn;
1632: PetscInt i,j,bs=mat->rmap->bs;
1633: if ((m+n)*bs <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1634: iidxm = buf; iidxn = buf + m*bs;
1635: } else {
1636: PetscMalloc2(m*bs,PetscInt,&bufr,n*bs,PetscInt,&bufc);
1637: iidxm = bufr; iidxn = bufc;
1638: }
1639: for (i=0; i<m; i++)
1640: for (j=0; j<bs; j++)
1641: iidxm[i*bs+j] = bs*idxm[i] + j;
1642: for (i=0; i<n; i++)
1643: for (j=0; j<bs; j++)
1644: iidxn[i*bs+j] = bs*idxn[i] + j;
1645: MatSetValues(mat,m*bs,iidxm,n*bs,iidxn,v,addv);
1646: PetscFree2(bufr,bufc);
1647: }
1648: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1649: #if defined(PETSC_HAVE_CUSP)
1650: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1651: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1652: }
1653: #endif
1654: return(0);
1655: }
1659: /*@
1660: MatGetValues - Gets a block of values from a matrix.
1662: Not Collective; currently only returns a local block
1664: Input Parameters:
1665: + mat - the matrix
1666: . v - a logically two-dimensional array for storing the values
1667: . m, idxm - the number of rows and their global indices
1668: - n, idxn - the number of columns and their global indices
1670: Notes:
1671: The user must allocate space (m*n PetscScalars) for the values, v.
1672: The values, v, are then returned in a row-oriented format,
1673: analogous to that used by default in MatSetValues().
1675: MatGetValues() uses 0-based row and column numbers in
1676: Fortran as well as in C.
1678: MatGetValues() requires that the matrix has been assembled
1679: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1680: MatSetValues() and MatGetValues() CANNOT be made in succession
1681: without intermediate matrix assembly.
1683: Negative row or column indices will be ignored and those locations in v[] will be
1684: left unchanged.
1686: Level: advanced
1688: Concepts: matrices^accessing values
1690: .seealso: MatGetRow(), MatGetSubMatrices(), MatSetValues()
1691: @*/
1692: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1693: {
1699: if (!m || !n) return(0);
1703: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1704: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1705: if (!mat->ops->getvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1706: MatCheckPreallocated(mat,1);
1708: PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
1709: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1710: PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
1711: return(0);
1712: }
1716: /*@
1717: MatSetValuesBatch - Adds (ADD_VALUES) many blocks of values into a matrix at once. The blocks must all be square and
1718: the same size. Currently, this can only be called once and creates the given matrix.
1720: Not Collective
1722: Input Parameters:
1723: + mat - the matrix
1724: . nb - the number of blocks
1725: . bs - the number of rows (and columns) in each block
1726: . rows - a concatenation of the rows for each block
1727: - v - a concatenation of logically two-dimensional arrays of values
1729: Notes:
1730: In the future, we will extend this routine to handle rectangular blocks, and to allow multiple calls for a given matrix.
1732: Level: advanced
1734: Concepts: matrices^putting entries in
1736: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
1737: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
1738: @*/
1739: PetscErrorCode MatSetValuesBatch(Mat mat, PetscInt nb, PetscInt bs, PetscInt rows[], const PetscScalar v[])
1740: {
1748: #if defined(PETSC_USE_DEBUG)
1749: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1750: #endif
1752: PetscLogEventBegin(MAT_SetValuesBatch,mat,0,0,0);
1753: if (mat->ops->setvaluesbatch) {
1754: (*mat->ops->setvaluesbatch)(mat,nb,bs,rows,v);
1755: } else {
1756: PetscInt b;
1757: for(b = 0; b < nb; ++b) {
1758: MatSetValues(mat, bs, &rows[b*bs], bs, &rows[b*bs], &v[b*bs*bs], ADD_VALUES);
1759: }
1760: }
1761: PetscLogEventEnd(MAT_SetValuesBatch,mat,0,0,0);
1762: return(0);
1763: }
1767: /*@
1768: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1769: the routine MatSetValuesLocal() to allow users to insert matrix entries
1770: using a local (per-processor) numbering.
1772: Not Collective
1774: Input Parameters:
1775: + x - the matrix
1776: . rmapping - row mapping created with ISLocalToGlobalMappingCreate()
1777: or ISLocalToGlobalMappingCreateIS()
1778: - cmapping - column mapping
1780: Level: intermediate
1782: Concepts: matrices^local to global mapping
1783: Concepts: local to global mapping^for matrices
1785: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1786: @*/
1787: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1788: {
1795: MatCheckPreallocated(x,1);
1797: if (x->ops->setlocaltoglobalmapping) {
1798: (*x->ops->setlocaltoglobalmapping)(x,rmapping,cmapping);
1799: } else {
1800: PetscLayoutSetISLocalToGlobalMapping(x->rmap,rmapping);
1801: PetscLayoutSetISLocalToGlobalMapping(x->cmap,cmapping);
1802: }
1803: return(0);
1804: }
1808: /*@
1809: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
1810: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
1811: entries using a local (per-processor) numbering.
1813: Not Collective
1815: Input Parameters:
1816: + x - the matrix
1817: . rmapping - row mapping created with ISLocalToGlobalMappingCreate() or
1818: ISLocalToGlobalMappingCreateIS()
1819: - cmapping - column mapping
1821: Level: intermediate
1823: Concepts: matrices^local to global mapping blocked
1824: Concepts: local to global mapping^for matrices, blocked
1826: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
1827: MatSetValuesBlocked(), MatSetValuesLocal()
1828: @*/
1829: PetscErrorCode MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping)
1830: {
1837: MatCheckPreallocated(x,1);
1839: PetscLayoutSetISLocalToGlobalMappingBlock(x->rmap,rmapping);
1840: PetscLayoutSetISLocalToGlobalMappingBlock(x->cmap,cmapping);
1841: return(0);
1842: }
1846: /*@
1847: MatGetLocalToGlobalMapping - Gets the local-to-global numbering set by MatSetLocalToGlobalMapping()
1849: Not Collective
1851: Input Parameters:
1852: . A - the matrix
1854: Output Parameters:
1855: + rmapping - row mapping
1856: - cmapping - column mapping
1858: Level: advanced
1860: Concepts: matrices^local to global mapping
1861: Concepts: local to global mapping^for matrices
1863: .seealso: MatSetValuesLocal(), MatGetLocalToGlobalMappingBlock()
1864: @*/
1865: PetscErrorCode MatGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
1866: {
1872: if (rmapping) *rmapping = A->rmap->mapping;
1873: if (cmapping) *cmapping = A->cmap->mapping;
1874: return(0);
1875: }
1879: /*@
1880: MatGetLocalToGlobalMappingBlock - Gets the local-to-global numbering set by MatSetLocalToGlobalMappingBlock()
1882: Not Collective
1884: Input Parameters:
1885: . A - the matrix
1887: Output Parameters:
1888: + rmapping - row mapping
1889: - cmapping - column mapping
1891: Level: advanced
1893: Concepts: matrices^local to global mapping blocked
1894: Concepts: local to global mapping^for matrices, blocked
1896: .seealso: MatSetValuesBlockedLocal(), MatGetLocalToGlobalMapping()
1897: @*/
1898: PetscErrorCode MatGetLocalToGlobalMappingBlock(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping)
1899: {
1905: if (rmapping) *rmapping = A->rmap->bmapping;
1906: if (cmapping) *cmapping = A->cmap->bmapping;
1907: return(0);
1908: }
1912: /*@
1913: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
1914: using a local ordering of the nodes.
1916: Not Collective
1918: Input Parameters:
1919: + x - the matrix
1920: . nrow, irow - number of rows and their local indices
1921: . ncol, icol - number of columns and their local indices
1922: . y - a logically two-dimensional array of values
1923: - addv - either INSERT_VALUES or ADD_VALUES, where
1924: ADD_VALUES adds values to any existing entries, and
1925: INSERT_VALUES replaces existing entries with new values
1927: Notes:
1928: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call MatXXXXSetPreallocation() or
1929: MatSetUp() before using this routine
1931: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call MatSetLocalToGlobalMapping() before using this routine
1933: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
1934: options cannot be mixed without intervening calls to the assembly
1935: routines.
1937: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1938: MUST be called after all calls to MatSetValuesLocal() have been completed.
1940: Level: intermediate
1942: Concepts: matrices^putting entries in with local numbering
1944: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
1945: MatSetValueLocal()
1946: @*/
1947: PetscErrorCode MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1948: {
1954: MatCheckPreallocated(mat,1);
1955: if (!nrow || !ncol) return(0); /* no values to insert */
1959: if (mat->insertmode == NOT_SET_VALUES) {
1960: mat->insertmode = addv;
1961: }
1962: #if defined(PETSC_USE_DEBUG)
1963: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1964: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1965: if (!mat->ops->setvalueslocal && !mat->ops->setvalues) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1966: #endif
1968: if (mat->assembled) {
1969: mat->was_assembled = PETSC_TRUE;
1970: mat->assembled = PETSC_FALSE;
1971: }
1972: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1973: if (mat->ops->setvalueslocal) {
1974: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1975: } else {
1976: PetscInt buf[8192],*bufr=0,*bufc=0,*irowm,*icolm;
1977: if ((nrow+ncol) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
1978: irowm = buf; icolm = buf+nrow;
1979: } else {
1980: PetscMalloc2(nrow,PetscInt,&bufr,ncol,PetscInt,&bufc);
1981: irowm = bufr; icolm = bufc;
1982: }
1983: ISLocalToGlobalMappingApply(mat->rmap->mapping,nrow,irow,irowm);
1984: ISLocalToGlobalMappingApply(mat->cmap->mapping,ncol,icol,icolm);
1985: MatSetValues(mat,nrow,irowm,ncol,icolm,y,addv);
1986: PetscFree2(bufr,bufc);
1987: }
1988: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1989: #if defined(PETSC_HAVE_CUSP)
1990: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
1991: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
1992: }
1993: #endif
1994: return(0);
1995: }
1999: /*@
2000: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
2001: using a local ordering of the nodes a block at a time.
2003: Not Collective
2005: Input Parameters:
2006: + x - the matrix
2007: . nrow, irow - number of rows and their local indices
2008: . ncol, icol - number of columns and their local indices
2009: . y - a logically two-dimensional array of values
2010: - addv - either INSERT_VALUES or ADD_VALUES, where
2011: ADD_VALUES adds values to any existing entries, and
2012: INSERT_VALUES replaces existing entries with new values
2014: Notes:
2015: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call MatXXXXSetPreallocation() or
2016: MatSetUp() before using this routine
2018: If you create the matrix yourself (that is not with a call to DMCreateMatrix()) then you MUST call MatSetBlockSize() and MatSetLocalToGlobalMappingBlock()
2019: before using this routineBefore calling MatSetValuesLocal(), the user must first set the
2021: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
2022: options cannot be mixed without intervening calls to the assembly
2023: routines.
2025: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
2026: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
2028: Level: intermediate
2030: Concepts: matrices^putting blocked values in with local numbering
2032: .seealso: MatSetBlockSize(), MatSetLocalToGlobalMappingBlock(), MatAssemblyBegin(), MatAssemblyEnd(),
2033: MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
2034: @*/
2035: PetscErrorCode MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
2036: {
2042: MatCheckPreallocated(mat,1);
2043: if (!nrow || !ncol) return(0); /* no values to insert */
2047: if (mat->insertmode == NOT_SET_VALUES) {
2048: mat->insertmode = addv;
2049: }
2050: #if defined(PETSC_USE_DEBUG)
2051: else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
2052: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2053: 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);
2054: #endif
2056: if (mat->assembled) {
2057: mat->was_assembled = PETSC_TRUE;
2058: mat->assembled = PETSC_FALSE;
2059: }
2060: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
2061: if (mat->ops->setvaluesblockedlocal) {
2062: (*mat->ops->setvaluesblockedlocal)(mat,nrow,irow,ncol,icol,y,addv);
2063: } else {
2064: PetscInt buf[8192],*bufr=0,*bufc=0,*irowm,*icolm;
2065: if (mat->rmap->bmapping && mat->cmap->bmapping) {
2066: if ((nrow+ncol) <= (PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
2067: irowm = buf; icolm = buf + nrow;
2068: } else {
2069: PetscMalloc2(nrow,PetscInt,&bufr,ncol,PetscInt,&bufc);
2070: irowm = bufr; icolm = bufc;
2071: }
2072: ISLocalToGlobalMappingApply(mat->rmap->bmapping,nrow,irow,irowm);
2073: ISLocalToGlobalMappingApply(mat->cmap->bmapping,ncol,icol,icolm);
2074: MatSetValuesBlocked(mat,nrow,irowm,ncol,icolm,y,addv);
2075: PetscFree2(bufr,bufc);
2076: } else {
2077: PetscInt i,j,bs=mat->rmap->bs;
2078: if ((nrow+ncol)*bs <=(PetscInt)(sizeof(buf)/sizeof(PetscInt))) {
2079: irowm = buf; icolm = buf + nrow;
2080: } else {
2081: PetscMalloc2(nrow*bs,PetscInt,&bufr,ncol*bs,PetscInt,&bufc);
2082: irowm = bufr; icolm = bufc;
2083: }
2084: for (i=0; i<nrow; i++)
2085: for (j=0; j<bs; j++)
2086: irowm[i*bs+j] = irow[i]*bs+j;
2087: for (i=0; i<ncol; i++)
2088: for (j=0; j<bs; j++)
2089: icolm[i*bs+j] = icol[i]*bs+j;
2090: MatSetValuesLocal(mat,nrow*bs,irowm,ncol*bs,icolm,y,addv);
2091: PetscFree2(bufr,bufc);
2092: }
2093: }
2094: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
2095: #if defined(PETSC_HAVE_CUSP)
2096: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
2097: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
2098: }
2099: #endif
2100: return(0);
2101: }
2105: /*@
2106: MatMultDiagonalBlock - Computes the matrix-vector product, y = Dx. Where D is defined by the inode or block structure of the diagonal
2108: Collective on Mat and Vec
2110: Input Parameters:
2111: + mat - the matrix
2112: - x - the vector to be multiplied
2114: Output Parameters:
2115: . y - the result
2117: Notes:
2118: The vectors x and y cannot be the same. I.e., one cannot
2119: call MatMult(A,y,y).
2121: Level: developer
2123: Concepts: matrix-vector product
2125: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2126: @*/
2127: PetscErrorCode MatMultDiagonalBlock(Mat mat,Vec x,Vec y)
2128: {
2137: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2138: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2139: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2140: MatCheckPreallocated(mat,1);
2142: if (!mat->ops->multdiagonalblock) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
2143: (*mat->ops->multdiagonalblock)(mat,x,y);
2144: PetscObjectStateIncrease((PetscObject)y);
2145: return(0);
2146: }
2148: /* --------------------------------------------------------*/
2151: /*@
2152: MatMult - Computes the matrix-vector product, y = Ax.
2154: Neighbor-wise Collective on Mat and Vec
2156: Input Parameters:
2157: + mat - the matrix
2158: - x - the vector to be multiplied
2160: Output Parameters:
2161: . y - the result
2163: Notes:
2164: The vectors x and y cannot be the same. I.e., one cannot
2165: call MatMult(A,y,y).
2167: Level: beginner
2169: Concepts: matrix-vector product
2171: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2172: @*/
2173: PetscErrorCode MatMult(Mat mat,Vec x,Vec y)
2174: {
2182: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2183: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2184: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2185: #ifndef PETSC_HAVE_CONSTRAINTS
2186: 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);
2187: 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);
2188: 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);
2189: #endif
2190: VecValidValues(x,2,PETSC_TRUE);
2191: MatCheckPreallocated(mat,1);
2193: if (!mat->ops->mult) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
2194: PetscLogEventBegin(MAT_Mult,mat,x,y,0);
2195: (*mat->ops->mult)(mat,x,y);
2196: PetscLogEventEnd(MAT_Mult,mat,x,y,0);
2197: VecValidValues(y,3,PETSC_FALSE);
2198: return(0);
2199: }
2203: /*@
2204: MatMultTranspose - Computes matrix transpose times a vector.
2206: Neighbor-wise Collective on Mat and Vec
2208: Input Parameters:
2209: + mat - the matrix
2210: - x - the vector to be multilplied
2212: Output Parameters:
2213: . y - the result
2215: Notes:
2216: The vectors x and y cannot be the same. I.e., one cannot
2217: call MatMultTranspose(A,y,y).
2219: For complex numbers this does NOT compute the Hermitian (complex conjugate) transpose multiple,
2220: use MatMultHermitianTranspose()
2222: Level: beginner
2224: Concepts: matrix vector product^transpose
2226: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd(), MatMultHermitianTranspose(), MatTranspose()
2227: @*/
2228: PetscErrorCode MatMultTranspose(Mat mat,Vec x,Vec y)
2229: {
2238: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2239: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2240: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2241: #ifndef PETSC_HAVE_CONSTRAINTS
2242: 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);
2243: 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);
2244: #endif
2245: VecValidValues(x,2,PETSC_TRUE);
2246: MatCheckPreallocated(mat,1);
2248: if (!mat->ops->multtranspose) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
2249: PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
2250: (*mat->ops->multtranspose)(mat,x,y);
2251: PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
2252: PetscObjectStateIncrease((PetscObject)y);
2253: VecValidValues(y,3,PETSC_FALSE);
2254: return(0);
2255: }
2259: /*@
2260: MatMultHermitianTranspose - Computes matrix Hermitian transpose times a vector.
2262: Neighbor-wise Collective on Mat and Vec
2264: Input Parameters:
2265: + mat - the matrix
2266: - x - the vector to be multilplied
2268: Output Parameters:
2269: . y - the result
2271: Notes:
2272: The vectors x and y cannot be the same. I.e., one cannot
2273: call MatMultHermitianTranspose(A,y,y).
2275: Also called the conjugate transpose, complex conjugate transpose, or adjoint.
2277: For real numbers MatMultTranspose() and MatMultHermitianTranspose() are identical.
2279: Level: beginner
2281: Concepts: matrix vector product^transpose
2283: .seealso: MatMult(), MatMultAdd(), MatMultHermitianTransposeAdd(), MatMultTranspose()
2284: @*/
2285: PetscErrorCode MatMultHermitianTranspose(Mat mat,Vec x,Vec y)
2286: {
2295: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2296: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2297: if (x == y) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2298: #ifndef PETSC_HAVE_CONSTRAINTS
2299: 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);
2300: 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);
2301: #endif
2302: MatCheckPreallocated(mat,1);
2304: if (!mat->ops->multhermitiantranspose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2305: PetscLogEventBegin(MAT_MultHermitianTranspose,mat,x,y,0);
2306: (*mat->ops->multhermitiantranspose)(mat,x,y);
2307: PetscLogEventEnd(MAT_MultHermitianTranspose,mat,x,y,0);
2308: PetscObjectStateIncrease((PetscObject)y);
2309: return(0);
2310: }
2314: /*@
2315: MatMultAdd - Computes v3 = v2 + A * v1.
2317: Neighbor-wise Collective on Mat and Vec
2319: Input Parameters:
2320: + mat - the matrix
2321: - v1, v2 - the vectors
2323: Output Parameters:
2324: . v3 - the result
2326: Notes:
2327: The vectors v1 and v3 cannot be the same. I.e., one cannot
2328: call MatMultAdd(A,v1,v2,v1).
2330: Level: beginner
2332: Concepts: matrix vector product^addition
2334: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
2335: @*/
2336: PetscErrorCode MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2337: {
2347: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2348: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2349: 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);
2350: /* 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);
2351: 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); */
2352: 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);
2353: 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);
2354: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2355: MatCheckPreallocated(mat,1);
2357: if (!mat->ops->multadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"No MatMultAdd() for matrix type '%s'",((PetscObject)mat)->type_name);
2358: PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
2359: (*mat->ops->multadd)(mat,v1,v2,v3);
2360: PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
2361: PetscObjectStateIncrease((PetscObject)v3);
2362: return(0);
2363: }
2367: /*@
2368: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
2370: Neighbor-wise Collective on Mat and Vec
2372: Input Parameters:
2373: + mat - the matrix
2374: - v1, v2 - the vectors
2376: Output Parameters:
2377: . v3 - the result
2379: Notes:
2380: The vectors v1 and v3 cannot be the same. I.e., one cannot
2381: call MatMultTransposeAdd(A,v1,v2,v1).
2383: Level: beginner
2385: Concepts: matrix vector product^transpose and addition
2387: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
2388: @*/
2389: PetscErrorCode MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2390: {
2400: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2401: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2402: if (!mat->ops->multtransposeadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2403: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2404: 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);
2405: 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);
2406: 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);
2407: MatCheckPreallocated(mat,1);
2409: PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
2410: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
2411: PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
2412: PetscObjectStateIncrease((PetscObject)v3);
2413: return(0);
2414: }
2418: /*@
2419: MatMultHermitianTransposeAdd - Computes v3 = v2 + A^H * v1.
2421: Neighbor-wise Collective on Mat and Vec
2423: Input Parameters:
2424: + mat - the matrix
2425: - v1, v2 - the vectors
2427: Output Parameters:
2428: . v3 - the result
2430: Notes:
2431: The vectors v1 and v3 cannot be the same. I.e., one cannot
2432: call MatMultHermitianTransposeAdd(A,v1,v2,v1).
2434: Level: beginner
2436: Concepts: matrix vector product^transpose and addition
2438: .seealso: MatMultHermitianTranspose(), MatMultTranspose(), MatMultAdd(), MatMult()
2439: @*/
2440: PetscErrorCode MatMultHermitianTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
2441: {
2451: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2452: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2453: if (!mat->ops->multhermitiantransposeadd) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2454: if (v1 == v3) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
2455: 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);
2456: 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);
2457: 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);
2458: MatCheckPreallocated(mat,1);
2460: PetscLogEventBegin(MAT_MultHermitianTransposeAdd,mat,v1,v2,v3);
2461: (*mat->ops->multhermitiantransposeadd)(mat,v1,v2,v3);
2462: PetscLogEventEnd(MAT_MultHermitianTransposeAdd,mat,v1,v2,v3);
2463: PetscObjectStateIncrease((PetscObject)v3);
2464: return(0);
2465: }
2469: /*@
2470: MatMultConstrained - The inner multiplication routine for a
2471: constrained matrix P^T A P.
2473: Neighbor-wise Collective on Mat and Vec
2475: Input Parameters:
2476: + mat - the matrix
2477: - x - the vector to be multilplied
2479: Output Parameters:
2480: . y - the result
2482: Notes:
2483: The vectors x and y cannot be the same. I.e., one cannot
2484: call MatMult(A,y,y).
2486: Level: beginner
2488: .keywords: matrix, multiply, matrix-vector product, constraint
2489: .seealso: MatMult(), MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2490: @*/
2491: PetscErrorCode MatMultConstrained(Mat mat,Vec x,Vec y)
2492: {
2499: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2500: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2501: if (x == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2502: 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);
2503: 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);
2504: 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);
2506: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2507: (*mat->ops->multconstrained)(mat,x,y);
2508: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2509: PetscObjectStateIncrease((PetscObject)y);
2511: return(0);
2512: }
2516: /*@
2517: MatMultTransposeConstrained - The inner multiplication routine for a
2518: constrained matrix P^T A^T P.
2520: Neighbor-wise Collective on Mat and Vec
2522: Input Parameters:
2523: + mat - the matrix
2524: - x - the vector to be multilplied
2526: Output Parameters:
2527: . y - the result
2529: Notes:
2530: The vectors x and y cannot be the same. I.e., one cannot
2531: call MatMult(A,y,y).
2533: Level: beginner
2535: .keywords: matrix, multiply, matrix-vector product, constraint
2536: .seealso: MatMult(), MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
2537: @*/
2538: PetscErrorCode MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
2539: {
2546: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2547: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2548: if (x == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2549: 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);
2550: 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);
2552: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2553: (*mat->ops->multtransposeconstrained)(mat,x,y);
2554: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2555: PetscObjectStateIncrease((PetscObject)y);
2557: return(0);
2558: }
2562: /*@C
2563: MatGetFactorType - gets the type of factorization it is
2565: Note Collective
2566: as the flag
2568: Input Parameters:
2569: . mat - the matrix
2571: Output Parameters:
2572: . t - the type, one of MAT_FACTOR_NONE, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ILU, MAT_FACTOR_ICC,MAT_FACTOR_ILUDT
2574: Level: intermediate
2576: .seealso: MatFactorType, MatGetFactor()
2577: @*/
2578: PetscErrorCode MatGetFactorType(Mat mat,MatFactorType *t)
2579: {
2583: *t = mat->factortype;
2584: return(0);
2585: }
2587: /* ------------------------------------------------------------*/
2590: /*@C
2591: MatGetInfo - Returns information about matrix storage (number of
2592: nonzeros, memory, etc.).
2594: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used as the flag
2596: Input Parameters:
2597: . mat - the matrix
2599: Output Parameters:
2600: + flag - flag indicating the type of parameters to be returned
2601: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
2602: MAT_GLOBAL_SUM - sum over all processors)
2603: - info - matrix information context
2605: Notes:
2606: The MatInfo context contains a variety of matrix data, including
2607: number of nonzeros allocated and used, number of mallocs during
2608: matrix assembly, etc. Additional information for factored matrices
2609: is provided (such as the fill ratio, number of mallocs during
2610: factorization, etc.). Much of this info is printed to PETSC_STDOUT
2611: when using the runtime options
2612: $ -info -mat_view_info
2614: Example for C/C++ Users:
2615: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
2616: data within the MatInfo context. For example,
2617: .vb
2618: MatInfo info;
2619: Mat A;
2620: double mal, nz_a, nz_u;
2622: MatGetInfo(A,MAT_LOCAL,&info);
2623: mal = info.mallocs;
2624: nz_a = info.nz_allocated;
2625: .ve
2627: Example for Fortran Users:
2628: Fortran users should declare info as a double precision
2629: array of dimension MAT_INFO_SIZE, and then extract the parameters
2630: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
2631: a complete list of parameter names.
2632: .vb
2633: double precision info(MAT_INFO_SIZE)
2634: double precision mal, nz_a
2635: Mat A
2636: integer ierr
2638: call MatGetInfo(A,MAT_LOCAL,info,ierr)
2639: mal = info(MAT_INFO_MALLOCS)
2640: nz_a = info(MAT_INFO_NZ_ALLOCATED)
2641: .ve
2643: Level: intermediate
2645: Concepts: matrices^getting information on
2646:
2647: Developer Note: fortran interface is not autogenerated as the f90
2648: interface defintion cannot be generated correctly [due to MatInfo]
2650: .seealso: MatStashGetInfo()
2651:
2652: @*/
2653: PetscErrorCode MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
2654: {
2661: if (!mat->ops->getinfo) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2662: MatCheckPreallocated(mat,1);
2663: (*mat->ops->getinfo)(mat,flag,info);
2664: return(0);
2665: }
2667: /* ----------------------------------------------------------*/
2671: /*@C
2672: MatLUFactor - Performs in-place LU factorization of matrix.
2674: Collective on Mat
2676: Input Parameters:
2677: + mat - the matrix
2678: . row - row permutation
2679: . col - column permutation
2680: - info - options for factorization, includes
2681: $ fill - expected fill as ratio of original fill.
2682: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2683: $ Run with the option -info to determine an optimal value to use
2685: Notes:
2686: Most users should employ the simplified KSP interface for linear solvers
2687: instead of working directly with matrix algebra routines such as this.
2688: See, e.g., KSPCreate().
2690: This changes the state of the matrix to a factored matrix; it cannot be used
2691: for example with MatSetValues() unless one first calls MatSetUnfactored().
2693: Level: developer
2695: Concepts: matrices^LU factorization
2697: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
2698: MatGetOrdering(), MatSetUnfactored(), MatFactorInfo, MatGetFactor()
2700: Developer Note: fortran interface is not autogenerated as the f90
2701: interface defintion cannot be generated correctly [due to MatFactorInfo]
2703: @*/
2704: PetscErrorCode MatLUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2705: {
2707: MatFactorInfo tinfo;
2715: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2716: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2717: if (!mat->ops->lufactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2718: MatCheckPreallocated(mat,1);
2719: if (!info) {
2720: MatFactorInfoInitialize(&tinfo);
2721: info = &tinfo;
2722: }
2724: PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
2725: (*mat->ops->lufactor)(mat,row,col,info);
2726: PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
2727: PetscObjectStateIncrease((PetscObject)mat);
2728: return(0);
2729: }
2733: /*@C
2734: MatILUFactor - Performs in-place ILU factorization of matrix.
2736: Collective on Mat
2738: Input Parameters:
2739: + mat - the matrix
2740: . row - row permutation
2741: . col - column permutation
2742: - info - structure containing
2743: $ levels - number of levels of fill.
2744: $ expected fill - as ratio of original fill.
2745: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2746: missing diagonal entries)
2748: Notes:
2749: Probably really in-place only when level of fill is zero, otherwise allocates
2750: new space to store factored matrix and deletes previous memory.
2752: Most users should employ the simplified KSP interface for linear solvers
2753: instead of working directly with matrix algebra routines such as this.
2754: See, e.g., KSPCreate().
2756: Level: developer
2758: Concepts: matrices^ILU factorization
2760: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2762: Developer Note: fortran interface is not autogenerated as the f90
2763: interface defintion cannot be generated correctly [due to MatFactorInfo]
2765: @*/
2766: PetscErrorCode MatILUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2767: {
2776: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"matrix must be square");
2777: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2778: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2779: if (!mat->ops->ilufactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2780: MatCheckPreallocated(mat,1);
2782: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2783: (*mat->ops->ilufactor)(mat,row,col,info);
2784: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2785: PetscObjectStateIncrease((PetscObject)mat);
2786: return(0);
2787: }
2791: /*@C
2792: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
2793: Call this routine before calling MatLUFactorNumeric().
2795: Collective on Mat
2797: Input Parameters:
2798: + fact - the factor matrix obtained with MatGetFactor()
2799: . mat - the matrix
2800: . row, col - row and column permutations
2801: - info - options for factorization, includes
2802: $ fill - expected fill as ratio of original fill.
2803: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2804: $ Run with the option -info to determine an optimal value to use
2807: Notes:
2808: See the <a href="../../docs/manual.pdf">users manual</a> for additional information about
2809: choosing the fill factor for better efficiency.
2811: Most users should employ the simplified KSP interface for linear solvers
2812: instead of working directly with matrix algebra routines such as this.
2813: See, e.g., KSPCreate().
2815: Level: developer
2817: Concepts: matrices^LU symbolic factorization
2819: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2821: Developer Note: fortran interface is not autogenerated as the f90
2822: interface defintion cannot be generated correctly [due to MatFactorInfo]
2824: @*/
2825: PetscErrorCode MatLUFactorSymbolic(Mat fact,Mat mat,IS row,IS col,const MatFactorInfo *info)
2826: {
2836: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2837: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2838: if (!(fact)->ops->lufactorsymbolic) {
2839: const MatSolverPackage spackage;
2840: MatFactorGetSolverPackage(fact,&spackage);
2841: SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Matrix type %s symbolic LU using solver package %s",((PetscObject)mat)->type_name,spackage);
2842: }
2843: MatCheckPreallocated(mat,2);
2845: PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
2846: (fact->ops->lufactorsymbolic)(fact,mat,row,col,info);
2847: PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
2848: PetscObjectStateIncrease((PetscObject)fact);
2849: return(0);
2850: }
2854: /*@C
2855: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
2856: Call this routine after first calling MatLUFactorSymbolic().
2858: Collective on Mat
2860: Input Parameters:
2861: + fact - the factor matrix obtained with MatGetFactor()
2862: . mat - the matrix
2863: - info - options for factorization
2865: Notes:
2866: See MatLUFactor() for in-place factorization. See
2867: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
2869: Most users should employ the simplified KSP interface for linear solvers
2870: instead of working directly with matrix algebra routines such as this.
2871: See, e.g., KSPCreate().
2873: Level: developer
2875: Concepts: matrices^LU numeric factorization
2877: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
2879: Developer Note: fortran interface is not autogenerated as the f90
2880: interface defintion cannot be generated correctly [due to MatFactorInfo]
2882: @*/
2883: PetscErrorCode MatLUFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2884: {
2892: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2893: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2894: 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);
2895: }
2896: if (!(fact)->ops->lufactornumeric) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s numeric LU",((PetscObject)mat)->type_name);
2897: MatCheckPreallocated(mat,2);
2898: PetscLogEventBegin(MAT_LUFactorNumeric,mat,fact,0,0);
2899: (fact->ops->lufactornumeric)(fact,mat,info);
2900: PetscLogEventEnd(MAT_LUFactorNumeric,mat,fact,0,0);
2902: MatView_Private(fact);
2903: PetscObjectStateIncrease((PetscObject)fact);
2904: return(0);
2905: }
2909: /*@C
2910: MatCholeskyFactor - Performs in-place Cholesky factorization of a
2911: symmetric matrix.
2913: Collective on Mat
2915: Input Parameters:
2916: + mat - the matrix
2917: . perm - row and column permutations
2918: - f - expected fill as ratio of original fill
2920: Notes:
2921: See MatLUFactor() for the nonsymmetric case. See also
2922: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
2924: Most users should employ the simplified KSP interface for linear solvers
2925: instead of working directly with matrix algebra routines such as this.
2926: See, e.g., KSPCreate().
2928: Level: developer
2930: Concepts: matrices^Cholesky factorization
2932: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
2933: MatGetOrdering()
2935: Developer Note: fortran interface is not autogenerated as the f90
2936: interface defintion cannot be generated correctly [due to MatFactorInfo]
2938: @*/
2939: PetscErrorCode MatCholeskyFactor(Mat mat,IS perm,const MatFactorInfo *info)
2940: {
2948: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"Matrix must be square");
2949: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2950: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2951: if (!mat->ops->choleskyfactor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2952: MatCheckPreallocated(mat,1);
2954: PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
2955: (*mat->ops->choleskyfactor)(mat,perm,info);
2956: PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
2957: PetscObjectStateIncrease((PetscObject)mat);
2958: return(0);
2959: }
2963: /*@C
2964: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
2965: of a symmetric matrix.
2967: Collective on Mat
2969: Input Parameters:
2970: + fact - the factor matrix obtained with MatGetFactor()
2971: . mat - the matrix
2972: . perm - row and column permutations
2973: - info - options for factorization, includes
2974: $ fill - expected fill as ratio of original fill.
2975: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2976: $ Run with the option -info to determine an optimal value to use
2978: Notes:
2979: See MatLUFactorSymbolic() for the nonsymmetric case. See also
2980: MatCholeskyFactor() and MatCholeskyFactorNumeric().
2982: Most users should employ the simplified KSP interface for linear solvers
2983: instead of working directly with matrix algebra routines such as this.
2984: See, e.g., KSPCreate().
2986: Level: developer
2988: Concepts: matrices^Cholesky symbolic factorization
2990: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
2991: MatGetOrdering()
2993: Developer Note: fortran interface is not autogenerated as the f90
2994: interface defintion cannot be generated correctly [due to MatFactorInfo]
2996: @*/
2997: PetscErrorCode MatCholeskyFactorSymbolic(Mat fact,Mat mat,IS perm,const MatFactorInfo *info)
2998: {
3007: if (mat->rmap->N != mat->cmap->N) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"Matrix must be square");
3008: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3009: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3010: if (!(fact)->ops->choleskyfactorsymbolic) {
3011: const MatSolverPackage spackage;
3012: MatFactorGetSolverPackage(fact,&spackage);
3013: SETERRQ2(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s symbolic factor Cholesky using solver package %s",((PetscObject)mat)->type_name,spackage);
3014: }
3015: MatCheckPreallocated(mat,2);
3017: PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
3018: (fact->ops->choleskyfactorsymbolic)(fact,mat,perm,info);
3019: PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
3020: PetscObjectStateIncrease((PetscObject)fact);
3021: return(0);
3022: }
3026: /*@C
3027: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
3028: of a symmetric matrix. Call this routine after first calling
3029: MatCholeskyFactorSymbolic().
3031: Collective on Mat
3033: Input Parameters:
3034: + fact - the factor matrix obtained with MatGetFactor()
3035: . mat - the initial matrix
3036: . info - options for factorization
3037: - fact - the symbolic factor of mat
3040: Notes:
3041: Most users should employ the simplified KSP interface for linear solvers
3042: instead of working directly with matrix algebra routines such as this.
3043: See, e.g., KSPCreate().
3045: Level: developer
3047: Concepts: matrices^Cholesky numeric factorization
3049: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
3051: Developer Note: fortran interface is not autogenerated as the f90
3052: interface defintion cannot be generated correctly [due to MatFactorInfo]
3054: @*/
3055: PetscErrorCode MatCholeskyFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
3056: {
3064: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3065: if (!(fact)->ops->choleskyfactornumeric) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s numeric factor Cholesky",((PetscObject)mat)->type_name);
3066: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
3067: 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);
3068: }
3069: MatCheckPreallocated(mat,2);
3071: PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,fact,0,0);
3072: (fact->ops->choleskyfactornumeric)(fact,mat,info);
3073: PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,fact,0,0);
3075: MatView_Private(fact);
3076: PetscObjectStateIncrease((PetscObject)fact);
3077: return(0);
3078: }
3080: /* ----------------------------------------------------------------*/
3083: /*@
3084: MatSolve - Solves A x = b, given a factored matrix.
3086: Neighbor-wise Collective on Mat and Vec
3088: Input Parameters:
3089: + mat - the factored matrix
3090: - b - the right-hand-side vector
3092: Output Parameter:
3093: . x - the result vector
3095: Notes:
3096: The vectors b and x cannot be the same. I.e., one cannot
3097: call MatSolve(A,x,x).
3099: Notes:
3100: Most users should employ the simplified KSP interface for linear solvers
3101: instead of working directly with matrix algebra routines such as this.
3102: See, e.g., KSPCreate().
3104: Level: developer
3106: Concepts: matrices^triangular solves
3108: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
3109: @*/
3110: PetscErrorCode MatSolve(Mat mat,Vec b,Vec x)
3111: {
3121: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3122: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3123: 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);
3124: 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);
3125: 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);
3126: if (!mat->rmap->N && !mat->cmap->N) return(0);
3127: if (!mat->ops->solve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3128: MatCheckPreallocated(mat,1);
3130: PetscLogEventBegin(MAT_Solve,mat,b,x,0);
3131: (*mat->ops->solve)(mat,b,x);
3132: PetscLogEventEnd(MAT_Solve,mat,b,x,0);
3133: PetscObjectStateIncrease((PetscObject)x);
3134: return(0);
3135: }
3139: PetscErrorCode MatMatSolve_Basic(Mat A,Mat B,Mat X)
3140: {
3142: Vec b,x;
3143: PetscInt m,N,i;
3144: PetscScalar *bb,*xx;
3145: PetscBool flg;
3148: PetscTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
3149: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
3150: PetscTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
3151: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
3153: MatGetArray(B,&bb);
3154: MatGetArray(X,&xx);
3155: MatGetLocalSize(B,&m,PETSC_NULL); /* number local rows */
3156: MatGetSize(B,PETSC_NULL,&N); /* total columns in dense matrix */
3157: MatGetVecs(A,&x,&b);
3158: for (i=0; i<N; i++) {
3159: VecPlaceArray(b,bb + i*m);
3160: VecPlaceArray(x,xx + i*m);
3161: MatSolve(A,b,x);
3162: VecResetArray(x);
3163: VecResetArray(b);
3164: }
3165: VecDestroy(&b);
3166: VecDestroy(&x);
3167: MatRestoreArray(B,&bb);
3168: MatRestoreArray(X,&xx);
3169: return(0);
3170: }
3174: /*@
3175: MatMatSolve - Solves A X = B, given a factored matrix.
3177: Neighbor-wise Collective on Mat
3179: Input Parameters:
3180: + mat - the factored matrix
3181: - B - the right-hand-side matrix (dense matrix)
3183: Output Parameter:
3184: . X - the result matrix (dense matrix)
3186: Notes:
3187: The matrices b and x cannot be the same. I.e., one cannot
3188: call MatMatSolve(A,x,x).
3190: Notes:
3191: Most users should usually employ the simplified KSP interface for linear solvers
3192: instead of working directly with matrix algebra routines such as this.
3193: See, e.g., KSPCreate(). However KSP can only solve for one vector (column of X)
3194: at a time.
3196: When using SuperLU_Dist as a parallel solver PETSc will use the SuperLU_Dist functionality to solve multiple right hand sides simultaneously. For MUMPS
3197: it calls a separate solve for each right hand side since MUMPS does not yet support distributed right hand sides.
3199: Since the resulting matrix X must always be dense we do not support sparse representation of the matrix B.
3201: Level: developer
3203: Concepts: matrices^triangular solves
3205: .seealso: MatMatSolveAdd(), MatMatSolveTranspose(), MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()
3206: @*/
3207: PetscErrorCode MatMatSolve(Mat A,Mat B,Mat X)
3208: {
3218: if (X == B) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_IDN,"X and B must be different matrices");
3219: if (!A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3220: 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);
3221: 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);
3222: 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);
3223: if (!A->rmap->N && !A->cmap->N) return(0);
3224: MatCheckPreallocated(A,1);
3226: PetscLogEventBegin(MAT_MatSolve,A,B,X,0);
3227: if (!A->ops->matsolve) {
3228: PetscInfo1(A,"Mat type %s using basic MatMatSolve",((PetscObject)A)->type_name);
3229: MatMatSolve_Basic(A,B,X);
3230: } else {
3231: (*A->ops->matsolve)(A,B,X);
3232: }
3233: PetscLogEventEnd(MAT_MatSolve,A,B,X,0);
3234: PetscObjectStateIncrease((PetscObject)X);
3235: return(0);
3236: }
3241: /*@
3242: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU, or
3243: U^T*D^(1/2) x = b, given a factored symmetric matrix, A = U^T*D*U,
3245: Neighbor-wise Collective on Mat and Vec
3247: Input Parameters:
3248: + mat - the factored matrix
3249: - b - the right-hand-side vector
3251: Output Parameter:
3252: . x - the result vector
3254: Notes:
3255: MatSolve() should be used for most applications, as it performs
3256: a forward solve followed by a backward solve.
3258: The vectors b and x cannot be the same, i.e., one cannot
3259: call MatForwardSolve(A,x,x).
3261: For matrix in seqsbaij format with block size larger than 1,
3262: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
3263: MatForwardSolve() solves U^T*D y = b, and
3264: MatBackwardSolve() solves U x = y.
3265: Thus they do not provide a symmetric preconditioner.
3267: Most users should employ the simplified KSP interface for linear solvers
3268: instead of working directly with matrix algebra routines such as this.
3269: See, e.g., KSPCreate().
3271: Level: developer
3273: Concepts: matrices^forward solves
3275: .seealso: MatSolve(), MatBackwardSolve()
3276: @*/
3277: PetscErrorCode MatForwardSolve(Mat mat,Vec b,Vec x)
3278: {
3288: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3289: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3290: if (!mat->ops->forwardsolve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3291: 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);
3292: 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);
3293: 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);
3294: MatCheckPreallocated(mat,1);
3295: PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
3296: (*mat->ops->forwardsolve)(mat,b,x);
3297: PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
3298: PetscObjectStateIncrease((PetscObject)x);
3299: return(0);
3300: }
3304: /*@
3305: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
3306: D^(1/2) U x = b, given a factored symmetric matrix, A = U^T*D*U,
3308: Neighbor-wise Collective on Mat and Vec
3310: Input Parameters:
3311: + mat - the factored matrix
3312: - b - the right-hand-side vector
3314: Output Parameter:
3315: . x - the result vector
3317: Notes:
3318: MatSolve() should be used for most applications, as it performs
3319: a forward solve followed by a backward solve.
3321: The vectors b and x cannot be the same. I.e., one cannot
3322: call MatBackwardSolve(A,x,x).
3324: For matrix in seqsbaij format with block size larger than 1,
3325: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
3326: MatForwardSolve() solves U^T*D y = b, and
3327: MatBackwardSolve() solves U x = y.
3328: Thus they do not provide a symmetric preconditioner.
3330: Most users should employ the simplified KSP interface for linear solvers
3331: instead of working directly with matrix algebra routines such as this.
3332: See, e.g., KSPCreate().
3334: Level: developer
3336: Concepts: matrices^backward solves
3338: .seealso: MatSolve(), MatForwardSolve()
3339: @*/
3340: PetscErrorCode MatBackwardSolve(Mat mat,Vec b,Vec x)
3341: {
3351: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3352: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3353: if (!mat->ops->backwardsolve) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3354: 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);
3355: 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);
3356: 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);
3357: MatCheckPreallocated(mat,1);
3359: PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
3360: (*mat->ops->backwardsolve)(mat,b,x);
3361: PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
3362: PetscObjectStateIncrease((PetscObject)x);
3363: return(0);
3364: }
3368: /*@
3369: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
3371: Neighbor-wise Collective on Mat and Vec
3373: Input Parameters:
3374: + mat - the factored matrix
3375: . b - the right-hand-side vector
3376: - y - the vector to be added to
3378: Output Parameter:
3379: . x - the result vector
3381: Notes:
3382: The vectors b and x cannot be the same. I.e., one cannot
3383: call MatSolveAdd(A,x,y,x).
3385: Most users should employ the simplified KSP interface for linear solvers
3386: instead of working directly with matrix algebra routines such as this.
3387: See, e.g., KSPCreate().
3389: Level: developer
3391: Concepts: matrices^triangular solves
3393: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
3394: @*/
3395: PetscErrorCode MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
3396: {
3397: PetscScalar one = 1.0;
3398: Vec tmp;
3410: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3411: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3412: 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);
3413: 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);
3414: 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);
3415: 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);
3416: 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);
3417: MatCheckPreallocated(mat,1);
3419: PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
3420: if (mat->ops->solveadd) {
3421: (*mat->ops->solveadd)(mat,b,y,x);
3422: } else {
3423: /* do the solve then the add manually */
3424: if (x != y) {
3425: MatSolve(mat,b,x);
3426: VecAXPY(x,one,y);
3427: } else {
3428: VecDuplicate(x,&tmp);
3429: PetscLogObjectParent(mat,tmp);
3430: VecCopy(x,tmp);
3431: MatSolve(mat,b,x);
3432: VecAXPY(x,one,tmp);
3433: VecDestroy(&tmp);
3434: }
3435: }
3436: PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
3437: PetscObjectStateIncrease((PetscObject)x);
3438: return(0);
3439: }
3443: /*@
3444: MatSolveTranspose - Solves A' x = b, given a factored matrix.
3446: Neighbor-wise Collective on Mat and Vec
3448: Input Parameters:
3449: + mat - the factored matrix
3450: - b - the right-hand-side vector
3452: Output Parameter:
3453: . x - the result vector
3455: Notes:
3456: The vectors b and x cannot be the same. I.e., one cannot
3457: call MatSolveTranspose(A,x,x).
3459: Most users should employ the simplified KSP interface for linear solvers
3460: instead of working directly with matrix algebra routines such as this.
3461: See, e.g., KSPCreate().
3463: Level: developer
3465: Concepts: matrices^triangular solves
3467: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
3468: @*/
3469: PetscErrorCode MatSolveTranspose(Mat mat,Vec b,Vec x)
3470: {
3480: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3481: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3482: if (!mat->ops->solvetranspose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Matrix type %s",((PetscObject)mat)->type_name);
3483: 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);
3484: 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);
3485: MatCheckPreallocated(mat,1);
3486: PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
3487: (*mat->ops->solvetranspose)(mat,b,x);
3488: PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
3489: PetscObjectStateIncrease((PetscObject)x);
3490: return(0);
3491: }
3495: /*@
3496: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
3497: factored matrix.
3499: Neighbor-wise Collective on Mat and Vec
3501: Input Parameters:
3502: + mat - the factored matrix
3503: . b - the right-hand-side vector
3504: - y - the vector to be added to
3506: Output Parameter:
3507: . x - the result vector
3509: Notes:
3510: The vectors b and x cannot be the same. I.e., one cannot
3511: call MatSolveTransposeAdd(A,x,y,x).
3513: Most users should employ the simplified KSP interface for linear solvers
3514: instead of working directly with matrix algebra routines such as this.
3515: See, e.g., KSPCreate().
3517: Level: developer
3519: Concepts: matrices^triangular solves
3521: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
3522: @*/
3523: PetscErrorCode MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
3524: {
3525: PetscScalar one = 1.0;
3527: Vec tmp;
3538: if (x == b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3539: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3540: 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);
3541: 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);
3542: 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);
3543: 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);
3544: MatCheckPreallocated(mat,1);
3546: PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
3547: if (mat->ops->solvetransposeadd) {
3548: (*mat->ops->solvetransposeadd)(mat,b,y,x);
3549: } else {
3550: /* do the solve then the add manually */
3551: if (x != y) {
3552: MatSolveTranspose(mat,b,x);
3553: VecAXPY(x,one,y);
3554: } else {
3555: VecDuplicate(x,&tmp);
3556: PetscLogObjectParent(mat,tmp);
3557: VecCopy(x,tmp);
3558: MatSolveTranspose(mat,b,x);
3559: VecAXPY(x,one,tmp);
3560: VecDestroy(&tmp);
3561: }
3562: }
3563: PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
3564: PetscObjectStateIncrease((PetscObject)x);
3565: return(0);
3566: }
3567: /* ----------------------------------------------------------------*/
3571: /*@
3572: MatSOR - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3574: Neighbor-wise Collective on Mat and Vec
3576: Input Parameters:
3577: + mat - the matrix
3578: . b - the right hand side
3579: . omega - the relaxation factor
3580: . flag - flag indicating the type of SOR (see below)
3581: . shift - diagonal shift
3582: . its - the number of iterations
3583: - lits - the number of local iterations
3585: Output Parameters:
3586: . x - the solution (can contain an initial guess, use option SOR_ZERO_INITIAL_GUESS to indicate no guess)
3588: SOR Flags:
3589: . SOR_FORWARD_SWEEP - forward SOR
3590: . SOR_BACKWARD_SWEEP - backward SOR
3591: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
3592: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
3593: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
3594: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
3595: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
3596: upper/lower triangular part of matrix to
3597: vector (with omega)
3598: . SOR_ZERO_INITIAL_GUESS - zero initial guess
3600: Notes:
3601: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
3602: SOR_LOCAL_SYMMETRIC_SWEEP perform separate independent smoothings
3603: on each processor.
3605: Application programmers will not generally use MatSOR() directly,
3606: but instead will employ the KSP/PC interface.
3608: Notes: for BAIJ, SBAIJ, and AIJ matrices with Inodes this does a block SOR smoothing, otherwise it does a pointwise smoothing
3610: Notes for Advanced Users:
3611: The flags are implemented as bitwise inclusive or operations.
3612: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
3613: to specify a zero initial guess for SSOR.
3615: Most users should employ the simplified KSP interface for linear solvers
3616: instead of working directly with matrix algebra routines such as this.
3617: See, e.g., KSPCreate().
3619: Vectors x and b CANNOT be the same
3621: Developer Note: We should add block SOR support for AIJ matrices with block size set to great than one and no inodes
3623: Level: developer
3625: Concepts: matrices^relaxation
3626: Concepts: matrices^SOR
3627: Concepts: matrices^Gauss-Seidel
3629: @*/
3630: PetscErrorCode MatSOR(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3631: {
3641: if (!mat->ops->sor) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3642: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3643: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3644: 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);
3645: 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);
3646: 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);
3647: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
3648: if (lits <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires local its %D positive",lits);
3649: if (b == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"b and x vector cannot be the same");
3651: MatCheckPreallocated(mat,1);
3652: PetscLogEventBegin(MAT_SOR,mat,b,x,0);
3653: ierr =(*mat->ops->sor)(mat,b,omega,flag,shift,its,lits,x);
3654: PetscLogEventEnd(MAT_SOR,mat,b,x,0);
3655: PetscObjectStateIncrease((PetscObject)x);
3656: return(0);
3657: }
3661: /*
3662: Default matrix copy routine.
3663: */
3664: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
3665: {
3666: PetscErrorCode ierr;
3667: PetscInt i,rstart = 0,rend = 0,nz;
3668: const PetscInt *cwork;
3669: const PetscScalar *vwork;
3672: if (B->assembled) {
3673: MatZeroEntries(B);
3674: }
3675: MatGetOwnershipRange(A,&rstart,&rend);
3676: for (i=rstart; i<rend; i++) {
3677: MatGetRow(A,i,&nz,&cwork,&vwork);
3678: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
3679: MatRestoreRow(A,i,&nz,&cwork,&vwork);
3680: }
3681: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3682: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3683: PetscObjectStateIncrease((PetscObject)B);
3684: return(0);
3685: }
3689: /*@
3690: MatCopy - Copys a matrix to another matrix.
3692: Collective on Mat
3694: Input Parameters:
3695: + A - the matrix
3696: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
3698: Output Parameter:
3699: . B - where the copy is put
3701: Notes:
3702: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
3703: same nonzero pattern or the routine will crash.
3705: MatCopy() copies the matrix entries of a matrix to another existing
3706: matrix (after first zeroing the second matrix). A related routine is
3707: MatConvert(), which first creates a new matrix and then copies the data.
3709: Level: intermediate
3710:
3711: Concepts: matrices^copying
3713: .seealso: MatConvert(), MatDuplicate()
3715: @*/
3716: PetscErrorCode MatCopy(Mat A,Mat B,MatStructure str)
3717: {
3719: PetscInt i;
3727: MatCheckPreallocated(B,2);
3728: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3729: if (A->factortype) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3730: 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);
3731: MatCheckPreallocated(A,1);
3733: PetscLogEventBegin(MAT_Copy,A,B,0,0);
3734: if (A->ops->copy) {
3735: (*A->ops->copy)(A,B,str);
3736: } else { /* generic conversion */
3737: MatCopy_Basic(A,B,str);
3738: }
3740: B->stencil.dim = A->stencil.dim;
3741: B->stencil.noc = A->stencil.noc;
3742: for (i=0; i<=A->stencil.dim; i++) {
3743: B->stencil.dims[i] = A->stencil.dims[i];
3744: B->stencil.starts[i] = A->stencil.starts[i];
3745: }
3747: PetscLogEventEnd(MAT_Copy,A,B,0,0);
3748: PetscObjectStateIncrease((PetscObject)B);
3749: return(0);
3750: }
3754: /*@C
3755: MatConvert - Converts a matrix to another matrix, either of the same
3756: or different type.
3758: Collective on Mat
3760: Input Parameters:
3761: + mat - the matrix
3762: . newtype - new matrix type. Use MATSAME to create a new matrix of the
3763: same type as the original matrix.
3764: - reuse - denotes if the destination matrix is to be created or reused. Currently
3765: MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
3766: MAT_INITIAL_MATRIX.
3768: Output Parameter:
3769: . M - pointer to place new matrix
3771: Notes:
3772: MatConvert() first creates a new matrix and then copies the data from
3773: the first matrix. A related routine is MatCopy(), which copies the matrix
3774: entries of one matrix to another already existing matrix context.
3776: Cannot be used to convert a sequential matrix to parallel or parallel to sequential,
3777: the MPI communicator of the generated matrix is always the same as the communicator
3778: of the input matrix.
3780: Level: intermediate
3782: Concepts: matrices^converting between storage formats
3784: .seealso: MatCopy(), MatDuplicate()
3785: @*/
3786: PetscErrorCode MatConvert(Mat mat, const MatType newtype,MatReuse reuse,Mat *M)
3787: {
3788: PetscErrorCode ierr;
3789: PetscBool sametype,issame,flg;
3790: char convname[256],mtype[256];
3791: Mat B;
3797: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3798: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3799: MatCheckPreallocated(mat,1);
3800: MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
3802: PetscOptionsGetString(((PetscObject)mat)->prefix,"-matconvert_type",mtype,256,&flg);
3803: if (flg) {
3804: newtype = mtype;
3805: }
3806: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
3807: PetscStrcmp(newtype,"same",&issame);
3808: if ((reuse == MAT_REUSE_MATRIX) && (mat != *M)) {
3809: SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for in-place conversion currently");
3810: }
3812: if ((reuse == MAT_REUSE_MATRIX) && (issame || sametype)) return(0);
3813:
3814: if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
3815: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
3816: } else {
3817: PetscErrorCode (*conv)(Mat, const MatType,MatReuse,Mat*)=PETSC_NULL;
3818: const char *prefix[3] = {"seq","mpi",""};
3819: PetscInt i;
3820: /*
3821: Order of precedence:
3822: 1) See if a specialized converter is known to the current matrix.
3823: 2) See if a specialized converter is known to the desired matrix class.
3824: 3) See if a good general converter is registered for the desired class
3825: (as of 6/27/03 only MATMPIADJ falls into this category).
3826: 4) See if a good general converter is known for the current matrix.
3827: 5) Use a really basic converter.
3828: */
3829:
3830: /* 1) See if a specialized converter is known to the current matrix and the desired class */
3831: for (i=0; i<3; i++) {
3832: PetscStrcpy(convname,"MatConvert_");
3833: PetscStrcat(convname,((PetscObject)mat)->type_name);
3834: PetscStrcat(convname,"_");
3835: PetscStrcat(convname,prefix[i]);
3836: PetscStrcat(convname,issame?((PetscObject)mat)->type_name:newtype);
3837: PetscStrcat(convname,"_C");
3838: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3839: if (conv) goto foundconv;
3840: }
3842: /* 2) See if a specialized converter is known to the desired matrix class. */
3843: MatCreate(((PetscObject)mat)->comm,&B);
3844: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N);
3845: MatSetType(B,newtype);
3846: for (i=0; i<3; i++) {
3847: PetscStrcpy(convname,"MatConvert_");
3848: PetscStrcat(convname,((PetscObject)mat)->type_name);
3849: PetscStrcat(convname,"_");
3850: PetscStrcat(convname,prefix[i]);
3851: PetscStrcat(convname,newtype);
3852: PetscStrcat(convname,"_C");
3853: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
3854: if (conv) {
3855: MatDestroy(&B);
3856: goto foundconv;
3857: }
3858: }
3860: /* 3) See if a good general converter is registered for the desired class */
3861: conv = B->ops->convertfrom;
3862: MatDestroy(&B);
3863: if (conv) goto foundconv;
3865: /* 4) See if a good general converter is known for the current matrix */
3866: if (mat->ops->convert) {
3867: conv = mat->ops->convert;
3868: }
3869: if (conv) goto foundconv;
3871: /* 5) Use a really basic converter. */
3872: conv = MatConvert_Basic;
3874: foundconv:
3875: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3876: (*conv)(mat,newtype,reuse,M);
3877: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3878: }
3879: PetscObjectStateIncrease((PetscObject)*M);
3881: /* Copy Mat options */
3882: if (mat->symmetric){MatSetOption(*M,MAT_SYMMETRIC,PETSC_TRUE);}
3883: if (mat->hermitian){MatSetOption(*M,MAT_HERMITIAN,PETSC_TRUE);}
3884: return(0);
3885: }
3889: /*@C
3890: MatFactorGetSolverPackage - Returns name of the package providing the factorization routines
3892: Not Collective
3894: Input Parameter:
3895: . mat - the matrix, must be a factored matrix
3897: Output Parameter:
3898: . type - the string name of the package (do not free this string)
3900: Notes:
3901: In Fortran you pass in a empty string and the package name will be copied into it.
3902: (Make sure the string is long enough)
3904: Level: intermediate
3906: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable(), MatGetFactor()
3907: @*/
3908: PetscErrorCode MatFactorGetSolverPackage(Mat mat, const MatSolverPackage *type)
3909: {
3910: PetscErrorCode ierr;
3911: PetscErrorCode (*conv)(Mat,const MatSolverPackage*);
3916: if (!mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Only for factored matrix");
3917: PetscObjectQueryFunction((PetscObject)mat,"MatFactorGetSolverPackage_C",(void (**)(void))&conv);
3918: if (!conv) {
3919: *type = MATSOLVERPETSC;
3920: } else {
3921: (*conv)(mat,type);
3922: }
3923: return(0);
3924: }
3928: /*@C
3929: MatGetFactor - Returns a matrix suitable to calls to MatXXFactorSymbolic()
3931: Collective on Mat
3933: Input Parameters:
3934: + mat - the matrix
3935: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3936: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3938: Output Parameters:
3939: . f - the factor matrix used with MatXXFactorSymbolic() calls
3941: Notes:
3942: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3943: such as pastix, superlu, mumps, spooles etc.
3945: PETSc must have been ./configure to use the external solver, using the option --download-package
3947: Level: intermediate
3949: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable()
3950: @*/
3951: PetscErrorCode MatGetFactor(Mat mat, const MatSolverPackage type,MatFactorType ftype,Mat *f)
3952: {
3953: PetscErrorCode ierr,(*conv)(Mat,MatFactorType,Mat*);
3954: char convname[256];
3960: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3961: MatCheckPreallocated(mat,1);
3963: PetscStrcpy(convname,"MatGetFactor_");
3964: PetscStrcat(convname,type);
3965: PetscStrcat(convname,"_C");
3966: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3967: if (!conv) {
3968: PetscBool flag;
3969: MPI_Comm comm;
3971: PetscObjectGetComm((PetscObject)mat,&comm);
3972: PetscStrcasecmp(MATSOLVERPETSC,type,&flag);
3973: if (flag) {
3974: SETERRQ2(comm,PETSC_ERR_SUP,"Matrix format %s does not have a built-in PETSc %s",((PetscObject)mat)->type_name,MatFactorTypes[ftype]);
3975: } else {
3976: 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);
3977: }
3978: }
3979: (*conv)(mat,ftype,f);
3980: return(0);
3981: }
3985: /*@C
3986: MatGetFactorAvailable - Returns a a flag if matrix supports particular package and factor type
3988: Not Collective
3990: Input Parameters:
3991: + mat - the matrix
3992: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3993: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3995: Output Parameter:
3996: . flg - PETSC_TRUE if the factorization is available
3998: Notes:
3999: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
4000: such as pastix, superlu, mumps, spooles etc.
4002: PETSc must have been ./configure to use the external solver, using the option --download-package
4004: Level: intermediate
4006: .seealso: MatCopy(), MatDuplicate(), MatGetFactor()
4007: @*/
4008: PetscErrorCode MatGetFactorAvailable(Mat mat, const MatSolverPackage type,MatFactorType ftype,PetscBool *flg)
4009: {
4010: PetscErrorCode ierr;
4011: char convname[256];
4012: PetscErrorCode (*conv)(Mat,MatFactorType,PetscBool *);
4018: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4019: MatCheckPreallocated(mat,1);
4021: PetscStrcpy(convname,"MatGetFactorAvailable_");
4022: PetscStrcat(convname,type);
4023: PetscStrcat(convname,"_C");
4024: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
4025: if (!conv) {
4026: *flg = PETSC_FALSE;
4027: } else {
4028: (*conv)(mat,ftype,flg);
4029: }
4030: return(0);
4031: }
4036: /*@
4037: MatDuplicate - Duplicates a matrix including the non-zero structure.
4039: Collective on Mat
4041: Input Parameters:
4042: + mat - the matrix
4043: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy the numerical values in the matrix
4044: MAT_SHARE_NONZERO_PATTERN to share the nonzero patterns with the previous matrix and not copy them.
4046: Output Parameter:
4047: . M - pointer to place new matrix
4049: Level: intermediate
4051: Concepts: matrices^duplicating
4053: Notes: You cannot change the nonzero pattern for the parent or child matrix if you use MAT_SHARE_NONZERO_PATTERN.
4055: .seealso: MatCopy(), MatConvert()
4056: @*/
4057: PetscErrorCode MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
4058: {
4060: Mat B;
4061: PetscInt i;
4067: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4068: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4069: MatCheckPreallocated(mat,1);
4071: *M = 0;
4072: if (!mat->ops->duplicate) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Not written for this matrix type");
4073: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
4074: (*mat->ops->duplicate)(mat,op,M);
4075: B = *M;
4076:
4077: B->stencil.dim = mat->stencil.dim;
4078: B->stencil.noc = mat->stencil.noc;
4079: for (i=0; i<=mat->stencil.dim; i++) {
4080: B->stencil.dims[i] = mat->stencil.dims[i];
4081: B->stencil.starts[i] = mat->stencil.starts[i];
4082: }
4084: B->nooffproczerorows = mat->nooffproczerorows;
4085: B->nooffprocentries = mat->nooffprocentries;
4086: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
4087: PetscObjectStateIncrease((PetscObject)B);
4088: return(0);
4089: }
4093: /*@
4094: MatGetDiagonal - Gets the diagonal of a matrix.
4096: Logically Collective on Mat and Vec
4098: Input Parameters:
4099: + mat - the matrix
4100: - v - the vector for storing the diagonal
4102: Output Parameter:
4103: . v - the diagonal of the matrix
4105: Level: intermediate
4107: Note:
4108: Currently only correct in parallel for square matrices.
4110: Concepts: matrices^accessing diagonals
4112: .seealso: MatGetRow(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs()
4113: @*/
4114: PetscErrorCode MatGetDiagonal(Mat mat,Vec v)
4115: {
4122: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4123: if (!mat->ops->getdiagonal) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4124: MatCheckPreallocated(mat,1);
4126: (*mat->ops->getdiagonal)(mat,v);
4127: PetscObjectStateIncrease((PetscObject)v);
4128: return(0);
4129: }
4133: /*@
4134: MatGetRowMin - Gets the minimum value (of the real part) of each
4135: row of the matrix
4137: Logically Collective on Mat and Vec
4139: Input Parameters:
4140: . mat - the matrix
4142: Output Parameter:
4143: + v - the vector for storing the maximums
4144: - idx - the indices of the column found for each row (optional)
4146: Level: intermediate
4148: Notes: The result of this call are the same as if one converted the matrix to dense format
4149: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
4151: This code is only implemented for a couple of matrix formats.
4153: Concepts: matrices^getting row maximums
4155: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(),
4156: MatGetRowMax()
4157: @*/
4158: PetscErrorCode MatGetRowMin(Mat mat,Vec v,PetscInt idx[])
4159: {
4166: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4167: if (!mat->ops->getrowmax) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4168: MatCheckPreallocated(mat,1);
4170: (*mat->ops->getrowmin)(mat,v,idx);
4171: PetscObjectStateIncrease((PetscObject)v);
4172: return(0);
4173: }
4177: /*@
4178: MatGetRowMinAbs - Gets the minimum value (in absolute value) of each
4179: row of the matrix
4181: Logically Collective on Mat and Vec
4183: Input Parameters:
4184: . mat - the matrix
4186: Output Parameter:
4187: + v - the vector for storing the minimums
4188: - idx - the indices of the column found for each row (or PETSC_NULL if not needed)
4190: Level: intermediate
4192: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
4193: row is 0 (the first column).
4195: This code is only implemented for a couple of matrix formats.
4197: Concepts: matrices^getting row maximums
4199: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMaxAbs(), MatGetRowMin()
4200: @*/
4201: PetscErrorCode MatGetRowMinAbs(Mat mat,Vec v,PetscInt idx[])
4202: {
4209: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4210: if (!mat->ops->getrowminabs) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4211: MatCheckPreallocated(mat,1);
4212: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
4214: (*mat->ops->getrowminabs)(mat,v,idx);
4215: PetscObjectStateIncrease((PetscObject)v);
4216: return(0);
4217: }
4221: /*@
4222: MatGetRowMax - Gets the maximum value (of the real part) of each
4223: row of the matrix
4225: Logically Collective on Mat and Vec
4227: Input Parameters:
4228: . mat - the matrix
4230: Output Parameter:
4231: + v - the vector for storing the maximums
4232: - idx - the indices of the column found for each row (optional)
4234: Level: intermediate
4236: Notes: The result of this call are the same as if one converted the matrix to dense format
4237: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
4239: This code is only implemented for a couple of matrix formats.
4241: Concepts: matrices^getting row maximums
4243: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(), MatGetRowMin()
4244: @*/
4245: PetscErrorCode MatGetRowMax(Mat mat,Vec v,PetscInt idx[])
4246: {
4253: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4254: if (!mat->ops->getrowmax) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4255: MatCheckPreallocated(mat,1);
4257: (*mat->ops->getrowmax)(mat,v,idx);
4258: PetscObjectStateIncrease((PetscObject)v);
4259: return(0);
4260: }
4264: /*@
4265: MatGetRowMaxAbs - Gets the maximum value (in absolute value) of each
4266: row of the matrix
4268: Logically Collective on Mat and Vec
4270: Input Parameters:
4271: . mat - the matrix
4273: Output Parameter:
4274: + v - the vector for storing the maximums
4275: - idx - the indices of the column found for each row (or PETSC_NULL if not needed)
4277: Level: intermediate
4279: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
4280: row is 0 (the first column).
4282: This code is only implemented for a couple of matrix formats.
4284: Concepts: matrices^getting row maximums
4286: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
4287: @*/
4288: PetscErrorCode MatGetRowMaxAbs(Mat mat,Vec v,PetscInt idx[])
4289: {
4296: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4297: if (!mat->ops->getrowmaxabs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4298: MatCheckPreallocated(mat,1);
4299: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
4301: (*mat->ops->getrowmaxabs)(mat,v,idx);
4302: PetscObjectStateIncrease((PetscObject)v);
4303: return(0);
4304: }
4308: /*@
4309: MatGetRowSum - Gets the sum of each row of the matrix
4311: Logically Collective on Mat and Vec
4313: Input Parameters:
4314: . mat - the matrix
4316: Output Parameter:
4317: . v - the vector for storing the sum of rows
4319: Level: intermediate
4321: Notes: This code is slow since it is not currently specialized for different formats
4323: Concepts: matrices^getting row sums
4325: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
4326: @*/
4327: PetscErrorCode MatGetRowSum(Mat mat, Vec v)
4328: {
4329: PetscInt start = 0, end = 0, row;
4330: PetscScalar *array;
4337: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4338: MatCheckPreallocated(mat,1);
4339: MatGetOwnershipRange(mat, &start, &end);
4340: VecGetArray(v, &array);
4341: for(row = start; row < end; ++row) {
4342: PetscInt ncols, col;
4343: const PetscInt *cols;
4344: const PetscScalar *vals;
4346: array[row - start] = 0.0;
4347: MatGetRow(mat, row, &ncols, &cols, &vals);
4348: for(col = 0; col < ncols; col++) {
4349: array[row - start] += vals[col];
4350: }
4351: MatRestoreRow(mat, row, &ncols, &cols, &vals);
4352: }
4353: VecRestoreArray(v, &array);
4354: PetscObjectStateIncrease((PetscObject) v);
4355: return(0);
4356: }
4360: /*@
4361: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
4363: Collective on Mat
4365: Input Parameter:
4366: + mat - the matrix to transpose
4367: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4369: Output Parameters:
4370: . B - the transpose
4372: Notes:
4373: If you pass in &mat for B the transpose will be done in place, for example MatTranspose(mat,MAT_REUSE_MATRIX,&mat);
4375: Consider using MatCreateTranspose() instead if you only need a matrix that behaves like the transpose, but don't need the storage to be changed.
4377: Level: intermediate
4379: Concepts: matrices^transposing
4381: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose(), MatReuse
4382: @*/
4383: PetscErrorCode MatTranspose(Mat mat,MatReuse reuse,Mat *B)
4384: {
4390: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4391: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4392: if (!mat->ops->transpose) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4393: MatCheckPreallocated(mat,1);
4395: PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
4396: (*mat->ops->transpose)(mat,reuse,B);
4397: PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
4398: if (B) {PetscObjectStateIncrease((PetscObject)*B);}
4399: return(0);
4400: }
4404: /*@
4405: MatIsTranspose - Test whether a matrix is another one's transpose,
4406: or its own, in which case it tests symmetry.
4408: Collective on Mat
4410: Input Parameter:
4411: + A - the matrix to test
4412: - B - the matrix to test against, this can equal the first parameter
4414: Output Parameters:
4415: . flg - the result
4417: Notes:
4418: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
4419: has a running time of the order of the number of nonzeros; the parallel
4420: test involves parallel copies of the block-offdiagonal parts of the matrix.
4422: Level: intermediate
4424: Concepts: matrices^transposing, matrix^symmetry
4426: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
4427: @*/
4428: PetscErrorCode MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscBool *flg)
4429: {
4430: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscBool *),(*g)(Mat,Mat,PetscReal,PetscBool *);
4436: PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
4437: PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
4438: *flg = PETSC_FALSE;
4439: if (f && g) {
4440: if (f == g) {
4441: (*f)(A,B,tol,flg);
4442: } else {
4443: SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
4444: }
4445: } else {
4446: const MatType mattype;
4447: if (!f) {MatGetType(A,&mattype);}
4448: else {MatGetType(B,&mattype);}
4449: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix of type <%s> does not support checking for transpose",mattype);
4450: }
4451: return(0);
4452: }
4456: /*@
4457: MatHermitianTranspose - Computes an in-place or out-of-place transpose of a matrix in complex conjugate.
4459: Collective on Mat
4461: Input Parameter:
4462: + mat - the matrix to transpose and complex conjugate
4463: - reuse - store the transpose matrix in the provided B
4465: Output Parameters:
4466: . B - the Hermitian
4468: Notes:
4469: If you pass in &mat for B the Hermitian will be done in place
4471: Level: intermediate
4473: Concepts: matrices^transposing, complex conjugatex
4475: .seealso: MatTranspose(), MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose(), MatReuse
4476: @*/
4477: PetscErrorCode MatHermitianTranspose(Mat mat,MatReuse reuse,Mat *B)
4478: {
4482: MatTranspose(mat,reuse,B);
4483: #if defined(PETSC_USE_COMPLEX)
4484: MatConjugate(*B);
4485: #endif
4486: return(0);
4487: }
4491: /*@
4492: MatIsHermitianTranspose - Test whether a matrix is another one's Hermitian transpose,
4494: Collective on Mat
4496: Input Parameter:
4497: + A - the matrix to test
4498: - B - the matrix to test against, this can equal the first parameter
4500: Output Parameters:
4501: . flg - the result
4503: Notes:
4504: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
4505: has a running time of the order of the number of nonzeros; the parallel
4506: test involves parallel copies of the block-offdiagonal parts of the matrix.
4508: Level: intermediate
4510: Concepts: matrices^transposing, matrix^symmetry
4512: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian(), MatIsTranspose()
4513: @*/
4514: PetscErrorCode MatIsHermitianTranspose(Mat A,Mat B,PetscReal tol,PetscBool *flg)
4515: {
4516: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscBool *),(*g)(Mat,Mat,PetscReal,PetscBool *);
4522: PetscObjectQueryFunction((PetscObject)A,"MatIsHermitianTranspose_C",(void (**)(void))&f);
4523: PetscObjectQueryFunction((PetscObject)B,"MatIsHermitianTranspose_C",(void (**)(void))&g);
4524: if (f && g) {
4525: if (f==g) {
4526: (*f)(A,B,tol,flg);
4527: } else {
4528: SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for Hermitian test");
4529: }
4530: }
4531: return(0);
4532: }
4536: /*@
4537: MatPermute - Creates a new matrix with rows and columns permuted from the
4538: original.
4540: Collective on Mat
4542: Input Parameters:
4543: + mat - the matrix to permute
4544: . row - row permutation, each processor supplies only the permutation for its rows
4545: - col - column permutation, each processor needs the entire column permutation, that is
4546: this is the same size as the total number of columns in the matrix. It can often
4547: be obtained with ISAllGather() on the row permutation
4549: Output Parameters:
4550: . B - the permuted matrix
4552: Level: advanced
4554: Concepts: matrices^permuting
4556: .seealso: MatGetOrdering(), ISAllGather()
4558: @*/
4559: PetscErrorCode MatPermute(Mat mat,IS row,IS col,Mat *B)
4560: {
4569: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4570: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4571: if (!mat->ops->permute) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatPermute not available for Mat type %s",((PetscObject)mat)->type_name);
4572: MatCheckPreallocated(mat,1);
4574: (*mat->ops->permute)(mat,row,col,B);
4575: PetscObjectStateIncrease((PetscObject)*B);
4576: return(0);
4577: }
4581: /*@
4582: MatEqual - Compares two matrices.
4584: Collective on Mat
4586: Input Parameters:
4587: + A - the first matrix
4588: - B - the second matrix
4590: Output Parameter:
4591: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
4593: Level: intermediate
4595: Concepts: matrices^equality between
4596: @*/
4597: PetscErrorCode MatEqual(Mat A,Mat B,PetscBool *flg)
4598: {
4608: MatCheckPreallocated(B,2);
4609: if (!A->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4610: if (!B->assembled) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4611: 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);
4612: if (!A->ops->equal) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)A)->type_name);
4613: if (!B->ops->equal) SETERRQ1(((PetscObject)A)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)B)->type_name);
4614: 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);
4615: MatCheckPreallocated(A,1);
4617: (*A->ops->equal)(A,B,flg);
4618: return(0);
4619: }
4623: /*@
4624: MatDiagonalScale - Scales a matrix on the left and right by diagonal
4625: matrices that are stored as vectors. Either of the two scaling
4626: matrices can be PETSC_NULL.
4628: Collective on Mat
4630: Input Parameters:
4631: + mat - the matrix to be scaled
4632: . l - the left scaling vector (or PETSC_NULL)
4633: - r - the right scaling vector (or PETSC_NULL)
4635: Notes:
4636: MatDiagonalScale() computes A = LAR, where
4637: L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)
4638: The L scales the rows of the matrix, the R scales the columns of the matrix.
4640: Level: intermediate
4642: Concepts: matrices^diagonal scaling
4643: Concepts: diagonal scaling of matrices
4645: .seealso: MatScale()
4646: @*/
4647: PetscErrorCode MatDiagonalScale(Mat mat,Vec l,Vec r)
4648: {
4654: if (!mat->ops->diagonalscale) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4657: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4658: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4659: MatCheckPreallocated(mat,1);
4661: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4662: (*mat->ops->diagonalscale)(mat,l,r);
4663: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4664: PetscObjectStateIncrease((PetscObject)mat);
4665: #if defined(PETSC_HAVE_CUSP)
4666: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
4667: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
4668: }
4669: #endif
4670: return(0);
4671: }
4675: /*@
4676: MatScale - Scales all elements of a matrix by a given number.
4678: Logically Collective on Mat
4680: Input Parameters:
4681: + mat - the matrix to be scaled
4682: - a - the scaling value
4684: Output Parameter:
4685: . mat - the scaled matrix
4687: Level: intermediate
4689: Concepts: matrices^scaling all entries
4691: .seealso: MatDiagonalScale()
4692: @*/
4693: PetscErrorCode MatScale(Mat mat,PetscScalar a)
4694: {
4700: if (a != (PetscScalar)1.0 && !mat->ops->scale) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4701: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4702: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4704: MatCheckPreallocated(mat,1);
4706: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4707: if (a != (PetscScalar)1.0) {
4708: (*mat->ops->scale)(mat,a);
4709: PetscObjectStateIncrease((PetscObject)mat);
4710: }
4711: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4712: #if defined(PETSC_HAVE_CUSP)
4713: if (mat->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
4714: mat->valid_GPU_matrix = PETSC_CUSP_CPU;
4715: }
4716: #endif
4717: return(0);
4718: }
4722: /*@
4723: MatNorm - Calculates various norms of a matrix.
4725: Collective on Mat
4727: Input Parameters:
4728: + mat - the matrix
4729: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
4731: Output Parameters:
4732: . nrm - the resulting norm
4734: Level: intermediate
4736: Concepts: matrices^norm
4737: Concepts: norm^of matrix
4738: @*/
4739: PetscErrorCode MatNorm(Mat mat,NormType type,PetscReal *nrm)
4740: {
4748: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4749: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4750: if (!mat->ops->norm) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4751: MatCheckPreallocated(mat,1);
4753: (*mat->ops->norm)(mat,type,nrm);
4754: return(0);
4755: }
4757: /*
4758: This variable is used to prevent counting of MatAssemblyBegin() that
4759: are called from within a MatAssemblyEnd().
4760: */
4761: static PetscInt MatAssemblyEnd_InUse = 0;
4764: /*@
4765: MatAssemblyBegin - Begins assembling the matrix. This routine should
4766: be called after completing all calls to MatSetValues().
4768: Collective on Mat
4770: Input Parameters:
4771: + mat - the matrix
4772: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4773:
4774: Notes:
4775: MatSetValues() generally caches the values. The matrix is ready to
4776: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4777: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4778: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4779: using the matrix.
4781: Level: beginner
4783: Concepts: matrices^assembling
4785: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
4786: @*/
4787: PetscErrorCode MatAssemblyBegin(Mat mat,MatAssemblyType type)
4788: {
4794: MatCheckPreallocated(mat,1);
4795: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
4796: if (mat->assembled) {
4797: mat->was_assembled = PETSC_TRUE;
4798: mat->assembled = PETSC_FALSE;
4799: }
4800: if (!MatAssemblyEnd_InUse) {
4801: PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
4802: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4803: PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
4804: } else {
4805: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4806: }
4807: return(0);
4808: }
4812: /*@
4813: MatAssembled - Indicates if a matrix has been assembled and is ready for
4814: use; for example, in matrix-vector product.
4816: Not Collective
4818: Input Parameter:
4819: . mat - the matrix
4821: Output Parameter:
4822: . assembled - PETSC_TRUE or PETSC_FALSE
4824: Level: advanced
4826: Concepts: matrices^assembled?
4828: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
4829: @*/
4830: PetscErrorCode MatAssembled(Mat mat,PetscBool *assembled)
4831: {
4836: *assembled = mat->assembled;
4837: return(0);
4838: }
4842: /*
4843: Processes command line options to determine if/how a matrix
4844: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
4845: */
4846: PetscErrorCode MatView_Private(Mat mat)
4847: {
4848: PetscErrorCode ierr;
4849: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,flg3 = PETSC_FALSE,flg4 = PETSC_FALSE,flg6 = PETSC_FALSE,flg7 = PETSC_FALSE,flg8 = PETSC_FALSE;
4850: static PetscBool incall = PETSC_FALSE;
4851: #if defined(PETSC_USE_SOCKET_VIEWER)
4852: PetscBool flg5 = PETSC_FALSE;
4853: #endif
4856: if (incall) return(0);
4857: incall = PETSC_TRUE;
4858: PetscObjectOptionsBegin((PetscObject)mat);
4859: PetscOptionsBool("-mat_view_info","Information on matrix size","MatView",flg1,&flg1,PETSC_NULL);
4860: PetscOptionsBool("-mat_view_info_detailed","Nonzeros in the matrix","MatView",flg2,&flg2,PETSC_NULL);
4861: PetscOptionsBool("-mat_view","Print matrix to stdout","MatView",flg3,&flg3,PETSC_NULL);
4862: PetscOptionsBool("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",flg4,&flg4,PETSC_NULL);
4863: #if defined(PETSC_USE_SOCKET_VIEWER)
4864: PetscOptionsBool("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",flg5,&flg5,PETSC_NULL);
4865: #endif
4866: PetscOptionsBool("-mat_view_binary","Save matrix to file in binary format","MatView",flg6,&flg6,PETSC_NULL);
4867: PetscOptionsBool("-mat_view_draw","Draw the matrix nonzero structure","MatView",flg7,&flg7,PETSC_NULL);
4868: PetscOptionsEnd();
4870: if (flg1) {
4871: PetscViewer viewer;
4873: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4874: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
4875: MatView(mat,viewer);
4876: PetscViewerPopFormat(viewer);
4877: }
4878: if (flg2) {
4879: PetscViewer viewer;
4881: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4882: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
4883: MatView(mat,viewer);
4884: PetscViewerPopFormat(viewer);
4885: }
4886: if (flg3) {
4887: PetscViewer viewer;
4889: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4890: MatView(mat,viewer);
4891: }
4892: if (flg4) {
4893: PetscViewer viewer;
4895: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4896: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
4897: MatView(mat,viewer);
4898: PetscViewerPopFormat(viewer);
4899: }
4900: #if defined(PETSC_USE_SOCKET_VIEWER)
4901: if (flg5) {
4902: MatView(mat,PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4903: PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4904: }
4905: #endif
4906: if (flg6) {
4907: MatView(mat,PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4908: PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4909: }
4910: if (flg7) {
4911: PetscOptionsGetBool(((PetscObject)mat)->prefix,"-mat_view_contour",&flg8,PETSC_NULL);
4912: if (flg8) {
4913: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm),PETSC_VIEWER_DRAW_CONTOUR);
4914: }
4915: MatView(mat,PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4916: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4917: if (flg8) {
4918: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4919: }
4920: }
4921: incall = PETSC_FALSE;
4922: return(0);
4923: }
4927: /*@
4928: MatAssemblyEnd - Completes assembling the matrix. This routine should
4929: be called after MatAssemblyBegin().
4931: Collective on Mat
4933: Input Parameters:
4934: + mat - the matrix
4935: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4937: Options Database Keys:
4938: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
4939: . -mat_view_info_detailed - Prints more detailed info
4940: . -mat_view - Prints matrix in ASCII format
4941: . -mat_view_matlab - Prints matrix in Matlab format
4942: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
4943: . -display <name> - Sets display name (default is host)
4944: . -draw_pause <sec> - Sets number of seconds to pause after display
4945: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (See the <a href="../../docs/manual.pdf">users manual</a>)
4946: . -viewer_socket_machine <machine>
4947: . -viewer_socket_port <port>
4948: . -mat_view_binary - save matrix to file in binary format
4949: