Actual source code: matrix.c
1: #define PETSCMAT_DLL
3: /*
4: This is where the abstract matrix operations are defined
5: */
7: #include private/matimpl.h
8: #include private/vecimpl.h
10: /* Logging support */
11: PetscCookie MAT_COOKIE;
12: PetscCookie MAT_FDCOLORING_COOKIE;
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_Relax, 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_MatMult, MAT_MatMultSymbolic, MAT_MatMultNumeric;
23: PetscLogEvent MAT_PtAP, MAT_PtAPSymbolic, MAT_PtAPNumeric;
24: PetscLogEvent MAT_MatMultTranspose, MAT_MatMultTransposeSymbolic, MAT_MatMultTransposeNumeric;
25: PetscLogEvent MAT_Getsymtranspose, MAT_Getsymtransreduced, MAT_Transpose_SeqAIJ, MAT_GetBrowsOfAcols;
26: PetscLogEvent MAT_GetBrowsOfAocols, MAT_Getlocalmat, MAT_Getlocalmatcondensed, MAT_Seqstompi, MAT_Seqstompinum, MAT_Seqstompisym;
27: PetscLogEvent MAT_Applypapt, MAT_Applypapt_numeric, MAT_Applypapt_symbolic, MAT_GetSequentialNonzeroStructure;
29: /* nasty global values for MatSetValue() */
30: PetscInt MatSetValue_Row = 0;
31: PetscInt MatSetValue_Column = 0;
32: PetscScalar MatSetValue_Value = 0.0;
36: /*@
37: MatGetDiagonalBlock - Returns the part of the matrix associated with the on-process coupling
39: Not Collective
41: Input Parameters:
42: + mat - the matrix
43: - reuse - indicates you are passing in the a matrix and want it reused
45: Output Parameters:
46: + iscopy - indicates a copy of the diagonal matrix was created and you should use MatDestroy() on it
47: - a - the diagonal part (which is a SEQUENTIAL matrix)
49: Notes: see the manual page for MatCreateMPIAIJ() for more information on the "diagonal part" of the matrix
51: Level: advanced
53: @*/
54: PetscErrorCode MatGetDiagonalBlock(Mat A,PetscTruth *iscopy,MatReuse reuse,Mat *a)
55: {
56: PetscErrorCode ierr,(*f)(Mat,PetscTruth*,MatReuse,Mat*);
57: PetscMPIInt size;
60: MPI_Comm_size(((PetscObject)A)->comm,&size);
61: PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);
62: if (f) {
63: (*f)(A,iscopy,reuse,a);
64: } else if (size == 1) {
65: *a = A;
66: } else {
67: SETERRQ(PETSC_ERR_SUP,"Cannot get diagonal part for this matrix");
68: }
69: return(0);
70: }
74: /*@
75: MatRealPart - Zeros out the imaginary part of the matrix
77: Collective on Mat
79: Input Parameters:
80: . mat - the matrix
82: Level: advanced
85: .seealso: MatImaginaryPart()
86: @*/
87: PetscErrorCode MatRealPart(Mat mat)
88: {
94: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
95: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
96: if (!mat->ops->realpart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
97: MatPreallocated(mat);
98: (*mat->ops->realpart)(mat);
99: return(0);
100: }
105: /*@
106: MatImaginaryPart - Moves the imaginary part of the matrix to the real part and zeros the imaginary part
108: Collective on Mat
110: Input Parameters:
111: . mat - the matrix
113: Level: advanced
116: .seealso: MatRealPart()
117: @*/
118: PetscErrorCode MatImaginaryPart(Mat mat)
119: {
125: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
126: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
127: if (!mat->ops->imaginarypart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
128: MatPreallocated(mat);
129: (*mat->ops->imaginarypart)(mat);
130: return(0);
131: }
135: /*@
136: MatMissingDiagonal - Determine if sparse matrix is missing a diagonal entry (or block entry for BAIJ matrices)
138: Collective on Mat
140: Input Parameter:
141: . mat - the matrix
143: Output Parameters:
144: + missing - is any diagonal missing
145: - dd - first diagonal entry that is missing (optional)
147: Level: advanced
150: .seealso: MatRealPart()
151: @*/
152: PetscErrorCode MatMissingDiagonal(Mat mat,PetscTruth *missing,PetscInt *dd)
153: {
159: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
160: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
161: if (!mat->ops->missingdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
162: (*mat->ops->missingdiagonal)(mat,missing,dd);
163: return(0);
164: }
168: /*@C
169: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
170: for each row that you get to ensure that your application does
171: not bleed memory.
173: Not Collective
175: Input Parameters:
176: + mat - the matrix
177: - row - the row to get
179: Output Parameters:
180: + ncols - if not NULL, the number of nonzeros in the row
181: . cols - if not NULL, the column numbers
182: - vals - if not NULL, the values
184: Notes:
185: This routine is provided for people who need to have direct access
186: to the structure of a matrix. We hope that we provide enough
187: high-level matrix routines that few users will need it.
189: MatGetRow() always returns 0-based column indices, regardless of
190: whether the internal representation is 0-based (default) or 1-based.
192: For better efficiency, set cols and/or vals to PETSC_NULL if you do
193: not wish to extract these quantities.
195: The user can only examine the values extracted with MatGetRow();
196: the values cannot be altered. To change the matrix entries, one
197: must use MatSetValues().
199: You can only have one call to MatGetRow() outstanding for a particular
200: matrix at a time, per processor. MatGetRow() can only obtain rows
201: associated with the given processor, it cannot get rows from the
202: other processors; for that we suggest using MatGetSubMatrices(), then
203: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
204: is in the global number of rows.
206: Fortran Notes:
207: The calling sequence from Fortran is
208: .vb
209: MatGetRow(matrix,row,ncols,cols,values,ierr)
210: Mat matrix (input)
211: integer row (input)
212: integer ncols (output)
213: integer cols(maxcols) (output)
214: double precision (or double complex) values(maxcols) output
215: .ve
216: where maxcols >= maximum nonzeros in any row of the matrix.
219: Caution:
220: Do not try to change the contents of the output arrays (cols and vals).
221: In some cases, this may corrupt the matrix.
223: Level: advanced
225: Concepts: matrices^row access
227: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubMatrices(), MatGetDiagonal()
228: @*/
229: PetscErrorCode MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
230: {
232: PetscInt incols;
237: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
238: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
239: if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
240: MatPreallocated(mat);
241: PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
242: (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
243: if (ncols) *ncols = incols;
244: PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
245: return(0);
246: }
250: /*@
251: MatConjugate - replaces the matrix values with their complex conjugates
253: Collective on Mat
255: Input Parameters:
256: . mat - the matrix
258: Level: advanced
260: .seealso: VecConjugate()
261: @*/
262: PetscErrorCode MatConjugate(Mat mat)
263: {
268: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
269: if (!mat->ops->conjugate) SETERRQ(PETSC_ERR_SUP,"Not provided for this matrix format, send email to petsc-maint@mcs.anl.gov");
270: (*mat->ops->conjugate)(mat);
271: return(0);
272: }
276: /*@C
277: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
279: Not Collective
281: Input Parameters:
282: + mat - the matrix
283: . row - the row to get
284: . ncols, cols - the number of nonzeros and their columns
285: - vals - if nonzero the column values
287: Notes:
288: This routine should be called after you have finished examining the entries.
290: Fortran Notes:
291: The calling sequence from Fortran is
292: .vb
293: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
294: Mat matrix (input)
295: integer row (input)
296: integer ncols (output)
297: integer cols(maxcols) (output)
298: double precision (or double complex) values(maxcols) output
299: .ve
300: Where maxcols >= maximum nonzeros in any row of the matrix.
302: In Fortran MatRestoreRow() MUST be called after MatGetRow()
303: before another call to MatGetRow() can be made.
305: Level: advanced
307: .seealso: MatGetRow()
308: @*/
309: PetscErrorCode MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
310: {
316: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
317: if (!mat->ops->restorerow) return(0);
318: (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
319: return(0);
320: }
324: /*@
325: MatGetRowUpperTriangular - Sets a flag to enable calls to MatGetRow() for matrix in MATSBAIJ format.
326: You should call MatRestoreRowUpperTriangular() after calling MatGetRow/MatRestoreRow() to disable the flag.
328: Not Collective
330: Input Parameters:
331: + mat - the matrix
333: Notes:
334: 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.
336: Level: advanced
338: Concepts: matrices^row access
340: .seealso: MatRestoreRowRowUpperTriangular()
341: @*/
342: PetscErrorCode MatGetRowUpperTriangular(Mat mat)
343: {
349: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
350: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
351: if (!mat->ops->getrowuppertriangular) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
352: MatPreallocated(mat);
353: (*mat->ops->getrowuppertriangular)(mat);
354: return(0);
355: }
359: /*@
360: MatRestoreRowUpperTriangular - Disable calls to MatGetRow() for matrix in MATSBAIJ format.
362: Not Collective
364: Input Parameters:
365: + mat - the matrix
367: Notes:
368: This routine should be called after you have finished MatGetRow/MatRestoreRow().
371: Level: advanced
373: .seealso: MatGetRowUpperTriangular()
374: @*/
375: PetscErrorCode MatRestoreRowUpperTriangular(Mat mat)
376: {
381: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
382: if (!mat->ops->restorerowuppertriangular) return(0);
383: (*mat->ops->restorerowuppertriangular)(mat);
384: return(0);
385: }
389: /*@C
390: MatSetOptionsPrefix - Sets the prefix used for searching for all
391: Mat options in the database.
393: Collective on Mat
395: Input Parameter:
396: + A - the Mat context
397: - prefix - the prefix to prepend to all option names
399: Notes:
400: A hyphen (-) must NOT be given at the beginning of the prefix name.
401: The first character of all runtime options is AUTOMATICALLY the hyphen.
403: Level: advanced
405: .keywords: Mat, set, options, prefix, database
407: .seealso: MatSetFromOptions()
408: @*/
409: PetscErrorCode MatSetOptionsPrefix(Mat A,const char prefix[])
410: {
415: PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
416: return(0);
417: }
421: /*@C
422: MatAppendOptionsPrefix - Appends to the prefix used for searching for all
423: Mat options in the database.
425: Collective on Mat
427: Input Parameters:
428: + A - the Mat context
429: - prefix - the prefix to prepend to all option names
431: Notes:
432: A hyphen (-) must NOT be given at the beginning of the prefix name.
433: The first character of all runtime options is AUTOMATICALLY the hyphen.
435: Level: advanced
437: .keywords: Mat, append, options, prefix, database
439: .seealso: MatGetOptionsPrefix()
440: @*/
441: PetscErrorCode MatAppendOptionsPrefix(Mat A,const char prefix[])
442: {
444:
447: PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
448: return(0);
449: }
453: /*@C
454: MatGetOptionsPrefix - Sets the prefix used for searching for all
455: Mat options in the database.
457: Not Collective
459: Input Parameter:
460: . A - the Mat context
462: Output Parameter:
463: . prefix - pointer to the prefix string used
465: Notes: On the fortran side, the user should pass in a string 'prefix' of
466: sufficient length to hold the prefix.
468: Level: advanced
470: .keywords: Mat, get, options, prefix, database
472: .seealso: MatAppendOptionsPrefix()
473: @*/
474: PetscErrorCode MatGetOptionsPrefix(Mat A,const char *prefix[])
475: {
480: PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
481: return(0);
482: }
486: /*@
487: MatSetUp - Sets up the internal matrix data structures for the later use.
489: Collective on Mat
491: Input Parameters:
492: . A - the Mat context
494: Notes:
495: For basic use of the Mat classes the user need not explicitly call
496: MatSetUp(), since these actions will happen automatically.
498: Level: advanced
500: .keywords: Mat, setup
502: .seealso: MatCreate(), MatDestroy()
503: @*/
504: PetscErrorCode MatSetUp(Mat A)
505: {
506: PetscMPIInt size;
511: if (!((PetscObject)A)->type_name) {
512: MPI_Comm_size(((PetscObject)A)->comm, &size);
513: if (size == 1) {
514: MatSetType(A, MATSEQAIJ);
515: } else {
516: MatSetType(A, MATMPIAIJ);
517: }
518: }
519: MatSetUpPreallocation(A);
520: return(0);
521: }
525: /*@C
526: MatView - Visualizes a matrix object.
528: Collective on Mat
530: Input Parameters:
531: + mat - the matrix
532: - viewer - visualization context
534: Notes:
535: The available visualization contexts include
536: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
537: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
538: output where only the first processor opens
539: the file. All other processors send their
540: data to the first processor to print.
541: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
543: The user can open alternative visualization contexts with
544: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
545: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
546: specified file; corresponding input uses MatLoad()
547: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
548: an X window display
549: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
550: Currently only the sequential dense and AIJ
551: matrix types support the Socket viewer.
553: The user can call PetscViewerSetFormat() to specify the output
554: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
555: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
556: + PETSC_VIEWER_DEFAULT - default, prints matrix contents
557: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
558: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
559: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
560: format common among all matrix types
561: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
562: format (which is in many cases the same as the default)
563: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
564: size and structure (not the matrix entries)
565: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
566: the matrix structure
568: Options Database Keys:
569: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
570: . -mat_view_info_detailed - Prints more detailed info
571: . -mat_view - Prints matrix in ASCII format
572: . -mat_view_matlab - Prints matrix in Matlab format
573: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
574: . -display <name> - Sets display name (default is host)
575: . -draw_pause <sec> - Sets number of seconds to pause after display
576: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
577: . -viewer_socket_machine <machine>
578: . -viewer_socket_port <port>
579: . -mat_view_binary - save matrix to file in binary format
580: - -viewer_binary_filename <name>
581: Level: beginner
583: Notes: see the manual page for MatLoad() for the exact format of the binary file when the binary
584: viewer is used.
586: See bin/matlab/PetscBinaryRead.m for a Matlab code that can read in the binary file when the binary
587: viewer is used.
589: Concepts: matrices^viewing
590: Concepts: matrices^plotting
591: Concepts: matrices^printing
593: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
594: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
595: @*/
596: PetscErrorCode MatView(Mat mat,PetscViewer viewer)
597: {
598: PetscErrorCode ierr;
599: PetscInt rows,cols;
600: PetscTruth iascii;
601: const MatType cstr;
602: PetscViewerFormat format;
607: if (!viewer) {
608: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
609: }
612: if (!mat->assembled) SETERRQ(PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
613: MatPreallocated(mat);
615: PetscLogEventBegin(MAT_View,mat,viewer,0,0);
616: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
617: if (iascii) {
618: PetscViewerGetFormat(viewer,&format);
619: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
620: if (((PetscObject)mat)->prefix) {
621: PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)\n",((PetscObject)mat)->prefix);
622: } else {
623: PetscViewerASCIIPrintf(viewer,"Matrix Object:\n");
624: }
625: PetscViewerASCIIPushTab(viewer);
626: MatGetType(mat,&cstr);
627: MatGetSize(mat,&rows,&cols);
628: PetscViewerASCIIPrintf(viewer,"type=%s, rows=%D, cols=%D\n",cstr,rows,cols);
629: if (mat->factor) {
630: const MatSolverPackage solver;
631: MatFactorGetSolverPackage(mat,&solver);
632: PetscViewerASCIIPrintf(viewer,"package used to perform factorization: %s\n",solver);
633: }
634: if (mat->ops->getinfo) {
635: MatInfo info;
636: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
637: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",(PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
638: }
639: }
640: }
641: if (mat->ops->view) {
642: PetscViewerASCIIPushTab(viewer);
643: (*mat->ops->view)(mat,viewer);
644: PetscViewerASCIIPopTab(viewer);
645: } else if (!iascii) {
646: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
647: }
648: if (iascii) {
649: PetscViewerGetFormat(viewer,&format);
650: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
651: PetscViewerASCIIPopTab(viewer);
652: }
653: }
654: PetscLogEventEnd(MAT_View,mat,viewer,0,0);
655: return(0);
656: }
660: /*@
661: MatScaleSystem - Scale a vector solution and right hand side to
662: match the scaling of a scaled matrix.
663:
664: Collective on Mat
666: Input Parameter:
667: + mat - the matrix
668: . b - right hand side vector (or PETSC_NULL)
669: - x - solution vector (or PETSC_NULL)
672: Notes:
673: For AIJ, and BAIJ matrix formats, the matrices are not
674: internally scaled, so this does nothing. For MPIROWBS it
675: permutes and diagonally scales.
677: The KSP methods automatically call this routine when required
678: (via PCPreSolve()) so it is rarely used directly.
680: Level: Developer
682: Concepts: matrices^scaling
684: .seealso: MatUseScaledForm(), MatUnScaleSystem()
685: @*/
686: PetscErrorCode MatScaleSystem(Mat mat,Vec b,Vec x)
687: {
693: MatPreallocated(mat);
697: if (mat->ops->scalesystem) {
698: (*mat->ops->scalesystem)(mat,b,x);
699: }
700: return(0);
701: }
705: /*@
706: MatUnScaleSystem - Unscales a vector solution and right hand side to
707: match the original scaling of a scaled matrix.
708:
709: Collective on Mat
711: Input Parameter:
712: + mat - the matrix
713: . b - right hand side vector (or PETSC_NULL)
714: - x - solution vector (or PETSC_NULL)
717: Notes:
718: For AIJ and BAIJ matrix formats, the matrices are not
719: internally scaled, so this does nothing. For MPIROWBS it
720: permutes and diagonally scales.
722: The KSP methods automatically call this routine when required
723: (via PCPreSolve()) so it is rarely used directly.
725: Level: Developer
727: .seealso: MatUseScaledForm(), MatScaleSystem()
728: @*/
729: PetscErrorCode MatUnScaleSystem(Mat mat,Vec b,Vec x)
730: {
736: MatPreallocated(mat);
739: if (mat->ops->unscalesystem) {
740: (*mat->ops->unscalesystem)(mat,b,x);
741: }
742: return(0);
743: }
747: /*@
748: MatUseScaledForm - For matrix storage formats that scale the
749: matrix (for example MPIRowBS matrices are diagonally scaled on
750: assembly) indicates matrix operations (MatMult() etc) are
751: applied using the scaled matrix.
752:
753: Collective on Mat
755: Input Parameter:
756: + mat - the matrix
757: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
758: applying the original matrix
760: Notes:
761: For scaled matrix formats, applying the original, unscaled matrix
762: will be slightly more expensive
764: Level: Developer
766: .seealso: MatScaleSystem(), MatUnScaleSystem()
767: @*/
768: PetscErrorCode MatUseScaledForm(Mat mat,PetscTruth scaled)
769: {
775: MatPreallocated(mat);
776: if (mat->ops->usescaledform) {
777: (*mat->ops->usescaledform)(mat,scaled);
778: }
779: return(0);
780: }
784: /*@
785: MatDestroy - Frees space taken by a matrix.
786:
787: Collective on Mat
789: Input Parameter:
790: . A - the matrix
792: Level: beginner
794: @*/
795: PetscErrorCode MatDestroy(Mat A)
796: {
800: if (--((PetscObject)A)->refct > 0) return(0);
801: MatPreallocated(A);
802: /* if memory was published with AMS then destroy it */
803: PetscObjectDepublish(A);
804: if (A->ops->destroy) {
805: (*A->ops->destroy)(A);
806: }
807: if (A->mapping) {
808: ISLocalToGlobalMappingDestroy(A->mapping);
809: }
810: if (A->bmapping) {
811: ISLocalToGlobalMappingDestroy(A->bmapping);
812: }
814: if (A->spptr){PetscFree(A->spptr);}
815: PetscMapDestroy(A->rmap);
816: PetscMapDestroy(A->cmap);
817: PetscHeaderDestroy(A);
818: return(0);
819: }
823: /*@
824: MatValid - Checks whether a matrix object is valid.
826: Collective on Mat
828: Input Parameter:
829: . m - the matrix to check
831: Output Parameter:
832: flg - flag indicating matrix status, either
833: PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
835: Level: developer
837: Concepts: matrices^validity
838: @*/
839: PetscErrorCode MatValid(Mat m,PetscTruth *flg)
840: {
843: if (!m) *flg = PETSC_FALSE;
844: else if (((PetscObject)m)->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
845: else *flg = PETSC_TRUE;
846: return(0);
847: }
851: /*@
852: MatSetValues - Inserts or adds a block of values into a matrix.
853: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
854: MUST be called after all calls to MatSetValues() have been completed.
856: Not Collective
858: Input Parameters:
859: + mat - the matrix
860: . v - a logically two-dimensional array of values
861: . m, idxm - the number of rows and their global indices
862: . n, idxn - the number of columns and their global indices
863: - addv - either ADD_VALUES or INSERT_VALUES, where
864: ADD_VALUES adds values to any existing entries, and
865: INSERT_VALUES replaces existing entries with new values
867: Notes:
868: By default the values, v, are row-oriented and unsorted.
869: See MatSetOption() for other options.
871: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
872: options cannot be mixed without intervening calls to the assembly
873: routines.
875: MatSetValues() uses 0-based row and column numbers in Fortran
876: as well as in C.
878: Negative indices may be passed in idxm and idxn, these rows and columns are
879: simply ignored. This allows easily inserting element stiffness matrices
880: with homogeneous Dirchlet boundary conditions that you don't want represented
881: in the matrix.
883: Efficiency Alert:
884: The routine MatSetValuesBlocked() may offer much better efficiency
885: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
887: Level: beginner
889: Concepts: matrices^putting entries in
891: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
892: InsertMode, INSERT_VALUES, ADD_VALUES
893: @*/
894: PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
895: {
901: if (!m || !n) return(0); /* no values to insert */
904: MatPreallocated(mat);
905: if (mat->insertmode == NOT_SET_VALUES) {
906: mat->insertmode = addv;
907: }
908: #if defined(PETSC_USE_DEBUG)
909: else if (mat->insertmode != addv) {
910: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
911: }
912: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
913: #endif
915: if (mat->assembled) {
916: mat->was_assembled = PETSC_TRUE;
917: mat->assembled = PETSC_FALSE;
918: }
919: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
920: if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
921: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
922: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
923: return(0);
924: }
929: /*@
930: MatSetValuesRowLocal - Inserts a row (block row for BAIJ matrices) of nonzero
931: values into a matrix
933: Not Collective
935: Input Parameters:
936: + mat - the matrix
937: . row - the (block) row to set
938: - v - a logically two-dimensional array of values
940: Notes:
941: By the values, v, are column-oriented (for the block version) and sorted
943: All the nonzeros in the row must be provided
945: The matrix must have previously had its column indices set
947: The row must belong to this process
949: Level: intermediate
951: Concepts: matrices^putting entries in
953: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
954: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues(), MatSetValuesRow(), MatSetLocalToGlobalMapping()
955: @*/
956: PetscErrorCode MatSetValuesRowLocal(Mat mat,PetscInt row,const PetscScalar v[])
957: {
964: MatSetValuesRow(mat, mat->mapping->indices[row],v);
965: return(0);
966: }
970: /*@
971: MatSetValuesRow - Inserts a row (block row for BAIJ matrices) of nonzero
972: values into a matrix
974: Not Collective
976: Input Parameters:
977: + mat - the matrix
978: . row - the (block) row to set
979: - v - a logically two-dimensional array of values
981: Notes:
982: By the values, v, are column-oriented (for the block version) and sorted
984: All the nonzeros in the row must be provided
986: The matrix must have previously had its column indices set
988: The row must belong to this process
990: Level: intermediate
992: Concepts: matrices^putting entries in
994: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
995: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
996: @*/
997: PetscErrorCode MatSetValuesRow(Mat mat,PetscInt row,const PetscScalar v[])
998: {
1005: #if defined(PETSC_USE_DEBUG)
1006: if (mat->insertmode == ADD_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add and insert values");
1007: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1008: #endif
1009: mat->insertmode = INSERT_VALUES;
1011: if (mat->assembled) {
1012: mat->was_assembled = PETSC_TRUE;
1013: mat->assembled = PETSC_FALSE;
1014: }
1015: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1016: if (!mat->ops->setvaluesrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1017: (*mat->ops->setvaluesrow)(mat,row,v);
1018: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1019: return(0);
1020: }
1024: /*@
1025: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
1026: Using structured grid indexing
1028: Not Collective
1030: Input Parameters:
1031: + mat - the matrix
1032: . v - a logically two-dimensional array of values
1033: . m - number of rows being entered
1034: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
1035: . n - number of columns being entered
1036: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
1037: - addv - either ADD_VALUES or INSERT_VALUES, where
1038: ADD_VALUES adds values to any existing entries, and
1039: INSERT_VALUES replaces existing entries with new values
1041: Notes:
1042: By default the values, v, are row-oriented and unsorted.
1043: See MatSetOption() for other options.
1045: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
1046: options cannot be mixed without intervening calls to the assembly
1047: routines.
1049: The grid coordinates are across the entire grid, not just the local portion
1051: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
1052: as well as in C.
1054: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
1056: In order to use this routine you must either obtain the matrix with DAGetMatrix()
1057: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1059: The columns and rows in the stencil passed in MUST be contained within the
1060: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
1061: if you create a DA with an overlap of one grid level and on a particular process its first
1062: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1063: first i index you can use in your column and row indices in MatSetStencil() is 5.
1065: In Fortran idxm and idxn should be declared as
1066: $ MatStencil idxm(4,m),idxn(4,n)
1067: and the values inserted using
1068: $ idxm(MatStencil_i,1) = i
1069: $ idxm(MatStencil_j,1) = j
1070: $ idxm(MatStencil_k,1) = k
1071: $ idxm(MatStencil_c,1) = c
1072: etc
1074: For periodic boundary conditions use negative indices for values to the left (below 0; that are to be
1075: obtained by wrapping values from right edge). For values to the right of the last entry using that index plus one
1076: etc to obtain values that obtained by wrapping the values from the left edge. This does not work for the DA_NONPERIODIC
1077: wrap.
1079: 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
1080: a single value per point) you can skip filling those indices.
1082: Inspired by the structured grid interface to the HYPRE package
1083: (http://www.llnl.gov/CASC/hypre)
1085: Efficiency Alert:
1086: The routine MatSetValuesBlockedStencil() may offer much better efficiency
1087: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1089: Level: beginner
1091: Concepts: matrices^putting entries in
1093: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1094: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1095: @*/
1096: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1097: {
1099: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1100: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1103: if (!m || !n) return(0); /* no values to insert */
1110: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1111: if (n > 256) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1113: for (i=0; i<m; i++) {
1114: for (j=0; j<3-sdim; j++) dxm++;
1115: tmp = *dxm++ - starts[0];
1116: for (j=0; j<dim-1; j++) {
1117: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1118: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1119: }
1120: if (mat->stencil.noc) dxm++;
1121: jdxm[i] = tmp;
1122: }
1123: for (i=0; i<n; i++) {
1124: for (j=0; j<3-sdim; j++) dxn++;
1125: tmp = *dxn++ - starts[0];
1126: for (j=0; j<dim-1; j++) {
1127: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1128: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1129: }
1130: if (mat->stencil.noc) dxn++;
1131: jdxn[i] = tmp;
1132: }
1133: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1134: return(0);
1135: }
1139: /*@C
1140: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1141: Using structured grid indexing
1143: Not Collective
1145: Input Parameters:
1146: + mat - the matrix
1147: . v - a logically two-dimensional array of values
1148: . m - number of rows being entered
1149: . idxm - grid coordinates for matrix rows being entered
1150: . n - number of columns being entered
1151: . idxn - grid coordinates for matrix columns being entered
1152: - addv - either ADD_VALUES or INSERT_VALUES, where
1153: ADD_VALUES adds values to any existing entries, and
1154: INSERT_VALUES replaces existing entries with new values
1156: Notes:
1157: By default the values, v, are row-oriented and unsorted.
1158: See MatSetOption() for other options.
1160: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1161: options cannot be mixed without intervening calls to the assembly
1162: routines.
1164: The grid coordinates are across the entire grid, not just the local portion
1166: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1167: as well as in C.
1169: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
1171: In order to use this routine you must either obtain the matrix with DAGetMatrix()
1172: or call MatSetBlockSize(), MatSetLocalToGlobalMapping() and MatSetStencil() first.
1174: The columns and rows in the stencil passed in MUST be contained within the
1175: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
1176: if you create a DA with an overlap of one grid level and on a particular process its first
1177: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1178: first i index you can use in your column and row indices in MatSetStencil() is 5.
1180: In Fortran idxm and idxn should be declared as
1181: $ MatStencil idxm(4,m),idxn(4,n)
1182: and the values inserted using
1183: $ idxm(MatStencil_i,1) = i
1184: $ idxm(MatStencil_j,1) = j
1185: $ idxm(MatStencil_k,1) = k
1186: etc
1188: Negative indices may be passed in idxm and idxn, these rows and columns are
1189: simply ignored. This allows easily inserting element stiffness matrices
1190: with homogeneous Dirchlet boundary conditions that you don't want represented
1191: in the matrix.
1193: Inspired by the structured grid interface to the HYPRE package
1194: (http://www.llnl.gov/CASC/hypre)
1196: Level: beginner
1198: Concepts: matrices^putting entries in
1200: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1201: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil,
1202: MatSetBlockSize(), MatSetLocalToGlobalMapping()
1203: @*/
1204: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1205: {
1207: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1208: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1211: if (!m || !n) return(0); /* no values to insert */
1218: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1219: if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1221: for (i=0; i<m; i++) {
1222: for (j=0; j<3-sdim; j++) dxm++;
1223: tmp = *dxm++ - starts[0];
1224: for (j=0; j<sdim-1; j++) {
1225: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1226: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1227: }
1228: dxm++;
1229: jdxm[i] = tmp;
1230: }
1231: for (i=0; i<n; i++) {
1232: for (j=0; j<3-sdim; j++) dxn++;
1233: tmp = *dxn++ - starts[0];
1234: for (j=0; j<sdim-1; j++) {
1235: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1236: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1237: }
1238: dxn++;
1239: jdxn[i] = tmp;
1240: }
1241: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1242: return(0);
1243: }
1247: /*@
1248: MatSetStencil - Sets the grid information for setting values into a matrix via
1249: MatSetValuesStencil()
1251: Not Collective
1253: Input Parameters:
1254: + mat - the matrix
1255: . dim - dimension of the grid 1, 2, or 3
1256: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1257: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1258: - dof - number of degrees of freedom per node
1261: Inspired by the structured grid interface to the HYPRE package
1262: (www.llnl.gov/CASC/hyper)
1264: For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the
1265: user.
1266:
1267: Level: beginner
1269: Concepts: matrices^putting entries in
1271: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1272: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1273: @*/
1274: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1275: {
1276: PetscInt i;
1283: mat->stencil.dim = dim + (dof > 1);
1284: for (i=0; i<dim; i++) {
1285: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1286: mat->stencil.starts[i] = starts[dim-i-1];
1287: }
1288: mat->stencil.dims[dim] = dof;
1289: mat->stencil.starts[dim] = 0;
1290: mat->stencil.noc = (PetscTruth)(dof == 1);
1291: return(0);
1292: }
1296: /*@
1297: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1299: Not Collective
1301: Input Parameters:
1302: + mat - the matrix
1303: . v - a logically two-dimensional array of values
1304: . m, idxm - the number of block rows and their global block indices
1305: . n, idxn - the number of block columns and their global block indices
1306: - addv - either ADD_VALUES or INSERT_VALUES, where
1307: ADD_VALUES adds values to any existing entries, and
1308: INSERT_VALUES replaces existing entries with new values
1310: Notes:
1311: The m and n count the NUMBER of blocks in the row direction and column direction,
1312: NOT the total number of rows/columns; for example, if the block size is 2 and
1313: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1314: The values in idxm would be 1 2; that is the first index for each block divided by
1315: the block size.
1317: Note that you must call MatSetBlockSize() when constructing this matrix (and before
1318: preallocating it).
1320: By default the values, v, are row-oriented and unsorted. So the layout of
1321: v is the same as for MatSetValues(). See MatSetOption() for other options.
1323: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1324: options cannot be mixed without intervening calls to the assembly
1325: routines.
1327: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1328: as well as in C.
1330: Negative indices may be passed in idxm and idxn, these rows and columns are
1331: simply ignored. This allows easily inserting element stiffness matrices
1332: with homogeneous Dirchlet boundary conditions that you don't want represented
1333: in the matrix.
1335: Each time an entry is set within a sparse matrix via MatSetValues(),
1336: internal searching must be done to determine where to place the the
1337: data in the matrix storage space. By instead inserting blocks of
1338: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1339: reduced.
1341: Example:
1342: $ Suppose m=n=2 and block size(bs) = 2 The array is
1343: $
1344: $ 1 2 | 3 4
1345: $ 5 6 | 7 8
1346: $ - - - | - - -
1347: $ 9 10 | 11 12
1348: $ 13 14 | 15 16
1349: $
1350: $ v[] should be passed in like
1351: $ v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
1352: $
1353: $ If you are not using row oriented storage of v (that is you called MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE)) then
1354: $ v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
1356: Level: intermediate
1358: Concepts: matrices^putting entries in blocked
1360: .seealso: MatSetBlockSize(), MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1361: @*/
1362: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1363: {
1369: if (!m || !n) return(0); /* no values to insert */
1373: MatPreallocated(mat);
1374: if (mat->insertmode == NOT_SET_VALUES) {
1375: mat->insertmode = addv;
1376: }
1377: #if defined(PETSC_USE_DEBUG)
1378: else if (mat->insertmode != addv) {
1379: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1380: }
1381: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1382: #endif
1384: if (mat->assembled) {
1385: mat->was_assembled = PETSC_TRUE;
1386: mat->assembled = PETSC_FALSE;
1387: }
1388: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1389: if (mat->ops->setvaluesblocked) {
1390: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1391: } else {
1392: PetscInt buf[4096],*ibufm=0,*ibufn=0;
1393: PetscInt i,j,*iidxm,*iidxn,bs=mat->rmap->bs;
1394: if ((m+n)*bs <= 4096) {
1395: iidxm = buf; iidxn = buf + m*bs;
1396: } else {
1397: PetscMalloc2(m*bs,PetscInt,&ibufm,n*bs,PetscInt,&ibufn);
1398: iidxm = ibufm; iidxn = ibufn;
1399: }
1400: for (i=0; i<m; i++) {
1401: for (j=0; j<bs; j++) {
1402: iidxm[i*bs+j] = bs*idxm[i] + j;
1403: }
1404: }
1405: for (i=0; i<n; i++) {
1406: for (j=0; j<bs; j++) {
1407: iidxn[i*bs+j] = bs*idxn[i] + j;
1408: }
1409: }
1410: MatSetValues(mat,bs*m,iidxm,bs*n,iidxn,v,addv);
1411: PetscFree2(ibufm,ibufn);
1412: }
1413: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1414: return(0);
1415: }
1419: /*@
1420: MatGetValues - Gets a block of values from a matrix.
1422: Not Collective; currently only returns a local block
1424: Input Parameters:
1425: + mat - the matrix
1426: . v - a logically two-dimensional array for storing the values
1427: . m, idxm - the number of rows and their global indices
1428: - n, idxn - the number of columns and their global indices
1430: Notes:
1431: The user must allocate space (m*n PetscScalars) for the values, v.
1432: The values, v, are then returned in a row-oriented format,
1433: analogous to that used by default in MatSetValues().
1435: MatGetValues() uses 0-based row and column numbers in
1436: Fortran as well as in C.
1438: MatGetValues() requires that the matrix has been assembled
1439: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1440: MatSetValues() and MatGetValues() CANNOT be made in succession
1441: without intermediate matrix assembly.
1443: Negative row or column indices will be ignored and those locations in v[] will be
1444: left unchanged.
1446: Level: advanced
1448: Concepts: matrices^accessing values
1450: .seealso: MatGetRow(), MatGetSubMatrices(), MatSetValues()
1451: @*/
1452: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1453: {
1462: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1463: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1464: if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1465: MatPreallocated(mat);
1467: PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
1468: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1469: PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
1470: return(0);
1471: }
1475: /*@
1476: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1477: the routine MatSetValuesLocal() to allow users to insert matrix entries
1478: using a local (per-processor) numbering.
1480: Not Collective
1482: Input Parameters:
1483: + x - the matrix
1484: - mapping - mapping created with ISLocalToGlobalMappingCreate()
1485: or ISLocalToGlobalMappingCreateIS()
1487: Level: intermediate
1489: Concepts: matrices^local to global mapping
1490: Concepts: local to global mapping^for matrices
1492: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1493: @*/
1494: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
1495: {
1501: if (x->mapping) {
1502: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1503: }
1504: MatPreallocated(x);
1506: if (x->ops->setlocaltoglobalmapping) {
1507: (*x->ops->setlocaltoglobalmapping)(x,mapping);
1508: } else {
1509: PetscObjectReference((PetscObject)mapping);
1510: if (x->mapping) { ISLocalToGlobalMappingDestroy(x->mapping); }
1511: x->mapping = mapping;
1512: }
1513: return(0);
1514: }
1518: /*@
1519: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
1520: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
1521: entries using a local (per-processor) numbering.
1523: Not Collective
1525: Input Parameters:
1526: + x - the matrix
1527: - mapping - mapping created with ISLocalToGlobalMappingCreate() or
1528: ISLocalToGlobalMappingCreateIS()
1530: Level: intermediate
1532: Concepts: matrices^local to global mapping blocked
1533: Concepts: local to global mapping^for matrices, blocked
1535: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
1536: MatSetValuesBlocked(), MatSetValuesLocal()
1537: @*/
1538: PetscErrorCode MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
1539: {
1545: if (x->bmapping) {
1546: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1547: }
1548: PetscObjectReference((PetscObject)mapping);
1549: if (x->bmapping) { ISLocalToGlobalMappingDestroy(x->mapping); }
1550: x->bmapping = mapping;
1551: return(0);
1552: }
1556: /*@
1557: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
1558: using a local ordering of the nodes.
1560: Not Collective
1562: Input Parameters:
1563: + x - the matrix
1564: . nrow, irow - number of rows and their local indices
1565: . ncol, icol - number of columns and their local indices
1566: . y - a logically two-dimensional array of values
1567: - addv - either INSERT_VALUES or ADD_VALUES, where
1568: ADD_VALUES adds values to any existing entries, and
1569: INSERT_VALUES replaces existing entries with new values
1571: Notes:
1572: Before calling MatSetValuesLocal(), the user must first set the
1573: local-to-global mapping by calling MatSetLocalToGlobalMapping().
1575: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
1576: options cannot be mixed without intervening calls to the assembly
1577: routines.
1579: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1580: MUST be called after all calls to MatSetValuesLocal() have been completed.
1582: Level: intermediate
1584: Concepts: matrices^putting entries in with local numbering
1586: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
1587: MatSetValueLocal()
1588: @*/
1589: PetscErrorCode MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1590: {
1592: PetscInt irowm[2048],icolm[2048];
1597: if (!nrow || !ncol) return(0); /* no values to insert */
1601: MatPreallocated(mat);
1602: if (mat->insertmode == NOT_SET_VALUES) {
1603: mat->insertmode = addv;
1604: }
1605: #if defined(PETSC_USE_DEBUG)
1606: else if (mat->insertmode != addv) {
1607: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1608: }
1609: if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
1610: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1611: }
1612: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1613: #endif
1615: if (mat->assembled) {
1616: mat->was_assembled = PETSC_TRUE;
1617: mat->assembled = PETSC_FALSE;
1618: }
1619: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1620: if (!mat->ops->setvalueslocal) {
1621: ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
1622: ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
1623: (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
1624: } else {
1625: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1626: }
1627: mat->same_nonzero = PETSC_FALSE;
1628: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1629: return(0);
1630: }
1634: /*@
1635: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1636: using a local ordering of the nodes a block at a time.
1638: Not Collective
1640: Input Parameters:
1641: + x - the matrix
1642: . nrow, irow - number of rows and their local indices
1643: . ncol, icol - number of columns and their local indices
1644: . y - a logically two-dimensional array of values
1645: - addv - either INSERT_VALUES or ADD_VALUES, where
1646: ADD_VALUES adds values to any existing entries, and
1647: INSERT_VALUES replaces existing entries with new values
1649: Notes:
1650: Before calling MatSetValuesBlockedLocal(), the user must first set the
1651: block size using MatSetBlockSize(), and the local-to-global mapping by
1652: calling MatSetLocalToGlobalMappingBlock(), where the mapping MUST be
1653: set for matrix blocks, not for matrix elements.
1655: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1656: options cannot be mixed without intervening calls to the assembly
1657: routines.
1659: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1660: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
1662: Level: intermediate
1664: Concepts: matrices^putting blocked values in with local numbering
1666: .seealso: MatSetBlockSize(), MatSetLocalToGlobalMappingBlock(), MatAssemblyBegin(), MatAssemblyEnd(),
1667: MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
1668: @*/
1669: PetscErrorCode MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1670: {
1672: PetscInt irowm[2048],icolm[2048];
1677: if (!nrow || !ncol) return(0); /* no values to insert */
1681: MatPreallocated(mat);
1682: if (mat->insertmode == NOT_SET_VALUES) {
1683: mat->insertmode = addv;
1684: }
1685: #if defined(PETSC_USE_DEBUG)
1686: else if (mat->insertmode != addv) {
1687: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1688: }
1689: if (!mat->bmapping) {
1690: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
1691: }
1692: if (nrow > 2048 || ncol > 2048) {
1693: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1694: }
1695: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1696: #endif
1698: if (mat->assembled) {
1699: mat->was_assembled = PETSC_TRUE;
1700: mat->assembled = PETSC_FALSE;
1701: }
1702: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1703: ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
1704: ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
1705: if (mat->ops->setvaluesblocked) {
1706: (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
1707: } else {
1708: PetscInt buf[4096],*ibufm=0,*ibufn=0;
1709: PetscInt i,j,*iirowm,*iicolm,bs=mat->rmap->bs;
1710: if ((nrow+ncol)*bs <= 4096) {
1711: iirowm = buf; iicolm = buf + nrow*bs;
1712: } else {
1713: PetscMalloc2(nrow*bs,PetscInt,&ibufm,ncol*bs,PetscInt,&ibufn);
1714: iirowm = ibufm; iicolm = ibufn;
1715: }
1716: for (i=0; i<nrow; i++) {
1717: for (j=0; j<bs; j++) {
1718: iirowm[i*bs+j] = bs*irowm[i] + j;
1719: }
1720: }
1721: for (i=0; i<ncol; i++) {
1722: for (j=0; j<bs; j++) {
1723: iicolm[i*bs+j] = bs*icolm[i] + j;
1724: }
1725: }
1726: MatSetValues(mat,bs*nrow,iirowm,bs*ncol,iicolm,y,addv);
1727: PetscFree2(ibufm,ibufn);
1728: }
1729: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1730: return(0);
1731: }
1733: /* --------------------------------------------------------*/
1736: /*@
1737: MatMult - Computes the matrix-vector product, y = Ax.
1739: Collective on Mat and Vec
1741: Input Parameters:
1742: + mat - the matrix
1743: - x - the vector to be multiplied
1745: Output Parameters:
1746: . y - the result
1748: Notes:
1749: The vectors x and y cannot be the same. I.e., one cannot
1750: call MatMult(A,y,y).
1752: Level: beginner
1754: Concepts: matrix-vector product
1756: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
1757: @*/
1758: PetscErrorCode MatMult(Mat mat,Vec x,Vec y)
1759: {
1768: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1769: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1770: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1771: #ifndef PETSC_HAVE_CONSTRAINTS
1772: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
1773: if (mat->rmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
1774: if (mat->rmap->n != y->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->rmap->n,y->map->n);
1775: #endif
1776: MatPreallocated(mat);
1778: if (mat->nullsp) {
1779: MatNullSpaceRemove(mat->nullsp,x,&x);
1780: }
1782: if (!mat->ops->mult) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
1783: PetscLogEventBegin(MAT_Mult,mat,x,y,0);
1784: (*mat->ops->mult)(mat,x,y);
1785: PetscLogEventEnd(MAT_Mult,mat,x,y,0);
1787: if (mat->nullsp) {
1788: MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
1789: }
1790: PetscObjectStateIncrease((PetscObject)y);
1791: return(0);
1792: }
1796: /*@
1797: MatMultTranspose - Computes matrix transpose times a vector.
1799: Collective on Mat and Vec
1801: Input Parameters:
1802: + mat - the matrix
1803: - x - the vector to be multilplied
1805: Output Parameters:
1806: . y - the result
1808: Notes:
1809: The vectors x and y cannot be the same. I.e., one cannot
1810: call MatMultTranspose(A,y,y).
1812: Level: beginner
1814: Concepts: matrix vector product^transpose
1816: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
1817: @*/
1818: PetscErrorCode MatMultTranspose(Mat mat,Vec x,Vec y)
1819: {
1828: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1829: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1830: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1831: #ifndef PETSC_HAVE_CONSTRAINTS
1832: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
1833: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->cmap->N,y->map->N);
1834: #endif
1835: MatPreallocated(mat);
1837: if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
1838: PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
1839: (*mat->ops->multtranspose)(mat,x,y);
1840: PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
1841: PetscObjectStateIncrease((PetscObject)y);
1842: return(0);
1843: }
1847: /*@
1848: MatMultAdd - Computes v3 = v2 + A * v1.
1850: Collective on Mat and Vec
1852: Input Parameters:
1853: + mat - the matrix
1854: - v1, v2 - the vectors
1856: Output Parameters:
1857: . v3 - the result
1859: Notes:
1860: The vectors v1 and v3 cannot be the same. I.e., one cannot
1861: call MatMultAdd(A,v1,v2,v1).
1863: Level: beginner
1865: Concepts: matrix vector product^addition
1867: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1868: @*/
1869: PetscErrorCode MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1870: {
1880: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1881: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1882: if (mat->cmap->N != v1->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->cmap->N,v1->map->N);
1883: if (mat->rmap->N != v2->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->rmap->N,v2->map->N);
1884: if (mat->rmap->N != v3->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->rmap->N,v3->map->N);
1885: if (mat->rmap->n != v3->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %D %D",mat->rmap->n,v3->map->n);
1886: if (mat->rmap->n != v2->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %D %D",mat->rmap->n,v2->map->n);
1887: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1888: MatPreallocated(mat);
1890: PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1891: (*mat->ops->multadd)(mat,v1,v2,v3);
1892: PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1893: PetscObjectStateIncrease((PetscObject)v3);
1894: return(0);
1895: }
1899: /*@
1900: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1902: Collective on Mat and Vec
1904: Input Parameters:
1905: + mat - the matrix
1906: - v1, v2 - the vectors
1908: Output Parameters:
1909: . v3 - the result
1911: Notes:
1912: The vectors v1 and v3 cannot be the same. I.e., one cannot
1913: call MatMultTransposeAdd(A,v1,v2,v1).
1915: Level: beginner
1917: Concepts: matrix vector product^transpose and addition
1919: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1920: @*/
1921: PetscErrorCode MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1922: {
1932: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1933: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1934: if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1935: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1936: if (mat->rmap->N != v1->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->rmap->N,v1->map->N);
1937: if (mat->cmap->N != v2->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->cmap->N,v2->map->N);
1938: if (mat->cmap->N != v3->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->cmap->N,v3->map->N);
1939: MatPreallocated(mat);
1941: PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1942: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
1943: PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1944: PetscObjectStateIncrease((PetscObject)v3);
1945: return(0);
1946: }
1950: /*@
1951: MatMultConstrained - The inner multiplication routine for a
1952: constrained matrix P^T A P.
1954: Collective on Mat and Vec
1956: Input Parameters:
1957: + mat - the matrix
1958: - x - the vector to be multilplied
1960: Output Parameters:
1961: . y - the result
1963: Notes:
1964: The vectors x and y cannot be the same. I.e., one cannot
1965: call MatMult(A,y,y).
1967: Level: beginner
1969: .keywords: matrix, multiply, matrix-vector product, constraint
1970: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1971: @*/
1972: PetscErrorCode MatMultConstrained(Mat mat,Vec x,Vec y)
1973: {
1980: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1981: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1982: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1983: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
1984: if (mat->rmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
1985: if (mat->rmap->n != y->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->rmap->n,y->map->n);
1987: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1988: (*mat->ops->multconstrained)(mat,x,y);
1989: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
1990: PetscObjectStateIncrease((PetscObject)y);
1992: return(0);
1993: }
1997: /*@
1998: MatMultTransposeConstrained - The inner multiplication routine for a
1999: constrained matrix P^T A^T P.
2001: Collective on Mat and Vec
2003: Input Parameters:
2004: + mat - the matrix
2005: - x - the vector to be multilplied
2007: Output Parameters:
2008: . y - the result
2010: Notes:
2011: The vectors x and y cannot be the same. I.e., one cannot
2012: call MatMult(A,y,y).
2014: Level: beginner
2016: .keywords: matrix, multiply, matrix-vector product, constraint
2017: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
2018: @*/
2019: PetscErrorCode MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
2020: {
2027: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2028: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2029: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2030: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2031: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
2033: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2034: (*mat->ops->multtransposeconstrained)(mat,x,y);
2035: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2036: PetscObjectStateIncrease((PetscObject)y);
2038: return(0);
2039: }
2040: /* ------------------------------------------------------------*/
2043: /*@C
2044: MatGetInfo - Returns information about matrix storage (number of
2045: nonzeros, memory, etc.).
2047: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
2048: as the flag
2050: Input Parameters:
2051: . mat - the matrix
2053: Output Parameters:
2054: + flag - flag indicating the type of parameters to be returned
2055: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
2056: MAT_GLOBAL_SUM - sum over all processors)
2057: - info - matrix information context
2059: Notes:
2060: The MatInfo context contains a variety of matrix data, including
2061: number of nonzeros allocated and used, number of mallocs during
2062: matrix assembly, etc. Additional information for factored matrices
2063: is provided (such as the fill ratio, number of mallocs during
2064: factorization, etc.). Much of this info is printed to PETSC_STDOUT
2065: when using the runtime options
2066: $ -info -mat_view_info
2068: Example for C/C++ Users:
2069: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
2070: data within the MatInfo context. For example,
2071: .vb
2072: MatInfo info;
2073: Mat A;
2074: double mal, nz_a, nz_u;
2076: MatGetInfo(A,MAT_LOCAL,&info);
2077: mal = info.mallocs;
2078: nz_a = info.nz_allocated;
2079: .ve
2081: Example for Fortran Users:
2082: Fortran users should declare info as a double precision
2083: array of dimension MAT_INFO_SIZE, and then extract the parameters
2084: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
2085: a complete list of parameter names.
2086: .vb
2087: double precision info(MAT_INFO_SIZE)
2088: double precision mal, nz_a
2089: Mat A
2090: integer ierr
2092: call MatGetInfo(A,MAT_LOCAL,info,ierr)
2093: mal = info(MAT_INFO_MALLOCS)
2094: nz_a = info(MAT_INFO_NZ_ALLOCATED)
2095: .ve
2097: Level: intermediate
2099: Concepts: matrices^getting information on
2100:
2101: Developer Note: fortran interface is not autogenerated as the f90
2102: interface defintion cannot be generated correctly [due to MatInfo]
2103:
2104: @*/
2105: PetscErrorCode MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
2106: {
2113: if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2114: MatPreallocated(mat);
2115: (*mat->ops->getinfo)(mat,flag,info);
2116: return(0);
2117: }
2119: /* ----------------------------------------------------------*/
2122: /*@C
2123: MatILUDTFactor - Performs a drop tolerance ILU factorization.
2125: Collective on Mat
2127: Input Parameters:
2128: + mat - the matrix
2129: . row - row permutation
2130: . col - column permutation
2131: - info - information about the factorization to be done
2133: Output Parameters:
2134: . fact - the factored matrix
2136: Level: developer
2138: Notes:
2139: Most users should employ the simplified KSP interface for linear solvers
2140: instead of working directly with matrix algebra routines such as this.
2141: See, e.g., KSPCreate().
2143: This is currently only supported for the SeqAIJ matrix format using code
2144: from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or
2145: Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
2146: and thus can be distributed with PETSc.
2148: Concepts: matrices^ILUDT factorization
2150: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2151: @*/
2152: PetscErrorCode MatILUDTFactor(Mat mat,IS row,IS col,const MatFactorInfo *info,Mat *fact)
2153: {
2163: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2164: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2165: if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2166: MatPreallocated(mat);
2167: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2168: (*mat->ops->iludtfactor)(mat,row,col,info,fact);
2169: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2170: PetscObjectStateIncrease((PetscObject)*fact);
2172: return(0);
2173: }
2177: /*@C
2178: MatLUFactor - Performs in-place LU factorization of matrix.
2180: Collective on Mat
2182: Input Parameters:
2183: + mat - the matrix
2184: . row - row permutation
2185: . col - column permutation
2186: - info - options for factorization, includes
2187: $ fill - expected fill as ratio of original fill.
2188: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2189: $ Run with the option -info to determine an optimal value to use
2191: Notes:
2192: Most users should employ the simplified KSP interface for linear solvers
2193: instead of working directly with matrix algebra routines such as this.
2194: See, e.g., KSPCreate().
2196: This changes the state of the matrix to a factored matrix; it cannot be used
2197: for example with MatSetValues() unless one first calls MatSetUnfactored().
2199: Level: developer
2201: Concepts: matrices^LU factorization
2203: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
2204: MatGetOrdering(), MatSetUnfactored(), MatFactorInfo
2206: Developer Note: fortran interface is not autogenerated as the f90
2207: interface defintion cannot be generated correctly [due to MatFactorInfo]
2209: @*/
2210: PetscErrorCode MatLUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2211: {
2220: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2221: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2222: if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2223: MatPreallocated(mat);
2225: PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
2226: (*mat->ops->lufactor)(mat,row,col,info);
2227: PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
2228: PetscObjectStateIncrease((PetscObject)mat);
2229: return(0);
2230: }
2234: /*@C
2235: MatILUFactor - Performs in-place ILU factorization of matrix.
2237: Collective on Mat
2239: Input Parameters:
2240: + mat - the matrix
2241: . row - row permutation
2242: . col - column permutation
2243: - info - structure containing
2244: $ levels - number of levels of fill.
2245: $ expected fill - as ratio of original fill.
2246: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2247: missing diagonal entries)
2249: Notes:
2250: Probably really in-place only when level of fill is zero, otherwise allocates
2251: new space to store factored matrix and deletes previous memory.
2253: Most users should employ the simplified KSP interface for linear solvers
2254: instead of working directly with matrix algebra routines such as this.
2255: See, e.g., KSPCreate().
2257: Level: developer
2259: Concepts: matrices^ILU factorization
2261: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2263: Developer Note: fortran interface is not autogenerated as the f90
2264: interface defintion cannot be generated correctly [due to MatFactorInfo]
2266: @*/
2267: PetscErrorCode MatILUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2268: {
2277: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
2278: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2279: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2280: if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2281: MatPreallocated(mat);
2283: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2284: (*mat->ops->ilufactor)(mat,row,col,info);
2285: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2286: PetscObjectStateIncrease((PetscObject)mat);
2287: return(0);
2288: }
2292: /*@C
2293: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
2294: Call this routine before calling MatLUFactorNumeric().
2296: Collective on Mat
2298: Input Parameters:
2299: + fact - the factor matrix obtained with MatGetFactor()
2300: . mat - the matrix
2301: . row, col - row and column permutations
2302: - info - options for factorization, includes
2303: $ fill - expected fill as ratio of original fill.
2304: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2305: $ Run with the option -info to determine an optimal value to use
2308: Notes:
2309: See the users manual for additional information about
2310: choosing the fill factor for better efficiency.
2312: Most users should employ the simplified KSP interface for linear solvers
2313: instead of working directly with matrix algebra routines such as this.
2314: See, e.g., KSPCreate().
2316: Level: developer
2318: Concepts: matrices^LU symbolic factorization
2320: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2322: Developer Note: fortran interface is not autogenerated as the f90
2323: interface defintion cannot be generated correctly [due to MatFactorInfo]
2325: @*/
2326: PetscErrorCode MatLUFactorSymbolic(Mat fact,Mat mat,IS row,IS col,const MatFactorInfo *info)
2327: {
2337: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2338: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2339: if (!(fact)->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",((PetscObject)mat)->type_name);
2340: MatPreallocated(mat);
2342: PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
2343: (fact->ops->lufactorsymbolic)(fact,mat,row,col,info);
2344: PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
2345: PetscObjectStateIncrease((PetscObject)fact);
2346: return(0);
2347: }
2351: /*@C
2352: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
2353: Call this routine after first calling MatLUFactorSymbolic().
2355: Collective on Mat
2357: Input Parameters:
2358: + fact - the factor matrix obtained with MatGetFactor()
2359: . mat - the matrix
2360: - info - options for factorization
2362: Notes:
2363: See MatLUFactor() for in-place factorization. See
2364: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
2366: Most users should employ the simplified KSP interface for linear solvers
2367: instead of working directly with matrix algebra routines such as this.
2368: See, e.g., KSPCreate().
2370: Level: developer
2372: Concepts: matrices^LU numeric factorization
2374: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
2376: Developer Note: fortran interface is not autogenerated as the f90
2377: interface defintion cannot be generated correctly [due to MatFactorInfo]
2379: @*/
2380: PetscErrorCode MatLUFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2381: {
2389: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2390: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2391: SETERRQ4(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);
2392: }
2393: if (!(fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2394: MatPreallocated(mat);
2395: PetscLogEventBegin(MAT_LUFactorNumeric,mat,fact,0,0);
2396: (fact->ops->lufactornumeric)(fact,mat,info);
2397: PetscLogEventEnd(MAT_LUFactorNumeric,mat,fact,0,0);
2399: MatView_Private(fact);
2400: PetscObjectStateIncrease((PetscObject)fact);
2401: return(0);
2402: }
2406: /*@C
2407: MatCholeskyFactor - Performs in-place Cholesky factorization of a
2408: symmetric matrix.
2410: Collective on Mat
2412: Input Parameters:
2413: + mat - the matrix
2414: . perm - row and column permutations
2415: - f - expected fill as ratio of original fill
2417: Notes:
2418: See MatLUFactor() for the nonsymmetric case. See also
2419: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
2421: Most users should employ the simplified KSP interface for linear solvers
2422: instead of working directly with matrix algebra routines such as this.
2423: See, e.g., KSPCreate().
2425: Level: developer
2427: Concepts: matrices^Cholesky factorization
2429: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
2430: MatGetOrdering()
2432: Developer Note: fortran interface is not autogenerated as the f90
2433: interface defintion cannot be generated correctly [due to MatFactorInfo]
2435: @*/
2436: PetscErrorCode MatCholeskyFactor(Mat mat,IS perm,const MatFactorInfo *info)
2437: {
2445: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2446: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2447: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2448: if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2449: MatPreallocated(mat);
2451: PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
2452: (*mat->ops->choleskyfactor)(mat,perm,info);
2453: PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
2454: PetscObjectStateIncrease((PetscObject)mat);
2455: return(0);
2456: }
2460: /*@C
2461: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
2462: of a symmetric matrix.
2464: Collective on Mat
2466: Input Parameters:
2467: + fact - the factor matrix obtained with MatGetFactor()
2468: . mat - the matrix
2469: . perm - row and column permutations
2470: - info - options for factorization, includes
2471: $ fill - expected fill as ratio of original fill.
2472: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2473: $ Run with the option -info to determine an optimal value to use
2475: Notes:
2476: See MatLUFactorSymbolic() for the nonsymmetric case. See also
2477: MatCholeskyFactor() and MatCholeskyFactorNumeric().
2479: Most users should employ the simplified KSP interface for linear solvers
2480: instead of working directly with matrix algebra routines such as this.
2481: See, e.g., KSPCreate().
2483: Level: developer
2485: Concepts: matrices^Cholesky symbolic factorization
2487: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
2488: MatGetOrdering()
2490: Developer Note: fortran interface is not autogenerated as the f90
2491: interface defintion cannot be generated correctly [due to MatFactorInfo]
2493: @*/
2494: PetscErrorCode MatCholeskyFactorSymbolic(Mat fact,Mat mat,IS perm,const MatFactorInfo *info)
2495: {
2504: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2505: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2506: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2507: if (!(fact)->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2508: MatPreallocated(mat);
2510: PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2511: (fact->ops->choleskyfactorsymbolic)(fact,mat,perm,info);
2512: PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2513: PetscObjectStateIncrease((PetscObject)fact);
2514: return(0);
2515: }
2519: /*@C
2520: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
2521: of a symmetric matrix. Call this routine after first calling
2522: MatCholeskyFactorSymbolic().
2524: Collective on Mat
2526: Input Parameters:
2527: + fact - the factor matrix obtained with MatGetFactor()
2528: . mat - the initial matrix
2529: . info - options for factorization
2530: - fact - the symbolic factor of mat
2533: Notes:
2534: Most users should employ the simplified KSP interface for linear solvers
2535: instead of working directly with matrix algebra routines such as this.
2536: See, e.g., KSPCreate().
2538: Level: developer
2540: Concepts: matrices^Cholesky numeric factorization
2542: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
2544: Developer Note: fortran interface is not autogenerated as the f90
2545: interface defintion cannot be generated correctly [due to MatFactorInfo]
2547: @*/
2548: PetscErrorCode MatCholeskyFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2549: {
2557: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2558: if (!(fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2559: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2560: SETERRQ4(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);
2561: }
2562: MatPreallocated(mat);
2564: PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,fact,0,0);
2565: (fact->ops->choleskyfactornumeric)(fact,mat,info);
2566: PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,fact,0,0);
2568: MatView_Private(fact);
2569: PetscObjectStateIncrease((PetscObject)fact);
2570: return(0);
2571: }
2573: /* ----------------------------------------------------------------*/
2576: /*@
2577: MatSolve - Solves A x = b, given a factored matrix.
2579: Collective on Mat and Vec
2581: Input Parameters:
2582: + mat - the factored matrix
2583: - b - the right-hand-side vector
2585: Output Parameter:
2586: . x - the result vector
2588: Notes:
2589: The vectors b and x cannot be the same. I.e., one cannot
2590: call MatSolve(A,x,x).
2592: Notes:
2593: Most users should employ the simplified KSP interface for linear solvers
2594: instead of working directly with matrix algebra routines such as this.
2595: See, e.g., KSPCreate().
2597: Level: developer
2599: Concepts: matrices^triangular solves
2601: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
2602: @*/
2603: PetscErrorCode MatSolve(Mat mat,Vec b,Vec x)
2604: {
2614: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2615: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2616: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2617: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
2618: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
2619: if (!mat->rmap->N && !mat->cmap->N) return(0);
2620: if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2621: MatPreallocated(mat);
2623: PetscLogEventBegin(MAT_Solve,mat,b,x,0);
2624: (*mat->ops->solve)(mat,b,x);
2625: PetscLogEventEnd(MAT_Solve,mat,b,x,0);
2626: PetscObjectStateIncrease((PetscObject)x);
2627: return(0);
2628: }
2632: PetscErrorCode MatMatSolve_Basic(Mat A,Mat B,Mat X)
2633: {
2635: Vec b,x;
2636: PetscInt m,N,i;
2637: PetscScalar *bb,*xx;
2640: MatGetArray(B,&bb);
2641: MatGetArray(X,&xx);
2642: MatGetLocalSize(B,&m,PETSC_NULL); /* number local rows */
2643: MatGetSize(B,PETSC_NULL,&N); /* total columns in dense matrix */
2644: VecCreateMPIWithArray(((PetscObject)A)->comm,m,PETSC_DETERMINE,PETSC_NULL,&b);
2645: VecCreateMPIWithArray(((PetscObject)A)->comm,m,PETSC_DETERMINE,PETSC_NULL,&x);
2646: for (i=0; i<N; i++) {
2647: VecPlaceArray(b,bb + i*m);
2648: VecPlaceArray(x,xx + i*m);
2649: MatSolve(A,b,x);
2650: VecResetArray(x);
2651: VecResetArray(b);
2652: }
2653: VecDestroy(b);
2654: VecDestroy(x);
2655: MatRestoreArray(B,&bb);
2656: MatRestoreArray(X,&xx);
2657: return(0);
2658: }
2662: /*@
2663: MatMatSolve - Solves A X = B, given a factored matrix.
2665: Collective on Mat
2667: Input Parameters:
2668: + mat - the factored matrix
2669: - B - the right-hand-side matrix (dense matrix)
2671: Output Parameter:
2672: . X - the result matrix (dense matrix)
2674: Notes:
2675: The matrices b and x cannot be the same. I.e., one cannot
2676: call MatMatSolve(A,x,x).
2678: Notes:
2679: Most users should usually employ the simplified KSP interface for linear solvers
2680: instead of working directly with matrix algebra routines such as this.
2681: See, e.g., KSPCreate(). However KSP can only solve for one vector (column of X)
2682: at a time.
2684: Level: developer
2686: Concepts: matrices^triangular solves
2688: .seealso: MatMatSolveAdd(), MatMatSolveTranspose(), MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()
2689: @*/
2690: PetscErrorCode MatMatSolve(Mat A,Mat B,Mat X)
2691: {
2701: if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different matrices");
2702: if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2703: if (A->cmap->N != X->rmap->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat X: global dim %D %D",A->cmap->N,X->rmap->N);
2704: if (A->rmap->N != B->rmap->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %D %D",A->rmap->N,B->rmap->N);
2705: if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D",A->rmap->n,B->rmap->n);
2706: if (!A->rmap->N && !A->cmap->N) return(0);
2707: MatPreallocated(A);
2709: PetscLogEventBegin(MAT_MatSolve,A,B,X,0);
2710: if (!A->ops->matsolve) {
2711: PetscInfo1(A,"Mat type %s using basic MatMatSolve",((PetscObject)A)->type_name);
2712: MatMatSolve_Basic(A,B,X);
2713: } else {
2714: (*A->ops->matsolve)(A,B,X);
2715: }
2716: PetscLogEventEnd(MAT_MatSolve,A,B,X,0);
2717: PetscObjectStateIncrease((PetscObject)X);
2718: return(0);
2719: }
2724: /*@
2725: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU, or
2726: U^T*D^(1/2) x = b, given a factored symmetric matrix, A = U^T*D*U,
2728: Collective on Mat and Vec
2730: Input Parameters:
2731: + mat - the factored matrix
2732: - b - the right-hand-side vector
2734: Output Parameter:
2735: . x - the result vector
2737: Notes:
2738: MatSolve() should be used for most applications, as it performs
2739: a forward solve followed by a backward solve.
2741: The vectors b and x cannot be the same, i.e., one cannot
2742: call MatForwardSolve(A,x,x).
2744: For matrix in seqsbaij format with block size larger than 1,
2745: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2746: MatForwardSolve() solves U^T*D y = b, and
2747: MatBackwardSolve() solves U x = y.
2748: Thus they do not provide a symmetric preconditioner.
2750: Most users should employ the simplified KSP interface for linear solvers
2751: instead of working directly with matrix algebra routines such as this.
2752: See, e.g., KSPCreate().
2754: Level: developer
2756: Concepts: matrices^forward solves
2758: .seealso: MatSolve(), MatBackwardSolve()
2759: @*/
2760: PetscErrorCode MatForwardSolve(Mat mat,Vec b,Vec x)
2761: {
2771: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2772: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2773: if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2774: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2775: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
2776: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
2777: MatPreallocated(mat);
2778: PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
2779: (*mat->ops->forwardsolve)(mat,b,x);
2780: PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
2781: PetscObjectStateIncrease((PetscObject)x);
2782: return(0);
2783: }
2787: /*@
2788: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
2789: D^(1/2) U x = b, given a factored symmetric matrix, A = U^T*D*U,
2791: Collective on Mat and Vec
2793: Input Parameters:
2794: + mat - the factored matrix
2795: - b - the right-hand-side vector
2797: Output Parameter:
2798: . x - the result vector
2800: Notes:
2801: MatSolve() should be used for most applications, as it performs
2802: a forward solve followed by a backward solve.
2804: The vectors b and x cannot be the same. I.e., one cannot
2805: call MatBackwardSolve(A,x,x).
2807: For matrix in seqsbaij format with block size larger than 1,
2808: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2809: MatForwardSolve() solves U^T*D y = b, and
2810: MatBackwardSolve() solves U x = y.
2811: Thus they do not provide a symmetric preconditioner.
2813: Most users should employ the simplified KSP interface for linear solvers
2814: instead of working directly with matrix algebra routines such as this.
2815: See, e.g., KSPCreate().
2817: Level: developer
2819: Concepts: matrices^backward solves
2821: .seealso: MatSolve(), MatForwardSolve()
2822: @*/
2823: PetscErrorCode MatBackwardSolve(Mat mat,Vec b,Vec x)
2824: {
2834: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2835: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2836: if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2837: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2838: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
2839: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
2840: MatPreallocated(mat);
2842: PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
2843: (*mat->ops->backwardsolve)(mat,b,x);
2844: PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
2845: PetscObjectStateIncrease((PetscObject)x);
2846: return(0);
2847: }
2851: /*@
2852: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
2854: Collective on Mat and Vec
2856: Input Parameters:
2857: + mat - the factored matrix
2858: . b - the right-hand-side vector
2859: - y - the vector to be added to
2861: Output Parameter:
2862: . x - the result vector
2864: Notes:
2865: The vectors b and x cannot be the same. I.e., one cannot
2866: call MatSolveAdd(A,x,y,x).
2868: Most users should employ the simplified KSP interface for linear solvers
2869: instead of working directly with matrix algebra routines such as this.
2870: See, e.g., KSPCreate().
2872: Level: developer
2874: Concepts: matrices^triangular solves
2876: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
2877: @*/
2878: PetscErrorCode MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
2879: {
2880: PetscScalar one = 1.0;
2881: Vec tmp;
2893: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2894: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2895: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
2896: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
2897: if (mat->rmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->rmap->N,y->map->N);
2898: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
2899: if (x->map->n != y->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->map->n,y->map->n);
2900: MatPreallocated(mat);
2902: PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
2903: if (mat->ops->solveadd) {
2904: (*mat->ops->solveadd)(mat,b,y,x);
2905: } else {
2906: /* do the solve then the add manually */
2907: if (x != y) {
2908: MatSolve(mat,b,x);
2909: VecAXPY(x,one,y);
2910: } else {
2911: VecDuplicate(x,&tmp);
2912: PetscLogObjectParent(mat,tmp);
2913: VecCopy(x,tmp);
2914: MatSolve(mat,b,x);
2915: VecAXPY(x,one,tmp);
2916: VecDestroy(tmp);
2917: }
2918: }
2919: PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
2920: PetscObjectStateIncrease((PetscObject)x);
2921: return(0);
2922: }
2926: /*@
2927: MatSolveTranspose - Solves A' x = b, given a factored matrix.
2929: Collective on Mat and Vec
2931: Input Parameters:
2932: + mat - the factored matrix
2933: - b - the right-hand-side vector
2935: Output Parameter:
2936: . x - the result vector
2938: Notes:
2939: The vectors b and x cannot be the same. I.e., one cannot
2940: call MatSolveTranspose(A,x,x).
2942: Most users should employ the simplified KSP interface for linear solvers
2943: instead of working directly with matrix algebra routines such as this.
2944: See, e.g., KSPCreate().
2946: Level: developer
2948: Concepts: matrices^triangular solves
2950: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
2951: @*/
2952: PetscErrorCode MatSolveTranspose(Mat mat,Vec b,Vec x)
2953: {
2963: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2964: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2965: if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",((PetscObject)mat)->type_name);
2966: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
2967: if (mat->cmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->cmap->N,b->map->N);
2968: MatPreallocated(mat);
2969: PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
2970: (*mat->ops->solvetranspose)(mat,b,x);
2971: PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
2972: PetscObjectStateIncrease((PetscObject)x);
2973: return(0);
2974: }
2978: /*@
2979: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
2980: factored matrix.
2982: Collective on Mat and Vec
2984: Input Parameters:
2985: + mat - the factored matrix
2986: . b - the right-hand-side vector
2987: - y - the vector to be added to
2989: Output Parameter:
2990: . x - the result vector
2992: Notes:
2993: The vectors b and x cannot be the same. I.e., one cannot
2994: call MatSolveTransposeAdd(A,x,y,x).
2996: Most users should employ the simplified KSP interface for linear solvers
2997: instead of working directly with matrix algebra routines such as this.
2998: See, e.g., KSPCreate().
3000: Level: developer
3002: Concepts: matrices^triangular solves
3004: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
3005: @*/
3006: PetscErrorCode MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
3007: {
3008: PetscScalar one = 1.0;
3010: Vec tmp;
3021: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
3022: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
3023: if (mat->rmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->rmap->N,x->map->N);
3024: if (mat->cmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->cmap->N,b->map->N);
3025: if (mat->cmap->N != y->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->cmap->N,y->map->N);
3026: if (x->map->n != y->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->map->n,y->map->n);
3027: MatPreallocated(mat);
3029: PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
3030: if (mat->ops->solvetransposeadd) {
3031: (*mat->ops->solvetransposeadd)(mat,b,y,x);
3032: } else {
3033: /* do the solve then the add manually */
3034: if (x != y) {
3035: MatSolveTranspose(mat,b,x);
3036: VecAXPY(x,one,y);
3037: } else {
3038: VecDuplicate(x,&tmp);
3039: PetscLogObjectParent(mat,tmp);
3040: VecCopy(x,tmp);
3041: MatSolveTranspose(mat,b,x);
3042: VecAXPY(x,one,tmp);
3043: VecDestroy(tmp);
3044: }
3045: }
3046: PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
3047: PetscObjectStateIncrease((PetscObject)x);
3048: return(0);
3049: }
3050: /* ----------------------------------------------------------------*/
3054: /*@
3055: MatRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3057: Collective on Mat and Vec
3059: Input Parameters:
3060: + mat - the matrix
3061: . b - the right hand side
3062: . omega - the relaxation factor
3063: . flag - flag indicating the type of SOR (see below)
3064: . shift - diagonal shift
3065: . its - the number of iterations
3066: - lits - the number of local iterations
3068: Output Parameters:
3069: . x - the solution (can contain an initial guess)
3071: SOR Flags:
3072: . SOR_FORWARD_SWEEP - forward SOR
3073: . SOR_BACKWARD_SWEEP - backward SOR
3074: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
3075: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
3076: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
3077: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
3078: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
3079: upper/lower triangular part of matrix to
3080: vector (with omega)
3081: . SOR_ZERO_INITIAL_GUESS - zero initial guess
3083: Notes:
3084: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
3085: SOR_LOCAL_SYMMETRIC_SWEEP perform separate independent smoothings
3086: on each processor.
3088: Application programmers will not generally use MatRelax() directly,
3089: but instead will employ the KSP/PC interface.
3091: Notes for Advanced Users:
3092: The flags are implemented as bitwise inclusive or operations.
3093: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
3094: to specify a zero initial guess for SSOR.
3096: Most users should employ the simplified KSP interface for linear solvers
3097: instead of working directly with matrix algebra routines such as this.
3098: See, e.g., KSPCreate().
3100: See also, MatPBRelax(). This routine will automatically call the point block
3101: version if the point version is not available.
3103: Level: developer
3105: Concepts: matrices^relaxation
3106: Concepts: matrices^SOR
3107: Concepts: matrices^Gauss-Seidel
3109: @*/
3110: PetscErrorCode MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3111: {
3121: if (!mat->ops->relax && !mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3122: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3123: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3124: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3125: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3126: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3127: if (its <= 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
3128: if (lits <= 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Relaxation requires local its %D positive",lits);
3130: MatPreallocated(mat);
3131: PetscLogEventBegin(MAT_Relax,mat,b,x,0);
3132: if (mat->ops->relax) {
3133: ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
3134: } else {
3135: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
3136: }
3137: PetscLogEventEnd(MAT_Relax,mat,b,x,0);
3138: PetscObjectStateIncrease((PetscObject)x);
3139: return(0);
3140: }
3144: /*@
3145: MatPBRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3147: Collective on Mat and Vec
3149: See MatRelax() for usage
3151: For multi-component PDEs where the Jacobian is stored in a point block format
3152: (with the PETSc BAIJ matrix formats) the relaxation is done one point block at
3153: a time. That is, the small (for example, 4 by 4) blocks along the diagonal are solved
3154: simultaneously (that is a 4 by 4 linear solve is done) to update all the values at a point.
3156: Level: developer
3158: @*/
3159: PetscErrorCode MatPBRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3160: {
3170: if (!mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3171: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3172: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3173: if (mat->cmap->N != x->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->cmap->N,x->map->N);
3174: if (mat->rmap->N != b->map->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->rmap->N,b->map->N);
3175: if (mat->rmap->n != b->map->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->rmap->n,b->map->n);
3176: MatPreallocated(mat);
3178: PetscLogEventBegin(MAT_Relax,mat,b,x,0);
3179: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
3180: PetscLogEventEnd(MAT_Relax,mat,b,x,0);
3181: PetscObjectStateIncrease((PetscObject)x);
3182: return(0);
3183: }
3187: /*
3188: Default matrix copy routine.
3189: */
3190: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
3191: {
3192: PetscErrorCode ierr;
3193: PetscInt i,rstart,rend,nz;
3194: const PetscInt *cwork;
3195: const PetscScalar *vwork;
3198: if (B->assembled) {
3199: MatZeroEntries(B);
3200: }
3201: MatGetOwnershipRange(A,&rstart,&rend);
3202: for (i=rstart; i<rend; i++) {
3203: MatGetRow(A,i,&nz,&cwork,&vwork);
3204: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
3205: MatRestoreRow(A,i,&nz,&cwork,&vwork);
3206: }
3207: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3208: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3209: PetscObjectStateIncrease((PetscObject)B);
3210: return(0);
3211: }
3215: /*@
3216: MatCopy - Copys a matrix to another matrix.
3218: Collective on Mat
3220: Input Parameters:
3221: + A - the matrix
3222: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
3224: Output Parameter:
3225: . B - where the copy is put
3227: Notes:
3228: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
3229: same nonzero pattern or the routine will crash.
3231: MatCopy() copies the matrix entries of a matrix to another existing
3232: matrix (after first zeroing the second matrix). A related routine is
3233: MatConvert(), which first creates a new matrix and then copies the data.
3235: Level: intermediate
3236:
3237: Concepts: matrices^copying
3239: .seealso: MatConvert(), MatDuplicate()
3241: @*/
3242: PetscErrorCode MatCopy(Mat A,Mat B,MatStructure str)
3243: {
3245: PetscInt i;
3253: MatPreallocated(B);
3254: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3255: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3256: if (A->rmap->N != B->rmap->N || A->cmap->N != B->cmap->N) SETERRQ4(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);
3257: MatPreallocated(A);
3259: PetscLogEventBegin(MAT_Copy,A,B,0,0);
3260: if (A->ops->copy) {
3261: (*A->ops->copy)(A,B,str);
3262: } else { /* generic conversion */
3263: MatCopy_Basic(A,B,str);
3264: }
3265: if (A->mapping) {
3266: if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
3267: MatSetLocalToGlobalMapping(B,A->mapping);
3268: }
3269: if (A->bmapping) {
3270: if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
3271: MatSetLocalToGlobalMappingBlock(B,A->mapping);
3272: }
3274: B->stencil.dim = A->stencil.dim;
3275: B->stencil.noc = A->stencil.noc;
3276: for (i=0; i<=A->stencil.dim; i++) {
3277: B->stencil.dims[i] = A->stencil.dims[i];
3278: B->stencil.starts[i] = A->stencil.starts[i];
3279: }
3281: PetscLogEventEnd(MAT_Copy,A,B,0,0);
3282: PetscObjectStateIncrease((PetscObject)B);
3283: return(0);
3284: }
3288: /*@C
3289: MatConvert - Converts a matrix to another matrix, either of the same
3290: or different type.
3292: Collective on Mat
3294: Input Parameters:
3295: + mat - the matrix
3296: . newtype - new matrix type. Use MATSAME to create a new matrix of the
3297: same type as the original matrix.
3298: - reuse - denotes if the destination matrix is to be created or reused. Currently
3299: MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
3300: MAT_INITIAL_MATRIX.
3302: Output Parameter:
3303: . M - pointer to place new matrix
3305: Notes:
3306: MatConvert() first creates a new matrix and then copies the data from
3307: the first matrix. A related routine is MatCopy(), which copies the matrix
3308: entries of one matrix to another already existing matrix context.
3310: Cannot be used to convert a sequential matrix to parallel or parallel to sequential,
3311: the MPI communicator of the generated matrix is always the same as the communicator
3312: of the input matrix.
3314: Level: intermediate
3316: Concepts: matrices^converting between storage formats
3318: .seealso: MatCopy(), MatDuplicate()
3319: @*/
3320: PetscErrorCode MatConvert(Mat mat, const MatType newtype,MatReuse reuse,Mat *M)
3321: {
3322: PetscErrorCode ierr;
3323: PetscTruth sametype,issame,flg;
3324: char convname[256],mtype[256];
3325: Mat B;
3331: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3332: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3333: MatPreallocated(mat);
3335: PetscOptionsGetString(((PetscObject)mat)->prefix,"-matconvert_type",mtype,256,&flg);
3336: if (flg) {
3337: newtype = mtype;
3338: }
3339: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
3340: PetscStrcmp(newtype,"same",&issame);
3341: if ((reuse == MAT_REUSE_MATRIX) && (mat != *M)) {
3342: SETERRQ(PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for in-place conversion currently");
3343: }
3345: if ((reuse == MAT_REUSE_MATRIX) && (issame || sametype)) return(0);
3346:
3347: if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
3348: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
3349: } else {
3350: PetscErrorCode (*conv)(Mat, const MatType,MatReuse,Mat*)=PETSC_NULL;
3351: const char *prefix[3] = {"seq","mpi",""};
3352: PetscInt i;
3353: /*
3354: Order of precedence:
3355: 1) See if a specialized converter is known to the current matrix.
3356: 2) See if a specialized converter is known to the desired matrix class.
3357: 3) See if a good general converter is registered for the desired class
3358: (as of 6/27/03 only MATMPIADJ falls into this category).
3359: 4) See if a good general converter is known for the current matrix.
3360: 5) Use a really basic converter.
3361: */
3362:
3363: /* 1) See if a specialized converter is known to the current matrix and the desired class */
3364: for (i=0; i<3; i++) {
3365: PetscStrcpy(convname,"MatConvert_");
3366: PetscStrcat(convname,((PetscObject)mat)->type_name);
3367: PetscStrcat(convname,"_");
3368: PetscStrcat(convname,prefix[i]);
3369: PetscStrcat(convname,newtype);
3370: PetscStrcat(convname,"_C");
3371: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3372: if (conv) goto foundconv;
3373: }
3375: /* 2) See if a specialized converter is known to the desired matrix class. */
3376: MatCreate(((PetscObject)mat)->comm,&B);
3377: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N);
3378: MatSetType(B,newtype);
3379: for (i=0; i<3; i++) {
3380: PetscStrcpy(convname,"MatConvert_");
3381: PetscStrcat(convname,((PetscObject)mat)->type_name);
3382: PetscStrcat(convname,"_");
3383: PetscStrcat(convname,prefix[i]);
3384: PetscStrcat(convname,newtype);
3385: PetscStrcat(convname,"_C");
3386: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
3387: if (conv) {
3388: MatDestroy(B);
3389: goto foundconv;
3390: }
3391: }
3393: /* 3) See if a good general converter is registered for the desired class */
3394: conv = B->ops->convertfrom;
3395: MatDestroy(B);
3396: if (conv) goto foundconv;
3398: /* 4) See if a good general converter is known for the current matrix */
3399: if (mat->ops->convert) {
3400: conv = mat->ops->convert;
3401: }
3402: if (conv) goto foundconv;
3404: /* 5) Use a really basic converter. */
3405: conv = MatConvert_Basic;
3407: foundconv:
3408: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3409: (*conv)(mat,newtype,reuse,M);
3410: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3411: }
3412: PetscObjectStateIncrease((PetscObject)*M);
3413: return(0);
3414: }
3418: /*@C
3419: MatFactorGetSolverPackage - Returns name of the package providing the factorization routines
3421: Not Collective
3423: Input Parameter:
3424: . mat - the matrix, must be a factored matrix
3426: Output Parameter:
3427: . type - the string name of the package (do not free this string)
3429: Notes:
3430: In Fortran you pass in a empty string and the package name will be copied into it.
3431: (Make sure the string is long enough)
3433: Level: intermediate
3435: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable(), MatGetFactor()
3436: @*/
3437: PetscErrorCode MatFactorGetSolverPackage(Mat mat, const MatSolverPackage *type)
3438: {
3439: PetscErrorCode ierr;
3440: PetscErrorCode (*conv)(Mat,const MatSolverPackage*);
3445: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Only for factored matrix");
3446: PetscObjectQueryFunction((PetscObject)mat,"MatFactorGetSolverPackage_C",(void (**)(void))&conv);
3447: if (!conv) {
3448: *type = MAT_SOLVER_PETSC;
3449: } else {
3450: (*conv)(mat,type);
3451: }
3452: return(0);
3453: }
3457: /*@C
3458: MatGetFactor - Returns a matrix suitable to calls to MatXXFactorSymbolic()
3460: Collective on Mat
3462: Input Parameters:
3463: + mat - the matrix
3464: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3465: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3467: Output Parameters:
3468: . f - the factor matrix used with MatXXFactorSymbolic() calls
3470: Notes:
3471: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3472: such as pastix, superlu, mumps, spooles etc.
3474: PETSc must have been config/configure.py to use the external solver, using the option --download-package
3476: Level: intermediate
3478: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable()
3479: @*/
3480: PetscErrorCode MatGetFactor(Mat mat, const MatSolverPackage type,MatFactorType ftype,Mat *f)
3481: {
3482: PetscErrorCode ierr;
3483: char convname[256];
3484: PetscErrorCode (*conv)(Mat,MatFactorType,Mat*);
3490: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3491: MatPreallocated(mat);
3493: PetscStrcpy(convname,"MatGetFactor_");
3494: PetscStrcat(convname,((PetscObject)mat)->type_name);
3495: PetscStrcat(convname,"_");
3496: PetscStrcat(convname,type);
3497: PetscStrcat(convname,"_C");
3498: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3499: if (!conv) {
3500: PetscTruth flag;
3501: PetscStrcasecmp(MAT_SOLVER_PETSC,type,&flag);
3502: if (flag) {
3503: SETERRQ1(PETSC_ERR_SUP,"Matrix format %s does not have a built-in PETSc direct solver",((PetscObject)mat)->type_name);
3504: } else {
3505: SETERRQ3(PETSC_ERR_SUP,"Matrix format %s does not have a solver %s. Perhaps you must config/configure.py with --download-%s",((PetscObject)mat)->type_name,type,type);
3506: }
3507: }
3508: (*conv)(mat,ftype,f);
3509: return(0);
3510: }
3514: /*@C
3515: MatGetFactorAvailable - Returns a a flag if matrix supports particular package and factor type
3517: Collective on Mat
3519: Input Parameters:
3520: + mat - the matrix
3521: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3522: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3524: Output Parameter:
3525: . flg - PETSC_TRUE if the factorization is available
3527: Notes:
3528: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3529: such as pastix, superlu, mumps, spooles etc.
3531: PETSc must have been config/configure.py to use the external solver, using the option --download-package
3533: Level: intermediate
3535: .seealso: MatCopy(), MatDuplicate(), MatGetFactor()
3536: @*/
3537: PetscErrorCode MatGetFactorAvailable(Mat mat, const MatSolverPackage type,MatFactorType ftype,PetscTruth *flg)
3538: {
3539: PetscErrorCode ierr;
3540: char convname[256];
3541: PetscErrorCode (*conv)(Mat,MatFactorType,PetscTruth*);
3547: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3548: MatPreallocated(mat);
3550: PetscStrcpy(convname,"MatGetFactorAvailable_");
3551: PetscStrcat(convname,((PetscObject)mat)->type_name);
3552: PetscStrcat(convname,"_");
3553: PetscStrcat(convname,type);
3554: PetscStrcat(convname,"_C");
3555: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3556: if (!conv) {
3557: *flg = PETSC_FALSE;
3558: } else {
3559: (*conv)(mat,ftype,flg);
3560: }
3561: return(0);
3562: }
3567: /*@
3568: MatDuplicate - Duplicates a matrix including the non-zero structure.
3570: Collective on Mat
3572: Input Parameters:
3573: + mat - the matrix
3574: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
3575: values as well or not
3577: Output Parameter:
3578: . M - pointer to place new matrix
3580: Level: intermediate
3582: Concepts: matrices^duplicating
3584: .seealso: MatCopy(), MatConvert()
3585: @*/
3586: PetscErrorCode MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
3587: {
3589: Mat B;
3590: PetscInt i;
3596: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3597: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3598: MatPreallocated(mat);
3600: *M = 0;
3601: if (!mat->ops->duplicate) {
3602: SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
3603: }
3604: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3605: (*mat->ops->duplicate)(mat,op,M);
3606: B = *M;
3607: if (mat->mapping) {
3608: MatSetLocalToGlobalMapping(B,mat->mapping);
3609: }
3610: if (mat->bmapping) {
3611: MatSetLocalToGlobalMappingBlock(B,mat->bmapping);
3612: }
3613: PetscMapCopy(((PetscObject)mat)->comm,mat->rmap,B->rmap);
3614: PetscMapCopy(((PetscObject)mat)->comm,mat->cmap,B->cmap);
3615:
3616: B->stencil.dim = mat->stencil.dim;
3617: B->stencil.noc = mat->stencil.noc;
3618: for (i=0; i<=mat->stencil.dim; i++) {
3619: B->stencil.dims[i] = mat->stencil.dims[i];
3620: B->stencil.starts[i] = mat->stencil.starts[i];
3621: }
3623: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3624: PetscObjectStateIncrease((PetscObject)B);
3625: return(0);
3626: }
3630: /*@
3631: MatGetDiagonal - Gets the diagonal of a matrix.
3633: Collective on Mat and Vec
3635: Input Parameters:
3636: + mat - the matrix
3637: - v - the vector for storing the diagonal
3639: Output Parameter:
3640: . v - the diagonal of the matrix
3642: Level: intermediate
3644: Concepts: matrices^accessing diagonals
3646: .seealso: MatGetRow(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs()
3647: @*/
3648: PetscErrorCode MatGetDiagonal(Mat mat,Vec v)
3649: {
3656: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3657: if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3658: MatPreallocated(mat);
3660: (*mat->ops->getdiagonal)(mat,v);
3661: PetscObjectStateIncrease((PetscObject)v);
3662: return(0);
3663: }
3667: /*@
3668: MatGetRowMin - Gets the minimum value (of the real part) of each
3669: row of the matrix
3671: Collective on Mat and Vec
3673: Input Parameters:
3674: . mat - the matrix
3676: Output Parameter:
3677: + v - the vector for storing the maximums
3678: - idx - the indices of the column found for each row (optional)
3680: Level: intermediate
3682: Notes: The result of this call are the same as if one converted the matrix to dense format
3683: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
3685: This code is only implemented for a couple of matrix formats.
3687: Concepts: matrices^getting row maximums
3689: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(),
3690: MatGetRowMax()
3691: @*/
3692: PetscErrorCode MatGetRowMin(Mat mat,Vec v,PetscInt idx[])
3693: {
3700: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3701: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3702: MatPreallocated(mat);
3704: (*mat->ops->getrowmin)(mat,v,idx);
3705: PetscObjectStateIncrease((PetscObject)v);
3706: return(0);
3707: }
3711: /*@
3712: MatGetRowMinAbs - Gets the minimum value (in absolute value) of each
3713: row of the matrix
3715: Collective on Mat and Vec
3717: Input Parameters:
3718: . mat - the matrix
3720: Output Parameter:
3721: + v - the vector for storing the minimums
3722: - idx - the indices of the column found for each row (optional)
3724: Level: intermediate
3726: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
3727: row is 0 (the first column).
3729: This code is only implemented for a couple of matrix formats.
3731: Concepts: matrices^getting row maximums
3733: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMaxAbs(), MatGetRowMin()
3734: @*/
3735: PetscErrorCode MatGetRowMinAbs(Mat mat,Vec v,PetscInt idx[])
3736: {
3743: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3744: if (!mat->ops->getrowminabs) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3745: MatPreallocated(mat);
3746: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
3748: (*mat->ops->getrowminabs)(mat,v,idx);
3749: PetscObjectStateIncrease((PetscObject)v);
3750: return(0);
3751: }
3755: /*@
3756: MatGetRowMax - Gets the maximum value (of the real part) of each
3757: row of the matrix
3759: Collective on Mat and Vec
3761: Input Parameters:
3762: . mat - the matrix
3764: Output Parameter:
3765: + v - the vector for storing the maximums
3766: - idx - the indices of the column found for each row (optional)
3768: Level: intermediate
3770: Notes: The result of this call are the same as if one converted the matrix to dense format
3771: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
3773: This code is only implemented for a couple of matrix formats.
3775: Concepts: matrices^getting row maximums
3777: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(), MatGetRowMin()
3778: @*/
3779: PetscErrorCode MatGetRowMax(Mat mat,Vec v,PetscInt idx[])
3780: {
3787: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3788: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3789: MatPreallocated(mat);
3791: (*mat->ops->getrowmax)(mat,v,idx);
3792: PetscObjectStateIncrease((PetscObject)v);
3793: return(0);
3794: }
3798: /*@
3799: MatGetRowMaxAbs - Gets the maximum value (in absolute value) of each
3800: row of the matrix
3802: Collective on Mat and Vec
3804: Input Parameters:
3805: . mat - the matrix
3807: Output Parameter:
3808: + v - the vector for storing the maximums
3809: - idx - the indices of the column found for each row (optional)
3811: Level: intermediate
3813: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
3814: row is 0 (the first column).
3816: This code is only implemented for a couple of matrix formats.
3818: Concepts: matrices^getting row maximums
3820: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
3821: @*/
3822: PetscErrorCode MatGetRowMaxAbs(Mat mat,Vec v,PetscInt idx[])
3823: {
3830: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3831: if (!mat->ops->getrowmaxabs) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3832: MatPreallocated(mat);
3833: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
3835: (*mat->ops->getrowmaxabs)(mat,v,idx);
3836: PetscObjectStateIncrease((PetscObject)v);
3837: return(0);
3838: }
3842: /*@
3843: MatGetRowSum - Gets the sum of each row of the matrix
3845: Collective on Mat and Vec
3847: Input Parameters:
3848: . mat - the matrix
3850: Output Parameter:
3851: . v - the vector for storing the maximums
3853: Level: intermediate
3855: Notes: This code is slow since it is not currently specialized for different formats
3857: Concepts: matrices^getting row sums
3859: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
3860: @*/
3861: PetscErrorCode MatGetRowSum(Mat mat, Vec v)
3862: {
3863: PetscInt start, end, row;
3864: PetscScalar *array;
3871: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3872: MatPreallocated(mat);
3873: MatGetOwnershipRange(mat, &start, &end);
3874: VecGetArray(v, &array);
3875: for(row = start; row < end; ++row) {
3876: PetscInt ncols, col;
3877: const PetscInt *cols;
3878: const PetscScalar *vals;
3880: array[row - start] = 0.0;
3881: MatGetRow(mat, row, &ncols, &cols, &vals);
3882: for(col = 0; col < ncols; col++) {
3883: array[row - start] += vals[col];
3884: }
3885: MatRestoreRow(mat, row, &ncols, &cols, &vals);
3886: }
3887: VecRestoreArray(v, &array);
3888: PetscObjectStateIncrease((PetscObject) v);
3889: return(0);
3890: }
3894: /*@
3895: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
3897: Collective on Mat
3899: Input Parameter:
3900: + mat - the matrix to transpose
3901: - reuse - store the transpose matrix in the provided B
3903: Output Parameters:
3904: . B - the transpose
3906: Notes:
3907: If you pass in &mat for B the transpose will be done in place
3909: Level: intermediate
3911: Concepts: matrices^transposing
3913: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose()
3914: @*/
3915: PetscErrorCode MatTranspose(Mat mat,MatReuse reuse,Mat *B)
3916: {
3922: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3923: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3924: if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3925: MatPreallocated(mat);
3927: PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
3928: (*mat->ops->transpose)(mat,reuse,B);
3929: PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
3930: if (B) {PetscObjectStateIncrease((PetscObject)*B);}
3931: return(0);
3932: }
3936: /*@
3937: MatIsTranspose - Test whether a matrix is another one's transpose,
3938: or its own, in which case it tests symmetry.
3940: Collective on Mat
3942: Input Parameter:
3943: + A - the matrix to test
3944: - B - the matrix to test against, this can equal the first parameter
3946: Output Parameters:
3947: . flg - the result
3949: Notes:
3950: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
3951: has a running time of the order of the number of nonzeros; the parallel
3952: test involves parallel copies of the block-offdiagonal parts of the matrix.
3954: Level: intermediate
3956: Concepts: matrices^transposing, matrix^symmetry
3958: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
3959: @*/
3960: PetscErrorCode MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
3961: {
3962: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);
3968: PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
3969: PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
3970: if (f && g) {
3971: if (f==g) {
3972: (*f)(A,B,tol,flg);
3973: } else {
3974: SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
3975: }
3976: }
3977: return(0);
3978: }
3982: /*@
3983: MatIsHermitianTranspose - Test whether a matrix is another one's Hermitian transpose,
3985: Collective on Mat
3987: Input Parameter:
3988: + A - the matrix to test
3989: - B - the matrix to test against, this can equal the first parameter
3991: Output Parameters:
3992: . flg - the result
3994: Notes:
3995: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
3996: has a running time of the order of the number of nonzeros; the parallel
3997: test involves parallel copies of the block-offdiagonal parts of the matrix.
3999: Level: intermediate
4001: Concepts: matrices^transposing, matrix^symmetry
4003: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian(), MatIsTranspose()
4004: @*/
4005: PetscErrorCode MatIsHermitianTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
4006: {
4007: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);
4013: PetscObjectQueryFunction((PetscObject)A,"MatIsHermitianTranspose_C",(void (**)(void))&f);
4014: PetscObjectQueryFunction((PetscObject)B,"MatIsHermitianTranspose_C",(void (**)(void))&g);
4015: if (f && g) {
4016: if (f==g) {
4017: (*f)(A,B,tol,flg);
4018: } else {
4019: SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for Hermitian test");
4020: }
4021: }
4022: return(0);
4023: }
4027: /*@
4028: MatPermute - Creates a new matrix with rows and columns permuted from the
4029: original.
4031: Collective on Mat
4033: Input Parameters:
4034: + mat - the matrix to permute
4035: . row - row permutation, each processor supplies only the permutation for its rows
4036: - col - column permutation, each processor needs the entire column permutation, that is
4037: this is the same size as the total number of columns in the matrix
4039: Output Parameters:
4040: . B - the permuted matrix
4042: Level: advanced
4044: Concepts: matrices^permuting
4046: .seealso: MatGetOrdering()
4047: @*/
4048: PetscErrorCode MatPermute(Mat mat,IS row,IS col,Mat *B)
4049: {
4058: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4059: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4060: if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"MatPermute not available for Mat type %s",((PetscObject)mat)->type_name);
4061: MatPreallocated(mat);
4063: (*mat->ops->permute)(mat,row,col,B);
4064: PetscObjectStateIncrease((PetscObject)*B);
4065: return(0);
4066: }
4070: /*@
4071: MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the
4072: original and sparsified to the prescribed tolerance.
4074: Collective on Mat
4076: Input Parameters:
4077: + A - The matrix to permute
4078: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
4079: . frac - The half-bandwidth as a fraction of the total size, or 0.0
4080: . tol - The drop tolerance
4081: . rowp - The row permutation
4082: - colp - The column permutation
4084: Output Parameter:
4085: . B - The permuted, sparsified matrix
4087: Level: advanced
4089: Note:
4090: The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
4091: restrict the half-bandwidth of the resulting matrix to 5% of the
4092: total matrix size.
4094: .keywords: matrix, permute, sparsify
4096: .seealso: MatGetOrdering(), MatPermute()
4097: @*/
4098: PetscErrorCode MatPermuteSparsify(Mat A, PetscInt band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
4099: {
4100: IS irowp, icolp;
4101: const PetscInt *rows, *cols;
4102: PetscInt M, N, locRowStart, locRowEnd;
4103: PetscInt nz, newNz;
4104: const PetscInt *cwork;
4105: PetscInt *cnew;
4106: const PetscScalar *vwork;
4107: PetscScalar *vnew;
4108: PetscInt bw, issize;
4109: PetscInt row, locRow, newRow, col, newCol;
4110: PetscErrorCode ierr;
4117: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
4118: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
4119: if (!A->ops->permutesparsify) {
4120: MatGetSize(A, &M, &N);
4121: MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
4122: ISGetSize(rowp, &issize);
4123: if (issize != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for row permutation, should be %D", issize, M);
4124: ISGetSize(colp, &issize);
4125: if (issize != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for column permutation, should be %D", issize, N);
4126: ISInvertPermutation(rowp, 0, &irowp);
4127: ISGetIndices(irowp, &rows);
4128: ISInvertPermutation(colp, 0, &icolp);
4129: ISGetIndices(icolp, &cols);
4130: PetscMalloc(N * sizeof(PetscInt), &cnew);
4131: PetscMalloc(N * sizeof(PetscScalar), &vnew);
4133: /* Setup bandwidth to include */
4134: if (band == PETSC_DECIDE) {
4135: if (frac <= 0.0)
4136: bw = (PetscInt) (M * 0.05);
4137: else
4138: bw = (PetscInt) (M * frac);
4139: } else {
4140: if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
4141: bw = band;
4142: }
4144: /* Put values into new matrix */
4145: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
4146: for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
4147: MatGetRow(A, row, &nz, &cwork, &vwork);
4148: newRow = rows[locRow]+locRowStart;
4149: for(col = 0, newNz = 0; col < nz; col++) {
4150: newCol = cols[cwork[col]];
4151: if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
4152: cnew[newNz] = newCol;
4153: vnew[newNz] = vwork[col];
4154: newNz++;
4155: }
4156: }
4157: MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
4158: MatRestoreRow(A, row, &nz, &cwork, &vwork);
4159: }
4160: PetscFree(cnew);
4161: PetscFree(vnew);
4162: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
4163: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
4164: ISRestoreIndices(irowp, &rows);
4165: ISRestoreIndices(icolp, &cols);
4166: ISDestroy(irowp);
4167: ISDestroy(icolp);
4168: } else {
4169: (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
4170: }
4171: PetscObjectStateIncrease((PetscObject)*B);
4172: return(0);
4173: }
4177: /*@
4178: MatEqual - Compares two matrices.
4180: Collective on Mat
4182: Input Parameters:
4183: + A - the first matrix
4184: - B - the second matrix
4186: Output Parameter:
4187: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
4189: Level: intermediate
4191: Concepts: matrices^equality between
4192: @*/
4193: PetscErrorCode MatEqual(Mat A,Mat B,PetscTruth *flg)
4194: {
4204: MatPreallocated(B);
4205: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4206: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4207: if (A->rmap->N != B->rmap->N || A->cmap->N != B->cmap->N) SETERRQ4(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);
4208: if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)A)->type_name);
4209: if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)B)->type_name);
4210: if (A->ops->equal != B->ops->equal) SETERRQ2(PETSC_ERR_ARG_INCOMP,"A is type: %s\nB is type: %s",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
4211: MatPreallocated(A);
4213: (*A->ops->equal)(A,B,flg);
4214: return(0);
4215: }
4219: /*@
4220: MatDiagonalScale - Scales a matrix on the left and right by diagonal
4221: matrices that are stored as vectors. Either of the two scaling
4222: matrices can be PETSC_NULL.
4224: Collective on Mat
4226: Input Parameters:
4227: + mat - the matrix to be scaled
4228: . l - the left scaling vector (or PETSC_NULL)
4229: - r - the right scaling vector (or PETSC_NULL)
4231: Notes:
4232: MatDiagonalScale() computes A = LAR, where
4233: L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)
4235: Level: intermediate
4237: Concepts: matrices^diagonal scaling
4238: Concepts: diagonal scaling of matrices
4240: .seealso: MatScale()
4241: @*/
4242: PetscErrorCode MatDiagonalScale(Mat mat,Vec l,Vec r)
4243: {
4249: if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4252: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4253: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4254: MatPreallocated(mat);
4256: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4257: (*mat->ops->diagonalscale)(mat,l,r);
4258: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4259: PetscObjectStateIncrease((PetscObject)mat);
4260: return(0);
4261: }
4265: /*@
4266: MatScale - Scales all elements of a matrix by a given number.
4268: Collective on Mat
4270: Input Parameters:
4271: + mat - the matrix to be scaled
4272: - a - the scaling value
4274: Output Parameter:
4275: . mat - the scaled matrix
4277: Level: intermediate
4279: Concepts: matrices^scaling all entries
4281: .seealso: MatDiagonalScale()
4282: @*/
4283: PetscErrorCode MatScale(Mat mat,PetscScalar a)
4284: {
4290: if (a != 1.0 && !mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4291: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4292: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4293: MatPreallocated(mat);
4295: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4296: if (a != 1.0) {
4297: (*mat->ops->scale)(mat,a);
4298: PetscObjectStateIncrease((PetscObject)mat);
4299: }
4300: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4301: return(0);
4302: }
4306: /*@
4307: MatNorm - Calculates various norms of a matrix.
4309: Collective on Mat
4311: Input Parameters:
4312: + mat - the matrix
4313: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
4315: Output Parameters:
4316: . nrm - the resulting norm
4318: Level: intermediate
4320: Concepts: matrices^norm
4321: Concepts: norm^of matrix
4322: @*/
4323: PetscErrorCode MatNorm(Mat mat,NormType type,PetscReal *nrm)
4324: {
4332: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4333: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4334: if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4335: MatPreallocated(mat);
4337: (*mat->ops->norm)(mat,type,nrm);
4338: return(0);
4339: }
4341: /*
4342: This variable is used to prevent counting of MatAssemblyBegin() that
4343: are called from within a MatAssemblyEnd().
4344: */
4345: static PetscInt MatAssemblyEnd_InUse = 0;
4348: /*@
4349: MatAssemblyBegin - Begins assembling the matrix. This routine should
4350: be called after completing all calls to MatSetValues().
4352: Collective on Mat
4354: Input Parameters:
4355: + mat - the matrix
4356: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4357:
4358: Notes:
4359: MatSetValues() generally caches the values. The matrix is ready to
4360: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4361: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4362: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4363: using the matrix.
4365: Level: beginner
4367: Concepts: matrices^assembling
4369: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
4370: @*/
4371: PetscErrorCode MatAssemblyBegin(Mat mat,MatAssemblyType type)
4372: {
4378: MatPreallocated(mat);
4379: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
4380: if (mat->assembled) {
4381: mat->was_assembled = PETSC_TRUE;
4382: mat->assembled = PETSC_FALSE;
4383: }
4384: if (!MatAssemblyEnd_InUse) {
4385: PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
4386: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4387: PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
4388: } else {
4389: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4390: }
4391: return(0);
4392: }
4396: /*@
4397: MatAssembled - Indicates if a matrix has been assembled and is ready for
4398: use; for example, in matrix-vector product.
4400: Collective on Mat
4402: Input Parameter:
4403: . mat - the matrix
4405: Output Parameter:
4406: . assembled - PETSC_TRUE or PETSC_FALSE
4408: Level: advanced
4410: Concepts: matrices^assembled?
4412: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
4413: @*/
4414: PetscErrorCode MatAssembled(Mat mat,PetscTruth *assembled)
4415: {
4420: *assembled = mat->assembled;
4421: return(0);
4422: }
4426: /*
4427: Processes command line options to determine if/how a matrix
4428: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
4429: */
4430: PetscErrorCode MatView_Private(Mat mat)
4431: {
4432: PetscErrorCode ierr;
4433: PetscTruth flg1,flg2,flg3,flg4,flg6,flg7,flg8;
4434: static PetscTruth incall = PETSC_FALSE;
4435: #if defined(PETSC_USE_SOCKET_VIEWER)
4436: PetscTruth flg5;
4437: #endif
4440: if (incall) return(0);
4441: incall = PETSC_TRUE;
4442: PetscOptionsBegin(((PetscObject)mat)->comm,((PetscObject)mat)->prefix,"Matrix Options","Mat");
4443: PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg1);
4444: PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg2);
4445: PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg3);
4446: PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg4);
4447: #if defined(PETSC_USE_SOCKET_VIEWER)
4448: PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg5);
4449: #endif
4450: PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg6);
4451: PetscOptionsName("-mat_view_draw","Draw the matrix nonzero structure","MatView",&flg7);
4452: PetscOptionsEnd();
4454: if (flg1) {
4455: PetscViewer viewer;
4457: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4458: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
4459: MatView(mat,viewer);
4460: PetscViewerPopFormat(viewer);
4461: }
4462: if (flg2) {
4463: PetscViewer viewer;
4465: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4466: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
4467: MatView(mat,viewer);
4468: PetscViewerPopFormat(viewer);
4469: }
4470: if (flg3) {
4471: PetscViewer viewer;
4473: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4474: MatView(mat,viewer);
4475: }
4476: if (flg4) {
4477: PetscViewer viewer;
4479: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4480: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
4481: MatView(mat,viewer);
4482: PetscViewerPopFormat(viewer);
4483: }
4484: #if defined(PETSC_USE_SOCKET_VIEWER)
4485: if (flg5) {
4486: MatView(mat,PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4487: PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4488: }
4489: #endif
4490: if (flg6) {
4491: MatView(mat,PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4492: PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4493: }
4494: if (flg7) {
4495: PetscOptionsHasName(((PetscObject)mat)->prefix,"-mat_view_contour",&flg8);
4496: if (flg8) {
4497: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm),PETSC_VIEWER_DRAW_CONTOUR);
4498: }
4499: MatView(mat,PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4500: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4501: if (flg8) {
4502: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4503: }
4504: }
4505: incall = PETSC_FALSE;
4506: return(0);
4507: }
4511: /*@
4512: MatAssemblyEnd - Completes assembling the matrix. This routine should
4513: be called after MatAssemblyBegin().
4515: Collective on Mat
4517: Input Parameters:
4518: + mat - the matrix
4519: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4521: Options Database Keys:
4522: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
4523: . -mat_view_info_detailed - Prints more detailed info
4524: . -mat_view - Prints matrix in ASCII format
4525: . -mat_view_matlab - Prints matrix in Matlab format
4526: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
4527: . -display <name> - Sets display name (default is host)
4528: . -draw_pause <sec> - Sets number of seconds to pause after display
4529: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
4530: . -viewer_socket_machine <machine>
4531: . -viewer_socket_port <port>
4532: . -mat_view_binary - save matrix to file in binary format
4533: - -viewer_binary_filename <name>
4535: Notes:
4536: MatSetValues() generally caches the values. The matrix is ready to
4537: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4538: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4539: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4540: using the matrix.
4542: Level: beginner
4544: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
4545: @*/
4546: PetscErrorCode MatAssemblyEnd(Mat mat,MatAssemblyType type)
4547: {
4548: PetscErrorCode ierr;
4549: static PetscInt inassm = 0;
4550: PetscTruth flg;
4556: inassm++;
4557: MatAssemblyEnd_InUse++;
4558: if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
4559: PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
4560: if (mat->ops->assemblyend) {
4561: (*mat->ops->assemblyend)(mat,type);
4562: }
4563: PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
4564: } else {
4565: if (mat->ops->assemblyend) {
4566: (*mat->ops->assemblyend)(mat,type);
4567: }
4568: }
4570: /* Flush assembly is not a true assembly */
4571: if (type != MAT_FLUSH_ASSEMBLY) {
4572: mat->assembled = PETSC_TRUE; mat->num_ass++;
4573: }
4574: mat->insertmode = NOT_SET_VALUES;
4575: MatAssemblyEnd_InUse--;
4576: PetscObjectStateIncrease((PetscObject)mat);
4577: if (!mat->symmetric_eternal) {
4578: mat->symmetric_set = PETSC_FALSE;
4579: mat->hermitian_set = PETSC_FALSE;
4580: mat->structurally_symmetric_set = PETSC_FALSE;
4581: }
4582: if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
4583: MatView_Private(mat);
4584: PetscOptionsHasName(((PetscObject)mat)->prefix,"-mat_is_symmetric",&flg);
4585: if (flg) {
4586: PetscReal tol = 0.0;
4587: PetscOptionsGetReal(((PetscObject)mat)->prefix,"-mat_is_symmetric",&tol,PETSC_NULL);
4588: MatIsSymmetric(mat,tol,&flg);
4589: if (flg) {
4590: PetscPrintf(((PetscObject)mat)->comm,"Matrix is symmetric (tolerance %G)\n",tol);
4591: } else {
4592: PetscPrintf(((PetscObject)mat)->comm,"Matrix is not symmetric (tolerance %G)\n",tol);
4593: }
4594: }
4595: }
4596: inassm--;
4597: return(0);
4598: }
4603: /*@
4604: MatCompress - Tries to store the matrix in as little space as
4605: possible. May fail if memory is already fully used, since it
4606: tries to allocate new space.
4608: Collective on Mat
4610: Input Parameters:
4611: . mat - the matrix
4613: Level: advanced
4615: @*/
4616: PetscErrorCode MatCompress(Mat mat)
4617: {
4623: MatPreallocated(mat);
4624: if (mat->ops->compress) {(*mat->ops->compress)(mat);}
4625: return(0);
4626: }
4630: /*@
4631: MatSetOption - Sets a parameter option for a matrix. Some options
4632: may be specific to certain storage formats. Some options
4633: determine how values will be inserted (or added). Sorted,
4634: row-oriented input will generally assemble the fastest. The default
4635: is row-oriented, nonsorted input.
4637: Collective on Mat
4639: Input Parameters:
4640: + mat - the matrix
4641: . option - the option, one of those listed below (and possibly others),
4642: - flg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE)
4644: Options Describing Matrix Structure:
4645: + MAT_SYMMETRIC - symmetric in terms of both structure and value
4646: . MAT_HERMITIAN - transpose is the complex conjugation
4647: . MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
4648: - MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag
4649: you set to be kept with all future use of the matrix
4650: including after MatAssemblyBegin/End() which could
4651: potentially change the symmetry structure, i.e. you
4652: KNOW the matrix will ALWAYS have the property you set.
4655: Options For Use with MatSetValues():
4656: Insert a logically dense subblock, which can be
4657: . MAT_ROW_ORIENTED - row-oriented (default)
4659: Note these options reflect the data you pass in with MatSetValues(); it has
4660: nothing to do with how the data is stored internally in the matrix
4661: data structure.
4663: When (re)assembling a matrix, we can restrict the input for
4664: efficiency/debugging purposes. These options include
4665: + MAT_NEW_NONZERO_LOCATIONS - additional insertions will be
4666: allowed if they generate a new nonzero
4667: . MAT_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
4668: . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
4669: . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
4670: - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
4672: Notes:
4673: Some options are relevant only for particular matrix types and
4674: are thus ignored by others. Other options are not supported by
4675: certain matrix types and will generate an error message if set.
4677: If using a Fortran 77 module to compute a matrix, one may need to
4678: use the column-oriented option (or convert to the row-oriented
4679: format).
4681: MAT_NEW_NONZERO_LOCATIONS set to PETSC_FALSE indicates that any add or insertion
4682: that would generate a new entry in the nonzero structure is instead
4683: ignored. Thus, if memory has not alredy been allocated for this particular
4684: data, then the insertion is ignored. For dense matrices, in which
4685: the entire array is allocated, no entries are ever ignored.
4686: Set after the first MatAssemblyEnd()
4688: MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
4689: that would generate a new entry in the nonzero structure instead produces
4690: an error. (Currently supported for AIJ and BAIJ formats only.)
4691: This is a useful flag when using SAME_NONZERO_PATTERN in calling
4692: KSPSetOperators() to ensure that the nonzero pattern truely does
4693: remain unchanged. Set after the first MatAssemblyEnd()
4695: MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
4696: that would generate a new entry that has not been preallocated will
4697: instead produce an error. (Currently supported for AIJ and BAIJ formats
4698: only.) This is a useful flag when debugging matrix memory preallocation.
4700: MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
4701: other processors should be dropped, rather than stashed.
4702: This is useful if you know that the "owning" processor is also
4703: always generating the correct matrix entries, so that PETSc need
4704: not transfer duplicate entries generated on another processor.
4705:
4706: MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
4707: searches during matrix assembly. When this flag is set, the hash table
4708: is created during the first Matrix Assembly. This hash table is
4709: used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
4710: to improve the searching of indices. MAT_NEW_NONZERO_LOCATIONS flag
4711: should be used with MAT_USE_HASH_TABLE flag. This option is currently
4712: supported by MATMPIBAIJ format only.
4714: MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
4715: are kept in the nonzero structure
4717: MAT_IGNORE_ZERO_ENTRIES - for AIJ/IS matrices this will stop zero values from creating
4718: a zero location in the matrix
4720: MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
4721: ROWBS matrix types
4723: Level: intermediate
4725: Concepts: matrices^setting options
4727: @*/
4728: PetscErrorCode MatSetOption(Mat mat,MatOption op,PetscTruth flg)
4729: {
4735: if (((int) op) < 0 || ((int) op) >= NUM_MAT_OPTIONS) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Options %d is out of range",(int)op);
4736: MatPreallocated(mat);
4737: switch (op) {
4738: case MAT_SYMMETRIC:
4739: mat->symmetric = flg;
4740: if (flg) mat->structurally_symmetric = PETSC_TRUE;
4741: mat->symmetric_set = PETSC_TRUE;
4742: mat->structurally_symmetric_set = flg;
4743: break;
4744: case MAT_HERMITIAN:
4745: mat->hermitian = flg;
4746: if (flg) mat->structurally_symmetric = PETSC_TRUE;
4747: mat->hermitian_set = PETSC_TRUE;
4748: mat->structurally_symmetric_set = flg;
4749: break;
4750: case MAT_STRUCTURALLY_SYMMETRIC:
4751: mat->structurally_symmetric = flg;
4752: mat->structurally_symmetric_set = PETSC_TRUE;
4753: break;
4754: case MAT_SYMMETRY_ETERNAL:
4755: mat->symmetric_eternal = flg;
4756: break;
4757: default:
4758: break;
4759: }
4760: if (mat->ops->setoption) {
4761: (*mat->ops->setoption)(mat,op,flg);
4762: }
4763: return(0);
4764: }
4768: /*@
4769: MatZeroEntries - Zeros all entries of a matrix. For sparse matrices
4770: this routine retains the old nonzero structure.
4772: Collective on Mat
4774: Input Parameters:
4775: . mat - the matrix
4777: Level: intermediate
4779: Concepts: matrices^zeroing
4781: .seealso: MatZeroRows()
4782: @*/
4783: PetscErrorCode MatZeroEntries(Mat mat)
4784: {
4790: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4791: if (mat->insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for matrices where you have set values but not yet assembled");
4792: if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4793: MatPreallocated(mat);
4795: PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
4796: (*mat->ops->zeroentries)(mat);
4797: PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
4798: PetscObjectStateIncrease((PetscObject)mat);
4799: return(0);
4800: }
4804: /*@C
4805: MatZeroRows - Zeros all entries (except possibly the main diagonal)
4806: of a set of rows of a matrix.
4808: Collective on Mat
4810: Input Parameters:
4811: + mat - the matrix
4812: . numRows - the number of rows to remove
4813: . rows - the global row indices
4814: - diag - value put in all diagonals of eliminated rows (0.0 will even eliminate diagonal entry)
4816: Notes:
4817: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4818: but does not release memory. For the dense and block diagonal
4819: formats this does not alter the nonzero structure.
4821: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4822: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4823: merely zeroed.
4825: The user can set a value in the diagonal entry (or for the AIJ and
4826: row formats can optionally remove the main diagonal entry from the
4827: nonzero structure as well, by passing 0.0 as the final argument).
4829: For the parallel case, all processes that share the matrix (i.e.,
4830: those in the communicator used for matrix creation) MUST call this
4831: routine, regardless of whether any rows being zeroed are owned by
4832: them.
4834: Each processor can indicate any rows in the entire matrix to be zeroed (i.e. each process does NOT have to
4835: list only rows local to itself).
4837: Level: intermediate
4839: Concepts: matrices^zeroing rows
4841: .seealso: MatZeroRowsIS(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4842: @*/
4843: PetscErrorCode MatZeroRows(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4844: {
4851: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4852: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4853: if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4854: MatPreallocated(mat);
4856: (*mat->ops->zerorows)(mat,numRows,rows,diag);
4857: MatView_Private(mat);
4858: PetscObjectStateIncrease((PetscObject)mat);
4859: return(0);
4860: }
4864: /*@C
4865: MatZeroRowsIS - Zeros all entries (except possibly the main diagonal)
4866: of a set of rows of a matrix.
4868: Collective on Mat
4870: Input Parameters:
4871: + mat - the matrix
4872: . is - index set of rows to remove
4873: - diag - value put in all diagonals of eliminated rows
4875: Notes:
4876: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4877: but does not release memory. For the dense and block diagonal
4878: formats this does not alter the nonzero structure.
4880: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4881: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4882: merely zeroed.
4884: The user can set a value in the diagonal entry (or for the AIJ and
4885: row formats can optionally remove the main diagonal entry from the
4886: nonzero structure as well, by passing 0.0 as the final argument).
4888: For the parallel case, all processes that share the matrix (i.e.,
4889: those in the communicator used for matrix creation) MUST call this
4890: routine, regardless of whether any rows being zeroed are owned by
4891: them.
4893: Each processor should list the rows that IT wants zeroed
4895: Level: intermediate
4897: Concepts: matrices^zeroing rows
4899: .seealso: MatZeroRows(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4900: @*/
4901: PetscErrorCode MatZeroRowsIS(Mat mat,IS is,PetscScalar diag)
4902: {
4903: PetscInt numRows;
4904: const PetscInt *rows;
4911: ISGetLocalSize(is,&numRows);
4912: ISGetIndices(is,&rows);
4913: MatZeroRows(mat,numRows,rows,diag);
4914: ISRestoreIndices(is,&rows);
4915: return(0);
4916: }
4920: /*@C
4921: MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
4922: of a set of rows of a matrix; using local numbering of rows.
4924: Collective on Mat
4926: Input Parameters:
4927: + mat - the matrix
4928: . numRows - the number of rows to remove
4929: . rows - the global row indices
4930: - diag - value put in all diagonals of eliminated rows
4932: Notes:
4933: Before calling MatZeroRowsLocal(), the user must first set the
4934: local-to-global mapping by calling MatSetLocalToGlobalMapping().
4936: For the AIJ matrix formats this removes the old nonzero structure,
4937: but does not release memory. For the dense and block diagonal
4938: formats this does not alter the nonzero structure.
4940: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4941: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4942: merely zeroed.
4944: The user can set a value in the diagonal entry (or for the AIJ and
4945: row formats can optionally remove the main diagonal entry from the
4946: nonzero structure as well, by passing 0.0 as the final argument).
4948: Level: intermediate
4950: Concepts: matrices^zeroing
4952: .seealso: MatZeroRows(), MatZeroRowsLocalIS(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
4953: @*/
4954: PetscErrorCode MatZeroRowsLocal(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4955: {
4962: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4963: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4964: MatPreallocated(mat);
4966: if (mat->ops->zerorowslocal) {
4967: (*mat->ops->zerorowslocal)(mat,numRows,rows,diag);
4968: } else {
4969: IS is, newis;
4970: const PetscInt *newRows;
4972: if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
4973: ISCreateGeneral(PETSC_COMM_SELF,numRows,rows,&is);
4974: ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
4975: ISGetIndices(newis,&newRows);
4976: (*mat->ops->zerorows)(mat,numRows,newRows,diag);
4977: ISRestoreIndices(newis,&newRows);
4978: ISDestroy(newis);
4979: ISDestroy(is);
4980: }
4981: PetscObjectStateIncrease((PetscObject)mat);
4982: return(0);
4983: }
4987: /*@C
4988: MatZeroRowsLocalIS - Zeros all entries (except possibly the main diagonal)
4989: of a set of rows of a matrix; using local numbering of rows.
4991: Collective on Mat
4993: Input Parameters:
4994: + mat - the matrix
4995: . is - index set of rows to remove
4996: - diag - value put in all diagonals of eliminated rows
4998: Notes:
4999: Before calling MatZeroRowsLocalIS(), the user must first set the
5000: local-to-global mapping by calling MatSetLocalToGlobalMapping().
5002: For the AIJ matrix formats this removes the old nonzero structure,
5003: but does not release memory. For the dense and block diagonal
5004: formats this does not alter the nonzero structure.
5006: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
5007: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
5008: merely zeroed.
5010: The user can set a value in the diagonal entry (or for the AIJ and
5011: row formats can optionally remove the main diagonal entry from the
5012: nonzero structure as well, by passing 0.0 as the final argument).
5014: Level: intermediate
5016: Concepts: matrices^zeroing
5018: .seealso: MatZeroRows(), MatZeroRowsLocal(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
5019: @*/
5020: PetscErrorCode MatZeroRowsLocalIS(Mat mat,IS is,PetscScalar diag)
5021: {
5023: PetscInt numRows;
5024: const PetscInt *rows;
5030: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
5031: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
5032: MatPreallocated(mat);
5034: ISGetLocalSize(is,&numRows);
5035: ISGetIndices(is,&rows);
5036: MatZeroRowsLocal(mat,numRows,rows,diag);
5037: ISRestoreIndices(is,&rows);
5038: return(0);
5039: }
5043: /*@
5044: MatGetSize - Returns the numbers of rows and columns in a matrix.
5046: Not Collective
5048: Input Parameter:
5049: . mat - the matrix
5051: Output Parameters:
5052: + m - the number of global rows
5053: - n - the number of global columns
5055: Note: both output parameters can be PETSC_NULL on input.
5057: Level: beginner
5059: Concepts: matrices^size
5061: .seealso: MatGetLocalSize()
5062: @*/
5063: PetscErrorCode MatGetSize(Mat mat,PetscInt *m,PetscInt* n)
5064: {
5067: if (m) *m = mat->rmap->N;
5068: if (n) *n = mat->cmap->N;
5069: return(0);
5070: }
5074: /*@
5075: MatGetLocalSize - Returns the number of rows and columns in a matrix
5076: stored locally. This information may be implementation dependent, so
5077: use with care.
5079: Not Collective
5081: Input Parameters:
5082: . mat - the matrix
5084: Output Parameters:
5085: + m - the number of local rows
5086: - n - the number of local columns
5088: Note: both output parameters can be PETSC_NULL on input.
5090: Level: beginner
5092: Concepts: matrices^local size
5094: .seealso: MatGetSize()
5095: @*/
5096: PetscErrorCode MatGetLocalSize(Mat mat,PetscInt *m,PetscInt* n)
5097: {
5102: if (m) *m = mat->rmap->n;
5103: if (n) *n = mat->cmap->n;
5104: return(0);
5105: }
5109: /*@
5110: MatGetOwnershipRangeColumn - Returns the range of matrix columns owned by
5111: this processor.
5113: Not Collective, unless matrix has not been allocated, then collective on Mat
5115: Input Parameters:
5116: . mat - the matrix
5118: Output Parameters:
5119: + m - the global index of the first local column
5120: - n - one more than the global index of the last local column
5122: Notes: both output parameters can be PETSC_NULL on input.
5124: Level: developer
5126: Concepts: matrices^column ownership
5128: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), MatGetOwnershipRangesColumn()
5130: @*/
5131: PetscErrorCode MatGetOwnershipRangeColumn(Mat mat,PetscInt *m,PetscInt* n)
5132: {
5140: MatPreallocated(mat);
5141: if (m) *m = mat->cmap->rstart;
5142: if (n) *n = mat->cmap->rend;
5143: return(0);
5144: }
5148: /*@
5149: MatGetOwnershipRange - Returns the range of matrix rows owned by
5150: this processor, assuming that the matrix is laid out with the first
5151: n1 rows on the first processor, the next n2 rows on the second, etc.
5152: For certain parallel layouts this range may not be well defined.
5154: Not Collective, unless matrix has not been allocated, then collective on Mat
5156: Input Parameters:
5157: . mat - the matrix
5159: Output Parameters:
5160: + m - the global index of the first local row
5161: - n - one more than the global index of the last local row
5163: Note: both output parameters can be PETSC_NULL on input.
5165: Level: beginner
5167: Concepts: matrices^row ownership
5169: .seealso: MatGetOwnershipRanges(), MatGetOwnershipRangeColumn(), MatGetOwnershipRangesColumn()
5171: @*/
5172: PetscErrorCode MatGetOwnershipRange(Mat mat,PetscInt *m,PetscInt* n)
5173: {
5181: MatPreallocated(mat);
5182: if (m) *m = mat->rmap->rstart;
5183: if (n) *n = mat->rmap->rend;
5184: return(0);
5185: }
5189: /*@C
5190: MatGetOwnershipRanges - Returns the range of matrix rows owned by
5191: each process
5193: Not Collective, unless matrix has not been allocated, then collective on Mat
5195: Input Parameters:
5196: . mat - the matrix
5198: Output Parameters:
5199: . ranges - start of each processors portion plus one more then the total length at the end
5201: Level: beginner
5203: Concepts: matrices^row ownership
5205: .seealso: MatGetOwnershipRange(), MatGetOwnershipRangeColumn(), MatGetOwnershipRangesColumn()
5207: @*/
5208: PetscErrorCode MatGetOwnershipRanges(Mat mat,const PetscInt **ranges)
5209: {
5215: MatPreallocated(mat);
5216: PetscMapGetRanges(mat->rmap,ranges);
5217: return(0);
5218: }
5222: /*@C
5223: MatGetOwnershipRangesColumn - Returns the range of local columns for each process
5225: Not Collective, unless matrix has not been allocated, then collective on Mat
5227: Input Parameters:
5228: . mat - the matrix
5230: Output Parameters:
5231: . ranges - start of each processors portion plus one more then the total length at the end
5233: Level: beginner
5235: Concepts: matrices^column ownership
5237: .seealso: MatGetOwnershipRange(), MatGetOwnershipRangeColumn(), MatGetOwnershipRanges()
5239: @*/
5240: PetscErrorCode MatGetOwnershipRangesColumn(Mat mat,const PetscInt **ranges)
5241: {
5247: MatPreallocated(mat);
5248: PetscMapGetRanges(mat->cmap,ranges);
5249: return(0);
5250: }
5254: /*@C
5255: MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
5256: Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
5257: to complete the factorization.
5259: Collective