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.
1078: 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
1079: a single value per point) you can skip filling those indices.
1081: Inspired by the structured grid interface to the HYPRE package
1082: (http://www.llnl.gov/CASC/hypre)
1084: Efficiency Alert:
1085: The routine MatSetValuesBlockedStencil() may offer much better efficiency
1086: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
1088: Level: beginner
1090: Concepts: matrices^putting entries in
1092: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1093: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1094: @*/
1095: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1096: {
1098: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1099: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1102: if (!m || !n) return(0); /* no values to insert */
1109: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1110: if (n > 256) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1112: for (i=0; i<m; i++) {
1113: for (j=0; j<3-sdim; j++) dxm++;
1114: tmp = *dxm++ - starts[0];
1115: for (j=0; j<dim-1; j++) {
1116: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1117: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1118: }
1119: if (mat->stencil.noc) dxm++;
1120: jdxm[i] = tmp;
1121: }
1122: for (i=0; i<n; i++) {
1123: for (j=0; j<3-sdim; j++) dxn++;
1124: tmp = *dxn++ - starts[0];
1125: for (j=0; j<dim-1; j++) {
1126: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1127: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1128: }
1129: if (mat->stencil.noc) dxn++;
1130: jdxn[i] = tmp;
1131: }
1132: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1133: return(0);
1134: }
1138: /*@C
1139: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1140: Using structured grid indexing
1142: Not Collective
1144: Input Parameters:
1145: + mat - the matrix
1146: . v - a logically two-dimensional array of values
1147: . m - number of rows being entered
1148: . idxm - grid coordinates for matrix rows being entered
1149: . n - number of columns being entered
1150: . idxn - grid coordinates for matrix columns being entered
1151: - addv - either ADD_VALUES or INSERT_VALUES, where
1152: ADD_VALUES adds values to any existing entries, and
1153: INSERT_VALUES replaces existing entries with new values
1155: Notes:
1156: By default the values, v, are row-oriented and unsorted.
1157: See MatSetOption() for other options.
1159: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1160: options cannot be mixed without intervening calls to the assembly
1161: routines.
1163: The grid coordinates are across the entire grid, not just the local portion
1165: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1166: as well as in C.
1168: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
1170: In order to use this routine you must either obtain the matrix with DAGetMatrix()
1171: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1173: The columns and rows in the stencil passed in MUST be contained within the
1174: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
1175: if you create a DA with an overlap of one grid level and on a particular process its first
1176: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1177: first i index you can use in your column and row indices in MatSetStencil() is 5.
1179: In Fortran idxm and idxn should be declared as
1180: $ MatStencil idxm(4,m),idxn(4,n)
1181: and the values inserted using
1182: $ idxm(MatStencil_i,1) = i
1183: $ idxm(MatStencil_j,1) = j
1184: $ idxm(MatStencil_k,1) = k
1185: etc
1187: Negative indices may be passed in idxm and idxn, these rows and columns are
1188: simply ignored. This allows easily inserting element stiffness matrices
1189: with homogeneous Dirchlet boundary conditions that you don't want represented
1190: in the matrix.
1192: Inspired by the structured grid interface to the HYPRE package
1193: (http://www.llnl.gov/CASC/hypre)
1195: Level: beginner
1197: Concepts: matrices^putting entries in
1199: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1200: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1201: @*/
1202: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1203: {
1205: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1206: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1209: if (!m || !n) return(0); /* no values to insert */
1216: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1217: if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1219: for (i=0; i<m; i++) {
1220: for (j=0; j<3-sdim; j++) dxm++;
1221: tmp = *dxm++ - starts[0];
1222: for (j=0; j<sdim-1; j++) {
1223: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1224: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1225: }
1226: dxm++;
1227: jdxm[i] = tmp;
1228: }
1229: for (i=0; i<n; i++) {
1230: for (j=0; j<3-sdim; j++) dxn++;
1231: tmp = *dxn++ - starts[0];
1232: for (j=0; j<sdim-1; j++) {
1233: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1234: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1235: }
1236: dxn++;
1237: jdxn[i] = tmp;
1238: }
1239: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1240: return(0);
1241: }
1245: /*@
1246: MatSetStencil - Sets the grid information for setting values into a matrix via
1247: MatSetValuesStencil()
1249: Not Collective
1251: Input Parameters:
1252: + mat - the matrix
1253: . dim - dimension of the grid 1, 2, or 3
1254: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1255: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1256: - dof - number of degrees of freedom per node
1259: Inspired by the structured grid interface to the HYPRE package
1260: (www.llnl.gov/CASC/hyper)
1262: For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the
1263: user.
1264:
1265: Level: beginner
1267: Concepts: matrices^putting entries in
1269: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1270: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1271: @*/
1272: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1273: {
1274: PetscInt i;
1281: mat->stencil.dim = dim + (dof > 1);
1282: for (i=0; i<dim; i++) {
1283: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1284: mat->stencil.starts[i] = starts[dim-i-1];
1285: }
1286: mat->stencil.dims[dim] = dof;
1287: mat->stencil.starts[dim] = 0;
1288: mat->stencil.noc = (PetscTruth)(dof == 1);
1289: return(0);
1290: }
1294: /*@
1295: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1297: Not Collective
1299: Input Parameters:
1300: + mat - the matrix
1301: . v - a logically two-dimensional array of values
1302: . m, idxm - the number of block rows and their global block indices
1303: . n, idxn - the number of block columns and their global block indices
1304: - addv - either ADD_VALUES or INSERT_VALUES, where
1305: ADD_VALUES adds values to any existing entries, and
1306: INSERT_VALUES replaces existing entries with new values
1308: Notes:
1309: The m and n count the NUMBER of blocks in the row direction and column direction,
1310: NOT the total number of rows/columns; for example, if the block size is 2 and
1311: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1312: The values in idxm would be 1 2; that is the first index for each block divided by
1313: the block size.
1315: By default the values, v, are row-oriented and unsorted. So the layout of
1316: v is the same as for MatSetValues(). See MatSetOption() for other options.
1318: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1319: options cannot be mixed without intervening calls to the assembly
1320: routines.
1322: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1323: as well as in C.
1325: Negative indices may be passed in idxm and idxn, these rows and columns are
1326: simply ignored. This allows easily inserting element stiffness matrices
1327: with homogeneous Dirchlet boundary conditions that you don't want represented
1328: in the matrix.
1330: Each time an entry is set within a sparse matrix via MatSetValues(),
1331: internal searching must be done to determine where to place the the
1332: data in the matrix storage space. By instead inserting blocks of
1333: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1334: reduced.
1336: Example:
1337: $ Suppose m=n=2 and block size(bs) = 2 The array is
1338: $
1339: $ 1 2 | 3 4
1340: $ 5 6 | 7 8
1341: $ - - - | - - -
1342: $ 9 10 | 11 12
1343: $ 13 14 | 15 16
1344: $
1345: $ v[] should be passed in like
1346: $ v[] = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
1347: $
1348: $ If you are not using row oriented storage of v (that is you called MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE)) then
1349: $ v[] = [1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
1351: Level: intermediate
1353: Concepts: matrices^putting entries in blocked
1355: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1356: @*/
1357: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1358: {
1364: if (!m || !n) return(0); /* no values to insert */
1368: MatPreallocated(mat);
1369: if (mat->insertmode == NOT_SET_VALUES) {
1370: mat->insertmode = addv;
1371: }
1372: #if defined(PETSC_USE_DEBUG)
1373: else if (mat->insertmode != addv) {
1374: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1375: }
1376: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1377: #endif
1379: if (mat->assembled) {
1380: mat->was_assembled = PETSC_TRUE;
1381: mat->assembled = PETSC_FALSE;
1382: }
1383: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1384: if (mat->ops->setvaluesblocked) {
1385: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1386: } else {
1387: PetscInt buf[4096],*ibufm=0,*ibufn=0;
1388: PetscInt i,j,*iidxm,*iidxn,bs=mat->rmap->bs;
1389: if ((m+n)*bs <= 4096) {
1390: iidxm = buf; iidxn = buf + m*bs;
1391: } else {
1392: PetscMalloc2(m*bs,PetscInt,&ibufm,n*bs,PetscInt,&ibufn);
1393: iidxm = ibufm; iidxn = ibufn;
1394: }
1395: for (i=0; i<m; i++) {
1396: for (j=0; j<bs; j++) {
1397: iidxm[i*bs+j] = bs*idxm[i] + j;
1398: }
1399: }
1400: for (i=0; i<n; i++) {
1401: for (j=0; j<bs; j++) {
1402: iidxn[i*bs+j] = bs*idxn[i] + j;
1403: }
1404: }
1405: MatSetValues(mat,bs*m,iidxm,bs*n,iidxn,v,addv);
1406: PetscFree2(ibufm,ibufn);
1407: }
1408: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1409: return(0);
1410: }
1414: /*@
1415: MatGetValues - Gets a block of values from a matrix.
1417: Not Collective; currently only returns a local block
1419: Input Parameters:
1420: + mat - the matrix
1421: . v - a logically two-dimensional array for storing the values
1422: . m, idxm - the number of rows and their global indices
1423: - n, idxn - the number of columns and their global indices
1425: Notes:
1426: The user must allocate space (m*n PetscScalars) for the values, v.
1427: The values, v, are then returned in a row-oriented format,
1428: analogous to that used by default in MatSetValues().
1430: MatGetValues() uses 0-based row and column numbers in
1431: Fortran as well as in C.
1433: MatGetValues() requires that the matrix has been assembled
1434: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1435: MatSetValues() and MatGetValues() CANNOT be made in succession
1436: without intermediate matrix assembly.
1438: Negative row or column indices will be ignored and those locations in v[] will be
1439: left unchanged.
1441: Level: advanced
1443: Concepts: matrices^accessing values
1445: .seealso: MatGetRow(), MatGetSubMatrices(), MatSetValues()
1446: @*/
1447: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1448: {
1457: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1458: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1459: if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1460: MatPreallocated(mat);
1462: PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
1463: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1464: PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
1465: return(0);
1466: }
1470: /*@
1471: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1472: the routine MatSetValuesLocal() to allow users to insert matrix entries
1473: using a local (per-processor) numbering.
1475: Not Collective
1477: Input Parameters:
1478: + x - the matrix
1479: - mapping - mapping created with ISLocalToGlobalMappingCreate()
1480: or ISLocalToGlobalMappingCreateIS()
1482: Level: intermediate
1484: Concepts: matrices^local to global mapping
1485: Concepts: local to global mapping^for matrices
1487: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1488: @*/
1489: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
1490: {
1496: if (x->mapping) {
1497: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1498: }
1499: MatPreallocated(x);
1501: if (x->ops->setlocaltoglobalmapping) {
1502: (*x->ops->setlocaltoglobalmapping)(x,mapping);
1503: } else {
1504: PetscObjectReference((PetscObject)mapping);
1505: if (x->mapping) { ISLocalToGlobalMappingDestroy(x->mapping); }
1506: x->mapping = mapping;
1507: }
1508: return(0);
1509: }
1513: /*@
1514: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
1515: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
1516: entries using a local (per-processor) numbering.
1518: Not Collective
1520: Input Parameters:
1521: + x - the matrix
1522: - mapping - mapping created with ISLocalToGlobalMappingCreate() or
1523: ISLocalToGlobalMappingCreateIS()
1525: Level: intermediate
1527: Concepts: matrices^local to global mapping blocked
1528: Concepts: local to global mapping^for matrices, blocked
1530: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
1531: MatSetValuesBlocked(), MatSetValuesLocal()
1532: @*/
1533: PetscErrorCode MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
1534: {
1540: if (x->bmapping) {
1541: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1542: }
1543: PetscObjectReference((PetscObject)mapping);
1544: if (x->bmapping) { ISLocalToGlobalMappingDestroy(x->mapping); }
1545: x->bmapping = mapping;
1546: return(0);
1547: }
1551: /*@
1552: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
1553: using a local ordering of the nodes.
1555: Not Collective
1557: Input Parameters:
1558: + x - the matrix
1559: . nrow, irow - number of rows and their local indices
1560: . ncol, icol - number of columns and their local indices
1561: . y - a logically two-dimensional array of values
1562: - addv - either INSERT_VALUES or ADD_VALUES, where
1563: ADD_VALUES adds values to any existing entries, and
1564: INSERT_VALUES replaces existing entries with new values
1566: Notes:
1567: Before calling MatSetValuesLocal(), the user must first set the
1568: local-to-global mapping by calling MatSetLocalToGlobalMapping().
1570: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
1571: options cannot be mixed without intervening calls to the assembly
1572: routines.
1574: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1575: MUST be called after all calls to MatSetValuesLocal() have been completed.
1577: Level: intermediate
1579: Concepts: matrices^putting entries in with local numbering
1581: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
1582: MatSetValueLocal()
1583: @*/
1584: PetscErrorCode MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1585: {
1587: PetscInt irowm[2048],icolm[2048];
1592: if (!nrow || !ncol) return(0); /* no values to insert */
1596: MatPreallocated(mat);
1597: if (mat->insertmode == NOT_SET_VALUES) {
1598: mat->insertmode = addv;
1599: }
1600: #if defined(PETSC_USE_DEBUG)
1601: else if (mat->insertmode != addv) {
1602: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1603: }
1604: if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
1605: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1606: }
1607: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1608: #endif
1610: if (mat->assembled) {
1611: mat->was_assembled = PETSC_TRUE;
1612: mat->assembled = PETSC_FALSE;
1613: }
1614: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1615: if (!mat->ops->setvalueslocal) {
1616: ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
1617: ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
1618: (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
1619: } else {
1620: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1621: }
1622: mat->same_nonzero = PETSC_FALSE;
1623: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1624: return(0);
1625: }
1629: /*@
1630: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1631: using a local ordering of the nodes a block at a time.
1633: Not Collective
1635: Input Parameters:
1636: + x - the matrix
1637: . nrow, irow - number of rows and their local indices
1638: . ncol, icol - number of columns and their local indices
1639: . y - a logically two-dimensional array of values
1640: - addv - either INSERT_VALUES or ADD_VALUES, where
1641: ADD_VALUES adds values to any existing entries, and
1642: INSERT_VALUES replaces existing entries with new values
1644: Notes:
1645: Before calling MatSetValuesBlockedLocal(), the user must first set the
1646: local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
1647: where the mapping MUST be set for matrix blocks, not for matrix elements.
1649: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1650: options cannot be mixed without intervening calls to the assembly
1651: routines.
1653: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1654: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
1656: Level: intermediate
1658: Concepts: matrices^putting blocked values in with local numbering
1660: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
1661: @*/
1662: PetscErrorCode MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1663: {
1665: PetscInt irowm[2048],icolm[2048];
1670: if (!nrow || !ncol) return(0); /* no values to insert */
1674: MatPreallocated(mat);
1675: if (mat->insertmode == NOT_SET_VALUES) {
1676: mat->insertmode = addv;
1677: }
1678: #if defined(PETSC_USE_DEBUG)
1679: else if (mat->insertmode != addv) {
1680: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1681: }
1682: if (!mat->bmapping) {
1683: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
1684: }
1685: if (nrow > 2048 || ncol > 2048) {
1686: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1687: }
1688: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1689: #endif
1691: if (mat->assembled) {
1692: mat->was_assembled = PETSC_TRUE;
1693: mat->assembled = PETSC_FALSE;
1694: }
1695: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1696: ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
1697: ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
1698: if (mat->ops->setvaluesblocked) {
1699: (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
1700: } else {
1701: PetscInt buf[4096],*ibufm=0,*ibufn=0;
1702: PetscInt i,j,*iirowm,*iicolm,bs=mat->rmap->bs;
1703: if ((nrow+ncol)*bs <= 4096) {
1704: iirowm = buf; iicolm = buf + nrow*bs;
1705: } else {
1706: PetscMalloc2(nrow*bs,PetscInt,&ibufm,ncol*bs,PetscInt,&ibufn);
1707: iirowm = ibufm; iicolm = ibufn;
1708: }
1709: for (i=0; i<nrow; i++) {
1710: for (j=0; j<bs; j++) {
1711: iirowm[i*bs+j] = bs*irowm[i] + j;
1712: }
1713: }
1714: for (i=0; i<ncol; i++) {
1715: for (j=0; j<bs; j++) {
1716: iicolm[i*bs+j] = bs*icolm[i] + j;
1717: }
1718: }
1719: MatSetValues(mat,bs*nrow,iirowm,bs*ncol,iicolm,y,addv);
1720: PetscFree2(ibufm,ibufn);
1721: }
1722: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1723: return(0);
1724: }
1726: /* --------------------------------------------------------*/
1729: /*@
1730: MatMult - Computes the matrix-vector product, y = Ax.
1732: Collective on Mat and Vec
1734: Input Parameters:
1735: + mat - the matrix
1736: - x - the vector to be multiplied
1738: Output Parameters:
1739: . y - the result
1741: Notes:
1742: The vectors x and y cannot be the same. I.e., one cannot
1743: call MatMult(A,y,y).
1745: Level: beginner
1747: Concepts: matrix-vector product
1749: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
1750: @*/
1751: PetscErrorCode MatMult(Mat mat,Vec x,Vec y)
1752: {
1761: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1762: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1763: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1764: #ifndef PETSC_HAVE_CONSTRAINTS
1765: 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);
1766: 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);
1767: 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);
1768: #endif
1769: MatPreallocated(mat);
1771: if (mat->nullsp) {
1772: MatNullSpaceRemove(mat->nullsp,x,&x);
1773: }
1775: if (!mat->ops->mult) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
1776: PetscLogEventBegin(MAT_Mult,mat,x,y,0);
1777: (*mat->ops->mult)(mat,x,y);
1778: PetscLogEventEnd(MAT_Mult,mat,x,y,0);
1780: if (mat->nullsp) {
1781: MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
1782: }
1783: PetscObjectStateIncrease((PetscObject)y);
1784: return(0);
1785: }
1789: /*@
1790: MatMultTranspose - Computes matrix transpose times a vector.
1792: Collective on Mat and Vec
1794: Input Parameters:
1795: + mat - the matrix
1796: - x - the vector to be multilplied
1798: Output Parameters:
1799: . y - the result
1801: Notes:
1802: The vectors x and y cannot be the same. I.e., one cannot
1803: call MatMultTranspose(A,y,y).
1805: Level: beginner
1807: Concepts: matrix vector product^transpose
1809: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
1810: @*/
1811: PetscErrorCode MatMultTranspose(Mat mat,Vec x,Vec y)
1812: {
1821: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1822: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1823: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1824: #ifndef PETSC_HAVE_CONSTRAINTS
1825: 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);
1826: 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);
1827: #endif
1828: MatPreallocated(mat);
1830: if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
1831: PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
1832: (*mat->ops->multtranspose)(mat,x,y);
1833: PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
1834: PetscObjectStateIncrease((PetscObject)y);
1835: return(0);
1836: }
1840: /*@
1841: MatMultAdd - Computes v3 = v2 + A * v1.
1843: Collective on Mat and Vec
1845: Input Parameters:
1846: + mat - the matrix
1847: - v1, v2 - the vectors
1849: Output Parameters:
1850: . v3 - the result
1852: Notes:
1853: The vectors v1 and v3 cannot be the same. I.e., one cannot
1854: call MatMultAdd(A,v1,v2,v1).
1856: Level: beginner
1858: Concepts: matrix vector product^addition
1860: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1861: @*/
1862: PetscErrorCode MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1863: {
1873: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1874: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1875: 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);
1876: 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);
1877: 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);
1878: 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);
1879: 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);
1880: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1881: MatPreallocated(mat);
1883: PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1884: (*mat->ops->multadd)(mat,v1,v2,v3);
1885: PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1886: PetscObjectStateIncrease((PetscObject)v3);
1887: return(0);
1888: }
1892: /*@
1893: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1895: Collective on Mat and Vec
1897: Input Parameters:
1898: + mat - the matrix
1899: - v1, v2 - the vectors
1901: Output Parameters:
1902: . v3 - the result
1904: Notes:
1905: The vectors v1 and v3 cannot be the same. I.e., one cannot
1906: call MatMultTransposeAdd(A,v1,v2,v1).
1908: Level: beginner
1910: Concepts: matrix vector product^transpose and addition
1912: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1913: @*/
1914: PetscErrorCode MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1915: {
1925: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1926: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1927: if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
1928: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1929: 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);
1930: 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);
1931: 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);
1932: MatPreallocated(mat);
1934: PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1935: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
1936: PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1937: PetscObjectStateIncrease((PetscObject)v3);
1938: return(0);
1939: }
1943: /*@
1944: MatMultConstrained - The inner multiplication routine for a
1945: constrained matrix P^T A P.
1947: Collective on Mat and Vec
1949: Input Parameters:
1950: + mat - the matrix
1951: - x - the vector to be multilplied
1953: Output Parameters:
1954: . y - the result
1956: Notes:
1957: The vectors x and y cannot be the same. I.e., one cannot
1958: call MatMult(A,y,y).
1960: Level: beginner
1962: .keywords: matrix, multiply, matrix-vector product, constraint
1963: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1964: @*/
1965: PetscErrorCode MatMultConstrained(Mat mat,Vec x,Vec y)
1966: {
1973: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1974: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1975: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1976: 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);
1977: 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);
1978: 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);
1980: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1981: (*mat->ops->multconstrained)(mat,x,y);
1982: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
1983: PetscObjectStateIncrease((PetscObject)y);
1985: return(0);
1986: }
1990: /*@
1991: MatMultTransposeConstrained - The inner multiplication routine for a
1992: constrained matrix P^T A^T P.
1994: Collective on Mat and Vec
1996: Input Parameters:
1997: + mat - the matrix
1998: - x - the vector to be multilplied
2000: Output Parameters:
2001: . y - the result
2003: Notes:
2004: The vectors x and y cannot be the same. I.e., one cannot
2005: call MatMult(A,y,y).
2007: Level: beginner
2009: .keywords: matrix, multiply, matrix-vector product, constraint
2010: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
2011: @*/
2012: PetscErrorCode MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
2013: {
2020: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2021: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2022: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
2023: 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);
2024: 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);
2026: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
2027: (*mat->ops->multtransposeconstrained)(mat,x,y);
2028: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
2029: PetscObjectStateIncrease((PetscObject)y);
2031: return(0);
2032: }
2033: /* ------------------------------------------------------------*/
2036: /*@
2037: MatGetInfo - Returns information about matrix storage (number of
2038: nonzeros, memory, etc.).
2040: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
2041: as the flag
2043: Input Parameters:
2044: . mat - the matrix
2046: Output Parameters:
2047: + flag - flag indicating the type of parameters to be returned
2048: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
2049: MAT_GLOBAL_SUM - sum over all processors)
2050: - info - matrix information context
2052: Notes:
2053: The MatInfo context contains a variety of matrix data, including
2054: number of nonzeros allocated and used, number of mallocs during
2055: matrix assembly, etc. Additional information for factored matrices
2056: is provided (such as the fill ratio, number of mallocs during
2057: factorization, etc.). Much of this info is printed to PETSC_STDOUT
2058: when using the runtime options
2059: $ -info -mat_view_info
2061: Example for C/C++ Users:
2062: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
2063: data within the MatInfo context. For example,
2064: .vb
2065: MatInfo info;
2066: Mat A;
2067: double mal, nz_a, nz_u;
2069: MatGetInfo(A,MAT_LOCAL,&info);
2070: mal = info.mallocs;
2071: nz_a = info.nz_allocated;
2072: .ve
2074: Example for Fortran Users:
2075: Fortran users should declare info as a double precision
2076: array of dimension MAT_INFO_SIZE, and then extract the parameters
2077: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
2078: a complete list of parameter names.
2079: .vb
2080: double precision info(MAT_INFO_SIZE)
2081: double precision mal, nz_a
2082: Mat A
2083: integer ierr
2085: call MatGetInfo(A,MAT_LOCAL,info,ierr)
2086: mal = info(MAT_INFO_MALLOCS)
2087: nz_a = info(MAT_INFO_NZ_ALLOCATED)
2088: .ve
2090: Level: intermediate
2092: Concepts: matrices^getting information on
2093:
2094: @*/
2095: PetscErrorCode MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
2096: {
2103: if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2104: MatPreallocated(mat);
2105: (*mat->ops->getinfo)(mat,flag,info);
2106: return(0);
2107: }
2109: /* ----------------------------------------------------------*/
2112: /*@C
2113: MatILUDTFactor - Performs a drop tolerance ILU factorization.
2115: Collective on Mat
2117: Input Parameters:
2118: + mat - the matrix
2119: . row - row permutation
2120: . col - column permutation
2121: - info - information about the factorization to be done
2123: Output Parameters:
2124: . fact - the factored matrix
2126: Level: developer
2128: Notes:
2129: Most users should employ the simplified KSP interface for linear solvers
2130: instead of working directly with matrix algebra routines such as this.
2131: See, e.g., KSPCreate().
2133: This is currently only supported for the SeqAIJ matrix format using code
2134: from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or
2135: Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
2136: and thus can be distributed with PETSc.
2138: Concepts: matrices^ILUDT factorization
2140: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2141: @*/
2142: PetscErrorCode MatILUDTFactor(Mat mat,IS row,IS col,const MatFactorInfo *info,Mat *fact)
2143: {
2153: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2154: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2155: if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2156: MatPreallocated(mat);
2157: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2158: (*mat->ops->iludtfactor)(mat,row,col,info,fact);
2159: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2160: PetscObjectStateIncrease((PetscObject)*fact);
2162: return(0);
2163: }
2167: /*@
2168: MatLUFactor - Performs in-place LU factorization of matrix.
2170: Collective on Mat
2172: Input Parameters:
2173: + mat - the matrix
2174: . row - row permutation
2175: . col - column permutation
2176: - info - options for factorization, includes
2177: $ fill - expected fill as ratio of original fill.
2178: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2179: $ Run with the option -info to determine an optimal value to use
2181: Notes:
2182: Most users should employ the simplified KSP interface for linear solvers
2183: instead of working directly with matrix algebra routines such as this.
2184: See, e.g., KSPCreate().
2186: This changes the state of the matrix to a factored matrix; it cannot be used
2187: for example with MatSetValues() unless one first calls MatSetUnfactored().
2189: Level: developer
2191: Concepts: matrices^LU factorization
2193: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
2194: MatGetOrdering(), MatSetUnfactored(), MatFactorInfo
2196: @*/
2197: PetscErrorCode MatLUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2198: {
2207: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2208: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2209: if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2210: MatPreallocated(mat);
2212: PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
2213: (*mat->ops->lufactor)(mat,row,col,info);
2214: PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
2215: PetscObjectStateIncrease((PetscObject)mat);
2216: return(0);
2217: }
2221: /*@
2222: MatILUFactor - Performs in-place ILU factorization of matrix.
2224: Collective on Mat
2226: Input Parameters:
2227: + mat - the matrix
2228: . row - row permutation
2229: . col - column permutation
2230: - info - structure containing
2231: $ levels - number of levels of fill.
2232: $ expected fill - as ratio of original fill.
2233: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2234: missing diagonal entries)
2236: Notes:
2237: Probably really in-place only when level of fill is zero, otherwise allocates
2238: new space to store factored matrix and deletes previous memory.
2240: Most users should employ the simplified KSP interface for linear solvers
2241: instead of working directly with matrix algebra routines such as this.
2242: See, e.g., KSPCreate().
2244: Level: developer
2246: Concepts: matrices^ILU factorization
2248: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2249: @*/
2250: PetscErrorCode MatILUFactor(Mat mat,IS row,IS col,const MatFactorInfo *info)
2251: {
2260: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
2261: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2262: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2263: if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2264: MatPreallocated(mat);
2266: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
2267: (*mat->ops->ilufactor)(mat,row,col,info);
2268: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
2269: PetscObjectStateIncrease((PetscObject)mat);
2270: return(0);
2271: }
2275: /*@
2276: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
2277: Call this routine before calling MatLUFactorNumeric().
2279: Collective on Mat
2281: Input Parameters:
2282: + fact - the factor matrix obtained with MatGetFactor()
2283: . mat - the matrix
2284: . row, col - row and column permutations
2285: - info - options for factorization, includes
2286: $ fill - expected fill as ratio of original fill.
2287: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2288: $ Run with the option -info to determine an optimal value to use
2291: Notes:
2292: See the users manual for additional information about
2293: choosing the fill factor for better efficiency.
2295: Most users should employ the simplified KSP interface for linear solvers
2296: instead of working directly with matrix algebra routines such as this.
2297: See, e.g., KSPCreate().
2299: Level: developer
2301: Concepts: matrices^LU symbolic factorization
2303: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2304: @*/
2305: PetscErrorCode MatLUFactorSymbolic(Mat fact,Mat mat,IS row,IS col,const MatFactorInfo *info)
2306: {
2316: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2317: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2318: if (!(fact)->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",((PetscObject)mat)->type_name);
2319: MatPreallocated(mat);
2321: PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
2322: (fact->ops->lufactorsymbolic)(fact,mat,row,col,info);
2323: PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
2324: PetscObjectStateIncrease((PetscObject)fact);
2325: return(0);
2326: }
2330: /*@
2331: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
2332: Call this routine after first calling MatLUFactorSymbolic().
2334: Collective on Mat
2336: Input Parameters:
2337: + fact - the factor matrix obtained with MatGetFactor()
2338: . mat - the matrix
2339: - info - options for factorization
2341: Notes:
2342: See MatLUFactor() for in-place factorization. See
2343: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
2345: Most users should employ the simplified KSP interface for linear solvers
2346: instead of working directly with matrix algebra routines such as this.
2347: See, e.g., KSPCreate().
2349: Level: developer
2351: Concepts: matrices^LU numeric factorization
2353: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
2354: @*/
2355: PetscErrorCode MatLUFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2356: {
2364: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2365: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2366: 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);
2367: }
2368: if (!(fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2369: MatPreallocated(mat);
2370: PetscLogEventBegin(MAT_LUFactorNumeric,mat,fact,0,0);
2371: (fact->ops->lufactornumeric)(fact,mat,info);
2372: PetscLogEventEnd(MAT_LUFactorNumeric,mat,fact,0,0);
2374: MatView_Private(fact);
2375: PetscObjectStateIncrease((PetscObject)fact);
2376: return(0);
2377: }
2381: /*@
2382: MatCholeskyFactor - Performs in-place Cholesky factorization of a
2383: symmetric matrix.
2385: Collective on Mat
2387: Input Parameters:
2388: + mat - the matrix
2389: . perm - row and column permutations
2390: - f - expected fill as ratio of original fill
2392: Notes:
2393: See MatLUFactor() for the nonsymmetric case. See also
2394: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
2396: Most users should employ the simplified KSP interface for linear solvers
2397: instead of working directly with matrix algebra routines such as this.
2398: See, e.g., KSPCreate().
2400: Level: developer
2402: Concepts: matrices^Cholesky factorization
2404: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
2405: MatGetOrdering()
2407: @*/
2408: PetscErrorCode MatCholeskyFactor(Mat mat,IS perm,const MatFactorInfo *info)
2409: {
2417: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2418: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2419: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2420: if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2421: MatPreallocated(mat);
2423: PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
2424: (*mat->ops->choleskyfactor)(mat,perm,info);
2425: PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
2426: PetscObjectStateIncrease((PetscObject)mat);
2427: return(0);
2428: }
2432: /*@
2433: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
2434: of a symmetric matrix.
2436: Collective on Mat
2438: Input Parameters:
2439: + fact - the factor matrix obtained with MatGetFactor()
2440: . mat - the matrix
2441: . perm - row and column permutations
2442: - info - options for factorization, includes
2443: $ fill - expected fill as ratio of original fill.
2444: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2445: $ Run with the option -info to determine an optimal value to use
2447: Notes:
2448: See MatLUFactorSymbolic() for the nonsymmetric case. See also
2449: MatCholeskyFactor() and MatCholeskyFactorNumeric().
2451: Most users should employ the simplified KSP interface for linear solvers
2452: instead of working directly with matrix algebra routines such as this.
2453: See, e.g., KSPCreate().
2455: Level: developer
2457: Concepts: matrices^Cholesky symbolic factorization
2459: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
2460: MatGetOrdering()
2462: @*/
2463: PetscErrorCode MatCholeskyFactorSymbolic(Mat fact,Mat mat,IS perm,const MatFactorInfo *info)
2464: {
2473: if (mat->rmap->N != mat->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2474: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2475: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2476: if (!(fact)->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2477: MatPreallocated(mat);
2479: PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2480: (fact->ops->choleskyfactorsymbolic)(fact,mat,perm,info);
2481: PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
2482: PetscObjectStateIncrease((PetscObject)fact);
2483: return(0);
2484: }
2488: /*@
2489: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
2490: of a symmetric matrix. Call this routine after first calling
2491: MatCholeskyFactorSymbolic().
2493: Collective on Mat
2495: Input Parameters:
2496: + fact - the factor matrix obtained with MatGetFactor()
2497: . mat - the initial matrix
2498: . info - options for factorization
2499: - fact - the symbolic factor of mat
2502: Notes:
2503: Most users should employ the simplified KSP interface for linear solvers
2504: instead of working directly with matrix algebra routines such as this.
2505: See, e.g., KSPCreate().
2507: Level: developer
2509: Concepts: matrices^Cholesky numeric factorization
2511: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
2512: @*/
2513: PetscErrorCode MatCholeskyFactorNumeric(Mat fact,Mat mat,const MatFactorInfo *info)
2514: {
2522: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2523: if (!(fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2524: if (mat->rmap->N != (fact)->rmap->N || mat->cmap->N != (fact)->cmap->N) {
2525: 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);
2526: }
2527: MatPreallocated(mat);
2529: PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,fact,0,0);
2530: (fact->ops->choleskyfactornumeric)(fact,mat,info);
2531: PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,fact,0,0);
2533: MatView_Private(fact);
2534: PetscObjectStateIncrease((PetscObject)fact);
2535: return(0);
2536: }
2538: /* ----------------------------------------------------------------*/
2541: /*@
2542: MatSolve - Solves A x = b, given a factored matrix.
2544: Collective on Mat and Vec
2546: Input Parameters:
2547: + mat - the factored matrix
2548: - b - the right-hand-side vector
2550: Output Parameter:
2551: . x - the result vector
2553: Notes:
2554: The vectors b and x cannot be the same. I.e., one cannot
2555: call MatSolve(A,x,x).
2557: Notes:
2558: Most users should employ the simplified KSP interface for linear solvers
2559: instead of working directly with matrix algebra routines such as this.
2560: See, e.g., KSPCreate().
2562: Level: developer
2564: Concepts: matrices^triangular solves
2566: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
2567: @*/
2568: PetscErrorCode MatSolve(Mat mat,Vec b,Vec x)
2569: {
2579: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2580: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2581: 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);
2582: 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);
2583: 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);
2584: if (!mat->rmap->N && !mat->cmap->N) return(0);
2585: if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2586: MatPreallocated(mat);
2588: PetscLogEventBegin(MAT_Solve,mat,b,x,0);
2589: (*mat->ops->solve)(mat,b,x);
2590: PetscLogEventEnd(MAT_Solve,mat,b,x,0);
2591: PetscObjectStateIncrease((PetscObject)x);
2592: return(0);
2593: }
2597: PetscErrorCode MatMatSolve_Basic(Mat A,Mat B,Mat X)
2598: {
2600: Vec b,x;
2601: PetscInt m,N,i;
2602: PetscScalar *bb,*xx;
2605: MatGetArray(B,&bb);
2606: MatGetArray(X,&xx);
2607: MatGetLocalSize(B,&m,PETSC_NULL); /* number local rows */
2608: MatGetSize(B,PETSC_NULL,&N); /* total columns in dense matrix */
2609: VecCreateMPIWithArray(((PetscObject)A)->comm,m,PETSC_DETERMINE,PETSC_NULL,&b);
2610: VecCreateMPIWithArray(((PetscObject)A)->comm,m,PETSC_DETERMINE,PETSC_NULL,&x);
2611: for (i=0; i<N; i++) {
2612: VecPlaceArray(b,bb + i*m);
2613: VecPlaceArray(x,xx + i*m);
2614: MatSolve(A,b,x);
2615: VecResetArray(x);
2616: VecResetArray(b);
2617: }
2618: VecDestroy(b);
2619: VecDestroy(x);
2620: MatRestoreArray(B,&bb);
2621: MatRestoreArray(X,&xx);
2622: return(0);
2623: }
2627: /*@
2628: MatMatSolve - Solves A X = B, given a factored matrix.
2630: Collective on Mat
2632: Input Parameters:
2633: + mat - the factored matrix
2634: - B - the right-hand-side matrix (dense matrix)
2636: Output Parameter:
2637: . X - the result matrix (dense matrix)
2639: Notes:
2640: The matrices b and x cannot be the same. I.e., one cannot
2641: call MatMatSolve(A,x,x).
2643: Notes:
2644: Most users should usually employ the simplified KSP interface for linear solvers
2645: instead of working directly with matrix algebra routines such as this.
2646: See, e.g., KSPCreate(). However KSP can only solve for one vector (column of X)
2647: at a time.
2649: Level: developer
2651: Concepts: matrices^triangular solves
2653: .seealso: MatMatSolveAdd(), MatMatSolveTranspose(), MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()
2654: @*/
2655: PetscErrorCode MatMatSolve(Mat A,Mat B,Mat X)
2656: {
2666: if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different matrices");
2667: if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2668: 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);
2669: 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);
2670: 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);
2671: if (!A->rmap->N && !A->cmap->N) return(0);
2672: MatPreallocated(A);
2674: PetscLogEventBegin(MAT_MatSolve,A,B,X,0);
2675: if (!A->ops->matsolve) {
2676: PetscInfo1(A,"Mat type %s using basic MatMatSolve",((PetscObject)A)->type_name);
2677: MatMatSolve_Basic(A,B,X);
2678: } else {
2679: (*A->ops->matsolve)(A,B,X);
2680: }
2681: PetscLogEventEnd(MAT_MatSolve,A,B,X,0);
2682: PetscObjectStateIncrease((PetscObject)X);
2683: return(0);
2684: }
2689: /* @
2690: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU, or
2691: U^T*D^(1/2) x = b, given a factored symmetric matrix, A = U^T*D*U,
2693: Collective on Mat and Vec
2695: Input Parameters:
2696: + mat - the factored matrix
2697: - b - the right-hand-side vector
2699: Output Parameter:
2700: . x - the result vector
2702: Notes:
2703: MatSolve() should be used for most applications, as it performs
2704: a forward solve followed by a backward solve.
2706: The vectors b and x cannot be the same, i.e., one cannot
2707: call MatForwardSolve(A,x,x).
2709: For matrix in seqsbaij format with block size larger than 1,
2710: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2711: MatForwardSolve() solves U^T*D y = b, and
2712: MatBackwardSolve() solves U x = y.
2713: Thus they do not provide a symmetric preconditioner.
2715: Most users should employ the simplified KSP interface for linear solvers
2716: instead of working directly with matrix algebra routines such as this.
2717: See, e.g., KSPCreate().
2719: Level: developer
2721: Concepts: matrices^forward solves
2723: .seealso: MatSolve(), MatBackwardSolve()
2724: @ */
2725: PetscErrorCode MatForwardSolve(Mat mat,Vec b,Vec x)
2726: {
2736: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2737: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2738: if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2739: 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);
2740: 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);
2741: 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);
2742: MatPreallocated(mat);
2743: PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
2744: (*mat->ops->forwardsolve)(mat,b,x);
2745: PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
2746: PetscObjectStateIncrease((PetscObject)x);
2747: return(0);
2748: }
2752: /* @
2753: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
2754: D^(1/2) U x = b, given a factored symmetric matrix, A = U^T*D*U,
2756: Collective on Mat and Vec
2758: Input Parameters:
2759: + mat - the factored matrix
2760: - b - the right-hand-side vector
2762: Output Parameter:
2763: . x - the result vector
2765: Notes:
2766: MatSolve() should be used for most applications, as it performs
2767: a forward solve followed by a backward solve.
2769: The vectors b and x cannot be the same. I.e., one cannot
2770: call MatBackwardSolve(A,x,x).
2772: For matrix in seqsbaij format with block size larger than 1,
2773: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2774: MatForwardSolve() solves U^T*D y = b, and
2775: MatBackwardSolve() solves U x = y.
2776: Thus they do not provide a symmetric preconditioner.
2778: Most users should employ the simplified KSP interface for linear solvers
2779: instead of working directly with matrix algebra routines such as this.
2780: See, e.g., KSPCreate().
2782: Level: developer
2784: Concepts: matrices^backward solves
2786: .seealso: MatSolve(), MatForwardSolve()
2787: @ */
2788: PetscErrorCode MatBackwardSolve(Mat mat,Vec b,Vec x)
2789: {
2799: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2800: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2801: if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
2802: 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);
2803: 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);
2804: 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);
2805: MatPreallocated(mat);
2807: PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
2808: (*mat->ops->backwardsolve)(mat,b,x);
2809: PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
2810: PetscObjectStateIncrease((PetscObject)x);
2811: return(0);
2812: }
2816: /*@
2817: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
2819: Collective on Mat and Vec
2821: Input Parameters:
2822: + mat - the factored matrix
2823: . b - the right-hand-side vector
2824: - y - the vector to be added to
2826: Output Parameter:
2827: . x - the result vector
2829: Notes:
2830: The vectors b and x cannot be the same. I.e., one cannot
2831: call MatSolveAdd(A,x,y,x).
2833: Most users should employ the simplified KSP interface for linear solvers
2834: instead of working directly with matrix algebra routines such as this.
2835: See, e.g., KSPCreate().
2837: Level: developer
2839: Concepts: matrices^triangular solves
2841: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
2842: @*/
2843: PetscErrorCode MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
2844: {
2845: PetscScalar one = 1.0;
2846: Vec tmp;
2858: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2859: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2860: 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);
2861: 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);
2862: 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);
2863: 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);
2864: 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);
2865: MatPreallocated(mat);
2867: PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
2868: if (mat->ops->solveadd) {
2869: (*mat->ops->solveadd)(mat,b,y,x);
2870: } else {
2871: /* do the solve then the add manually */
2872: if (x != y) {
2873: MatSolve(mat,b,x);
2874: VecAXPY(x,one,y);
2875: } else {
2876: VecDuplicate(x,&tmp);
2877: PetscLogObjectParent(mat,tmp);
2878: VecCopy(x,tmp);
2879: MatSolve(mat,b,x);
2880: VecAXPY(x,one,tmp);
2881: VecDestroy(tmp);
2882: }
2883: }
2884: PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
2885: PetscObjectStateIncrease((PetscObject)x);
2886: return(0);
2887: }
2891: /*@
2892: MatSolveTranspose - Solves A' x = b, given a factored matrix.
2894: Collective on Mat and Vec
2896: Input Parameters:
2897: + mat - the factored matrix
2898: - b - the right-hand-side vector
2900: Output Parameter:
2901: . x - the result vector
2903: Notes:
2904: The vectors b and x cannot be the same. I.e., one cannot
2905: call MatSolveTranspose(A,x,x).
2907: Most users should employ the simplified KSP interface for linear solvers
2908: instead of working directly with matrix algebra routines such as this.
2909: See, e.g., KSPCreate().
2911: Level: developer
2913: Concepts: matrices^triangular solves
2915: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
2916: @*/
2917: PetscErrorCode MatSolveTranspose(Mat mat,Vec b,Vec x)
2918: {
2928: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2929: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2930: if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",((PetscObject)mat)->type_name);
2931: 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);
2932: 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);
2933: MatPreallocated(mat);
2934: PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
2935: (*mat->ops->solvetranspose)(mat,b,x);
2936: PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
2937: PetscObjectStateIncrease((PetscObject)x);
2938: return(0);
2939: }
2943: /*@
2944: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
2945: factored matrix.
2947: Collective on Mat and Vec
2949: Input Parameters:
2950: + mat - the factored matrix
2951: . b - the right-hand-side vector
2952: - y - the vector to be added to
2954: Output Parameter:
2955: . x - the result vector
2957: Notes:
2958: The vectors b and x cannot be the same. I.e., one cannot
2959: call MatSolveTransposeAdd(A,x,y,x).
2961: Most users should employ the simplified KSP interface for linear solvers
2962: instead of working directly with matrix algebra routines such as this.
2963: See, e.g., KSPCreate().
2965: Level: developer
2967: Concepts: matrices^triangular solves
2969: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
2970: @*/
2971: PetscErrorCode MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
2972: {
2973: PetscScalar one = 1.0;
2975: Vec tmp;
2986: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2987: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2988: 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);
2989: 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);
2990: 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);
2991: 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);
2992: MatPreallocated(mat);
2994: PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
2995: if (mat->ops->solvetransposeadd) {
2996: (*mat->ops->solvetransposeadd)(mat,b,y,x);
2997: } else {
2998: /* do the solve then the add manually */
2999: if (x != y) {
3000: MatSolveTranspose(mat,b,x);
3001: VecAXPY(x,one,y);
3002: } else {
3003: VecDuplicate(x,&tmp);
3004: PetscLogObjectParent(mat,tmp);
3005: VecCopy(x,tmp);
3006: MatSolveTranspose(mat,b,x);
3007: VecAXPY(x,one,tmp);
3008: VecDestroy(tmp);
3009: }
3010: }
3011: PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
3012: PetscObjectStateIncrease((PetscObject)x);
3013: return(0);
3014: }
3015: /* ----------------------------------------------------------------*/
3019: /*@
3020: MatRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3022: Collective on Mat and Vec
3024: Input Parameters:
3025: + mat - the matrix
3026: . b - the right hand side
3027: . omega - the relaxation factor
3028: . flag - flag indicating the type of SOR (see below)
3029: . shift - diagonal shift
3030: . its - the number of iterations
3031: - lits - the number of local iterations
3033: Output Parameters:
3034: . x - the solution (can contain an initial guess)
3036: SOR Flags:
3037: . SOR_FORWARD_SWEEP - forward SOR
3038: . SOR_BACKWARD_SWEEP - backward SOR
3039: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
3040: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
3041: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
3042: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
3043: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
3044: upper/lower triangular part of matrix to
3045: vector (with omega)
3046: . SOR_ZERO_INITIAL_GUESS - zero initial guess
3048: Notes:
3049: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
3050: SOR_LOCAL_SYMMETRIC_SWEEP perform separate independent smoothings
3051: on each processor.
3053: Application programmers will not generally use MatRelax() directly,
3054: but instead will employ the KSP/PC interface.
3056: Notes for Advanced Users:
3057: The flags are implemented as bitwise inclusive or operations.
3058: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
3059: to specify a zero initial guess for SSOR.
3061: Most users should employ the simplified KSP interface for linear solvers
3062: instead of working directly with matrix algebra routines such as this.
3063: See, e.g., KSPCreate().
3065: See also, MatPBRelax(). This routine will automatically call the point block
3066: version if the point version is not available.
3068: Level: developer
3070: Concepts: matrices^relaxation
3071: Concepts: matrices^SOR
3072: Concepts: matrices^Gauss-Seidel
3074: @*/
3075: PetscErrorCode MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3076: {
3086: if (!mat->ops->relax && !mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3087: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3088: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3089: 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);
3090: 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);
3091: 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);
3092: if (its <= 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
3093: if (lits <= 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Relaxation requires local its %D positive",lits);
3095: MatPreallocated(mat);
3096: PetscLogEventBegin(MAT_Relax,mat,b,x,0);
3097: if (mat->ops->relax) {
3098: ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
3099: } else {
3100: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
3101: }
3102: PetscLogEventEnd(MAT_Relax,mat,b,x,0);
3103: PetscObjectStateIncrease((PetscObject)x);
3104: return(0);
3105: }
3109: /*@
3110: MatPBRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
3112: Collective on Mat and Vec
3114: See MatRelax() for usage
3116: For multi-component PDEs where the Jacobian is stored in a point block format
3117: (with the PETSc BAIJ matrix formats) the relaxation is done one point block at
3118: a time. That is, the small (for example, 4 by 4) blocks along the diagonal are solved
3119: simultaneously (that is a 4 by 4 linear solve is done) to update all the values at a point.
3121: Level: developer
3123: @*/
3124: PetscErrorCode MatPBRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
3125: {
3135: if (!mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3136: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3137: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3138: 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);
3139: 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);
3140: 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);
3141: MatPreallocated(mat);
3143: PetscLogEventBegin(MAT_Relax,mat,b,x,0);
3144: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
3145: PetscLogEventEnd(MAT_Relax,mat,b,x,0);
3146: PetscObjectStateIncrease((PetscObject)x);
3147: return(0);
3148: }
3152: /*
3153: Default matrix copy routine.
3154: */
3155: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
3156: {
3157: PetscErrorCode ierr;
3158: PetscInt i,rstart,rend,nz;
3159: const PetscInt *cwork;
3160: const PetscScalar *vwork;
3163: if (B->assembled) {
3164: MatZeroEntries(B);
3165: }
3166: MatGetOwnershipRange(A,&rstart,&rend);
3167: for (i=rstart; i<rend; i++) {
3168: MatGetRow(A,i,&nz,&cwork,&vwork);
3169: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
3170: MatRestoreRow(A,i,&nz,&cwork,&vwork);
3171: }
3172: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3173: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3174: PetscObjectStateIncrease((PetscObject)B);
3175: return(0);
3176: }
3180: /*@
3181: MatCopy - Copys a matrix to another matrix.
3183: Collective on Mat
3185: Input Parameters:
3186: + A - the matrix
3187: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
3189: Output Parameter:
3190: . B - where the copy is put
3192: Notes:
3193: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
3194: same nonzero pattern or the routine will crash.
3196: MatCopy() copies the matrix entries of a matrix to another existing
3197: matrix (after first zeroing the second matrix). A related routine is
3198: MatConvert(), which first creates a new matrix and then copies the data.
3200: Level: intermediate
3201:
3202: Concepts: matrices^copying
3204: .seealso: MatConvert(), MatDuplicate()
3206: @*/
3207: PetscErrorCode MatCopy(Mat A,Mat B,MatStructure str)
3208: {
3210: PetscInt i;
3218: MatPreallocated(B);
3219: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3220: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3221: 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);
3222: MatPreallocated(A);
3224: PetscLogEventBegin(MAT_Copy,A,B,0,0);
3225: if (A->ops->copy) {
3226: (*A->ops->copy)(A,B,str);
3227: } else { /* generic conversion */
3228: MatCopy_Basic(A,B,str);
3229: }
3230: if (A->mapping) {
3231: if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
3232: MatSetLocalToGlobalMapping(B,A->mapping);
3233: }
3234: if (A->bmapping) {
3235: if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
3236: MatSetLocalToGlobalMappingBlock(B,A->mapping);
3237: }
3239: B->stencil.dim = A->stencil.dim;
3240: B->stencil.noc = A->stencil.noc;
3241: for (i=0; i<=A->stencil.dim; i++) {
3242: B->stencil.dims[i] = A->stencil.dims[i];
3243: B->stencil.starts[i] = A->stencil.starts[i];
3244: }
3246: PetscLogEventEnd(MAT_Copy,A,B,0,0);
3247: PetscObjectStateIncrease((PetscObject)B);
3248: return(0);
3249: }
3253: /*@C
3254: MatConvert - Converts a matrix to another matrix, either of the same
3255: or different type.
3257: Collective on Mat
3259: Input Parameters:
3260: + mat - the matrix
3261: . newtype - new matrix type. Use MATSAME to create a new matrix of the
3262: same type as the original matrix.
3263: - reuse - denotes if the destination matrix is to be created or reused. Currently
3264: MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
3265: MAT_INITIAL_MATRIX.
3267: Output Parameter:
3268: . M - pointer to place new matrix
3270: Notes:
3271: MatConvert() first creates a new matrix and then copies the data from
3272: the first matrix. A related routine is MatCopy(), which copies the matrix
3273: entries of one matrix to another already existing matrix context.
3275: Cannot be used to convert a sequential matrix to parallel or parallel to sequential,
3276: the MPI communicator of the generated matrix is always the same as the communicator
3277: of the input matrix.
3279: Level: intermediate
3281: Concepts: matrices^converting between storage formats
3283: .seealso: MatCopy(), MatDuplicate()
3284: @*/
3285: PetscErrorCode MatConvert(Mat mat, const MatType newtype,MatReuse reuse,Mat *M)
3286: {
3287: PetscErrorCode ierr;
3288: PetscTruth sametype,issame,flg;
3289: char convname[256],mtype[256];
3290: Mat B;
3296: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3297: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3298: MatPreallocated(mat);
3300: PetscOptionsGetString(((PetscObject)mat)->prefix,"-matconvert_type",mtype,256,&flg);
3301: if (flg) {
3302: newtype = mtype;
3303: }
3304: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
3305: PetscStrcmp(newtype,"same",&issame);
3306: if ((reuse == MAT_REUSE_MATRIX) && (mat != *M)) {
3307: SETERRQ(PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for in-place conversion currently");
3308: }
3310: if ((reuse == MAT_REUSE_MATRIX) && (issame || sametype)) return(0);
3311:
3312: if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
3313: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
3314: } else {
3315: PetscErrorCode (*conv)(Mat, const MatType,MatReuse,Mat*)=PETSC_NULL;
3316: const char *prefix[3] = {"seq","mpi",""};
3317: PetscInt i;
3318: /*
3319: Order of precedence:
3320: 1) See if a specialized converter is known to the current matrix.
3321: 2) See if a specialized converter is known to the desired matrix class.
3322: 3) See if a good general converter is registered for the desired class
3323: (as of 6/27/03 only MATMPIADJ falls into this category).
3324: 4) See if a good general converter is known for the current matrix.
3325: 5) Use a really basic converter.
3326: */
3327:
3328: /* 1) See if a specialized converter is known to the current matrix and the desired class */
3329: for (i=0; i<3; i++) {
3330: PetscStrcpy(convname,"MatConvert_");
3331: PetscStrcat(convname,((PetscObject)mat)->type_name);
3332: PetscStrcat(convname,"_");
3333: PetscStrcat(convname,prefix[i]);
3334: PetscStrcat(convname,newtype);
3335: PetscStrcat(convname,"_C");
3336: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3337: if (conv) goto foundconv;
3338: }
3340: /* 2) See if a specialized converter is known to the desired matrix class. */
3341: MatCreate(((PetscObject)mat)->comm,&B);
3342: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N);
3343: MatSetType(B,newtype);
3344: for (i=0; i<3; i++) {
3345: PetscStrcpy(convname,"MatConvert_");
3346: PetscStrcat(convname,((PetscObject)mat)->type_name);
3347: PetscStrcat(convname,"_");
3348: PetscStrcat(convname,prefix[i]);
3349: PetscStrcat(convname,newtype);
3350: PetscStrcat(convname,"_C");
3351: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
3352: if (conv) {
3353: MatDestroy(B);
3354: goto foundconv;
3355: }
3356: }
3358: /* 3) See if a good general converter is registered for the desired class */
3359: conv = B->ops->convertfrom;
3360: MatDestroy(B);
3361: if (conv) goto foundconv;
3363: /* 4) See if a good general converter is known for the current matrix */
3364: if (mat->ops->convert) {
3365: conv = mat->ops->convert;
3366: }
3367: if (conv) goto foundconv;
3369: /* 5) Use a really basic converter. */
3370: conv = MatConvert_Basic;
3372: foundconv:
3373: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3374: (*conv)(mat,newtype,reuse,M);
3375: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3376: }
3377: PetscObjectStateIncrease((PetscObject)*M);
3378: return(0);
3379: }
3383: /*@C
3384: MatFactorGetSolverPackage - Returns name of the package providing the factorization routines
3386: Not Collective
3388: Input Parameter:
3389: . mat - the matrix, must be a factored matrix
3391: Output Parameter:
3392: . type - the string name of the package (do not free this string)
3394: Notes:
3395: In Fortran you pass in a empty string and the package name will be copied into it.
3396: (Make sure the string is long enough)
3398: Level: intermediate
3400: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable(), MatGetFactor()
3401: @*/
3402: PetscErrorCode MatFactorGetSolverPackage(Mat mat, const MatSolverPackage *type)
3403: {
3404: PetscErrorCode ierr;
3405: PetscErrorCode (*conv)(Mat,const MatSolverPackage*);
3410: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Only for factored matrix");
3411: PetscObjectQueryFunction((PetscObject)mat,"MatFactorGetSolverPackage_C",(void (**)(void))&conv);
3412: if (!conv) {
3413: *type = MAT_SOLVER_PETSC;
3414: } else {
3415: (*conv)(mat,type);
3416: }
3417: return(0);
3418: }
3422: /*@C
3423: MatGetFactor - Returns a matrix suitable to calls to MatXXFactorSymbolic()
3425: Collective on Mat
3427: Input Parameters:
3428: + mat - the matrix
3429: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3430: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3432: Output Parameters:
3433: . f - the factor matrix used with MatXXFactorSymbolic() calls
3435: Notes:
3436: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3437: such as pastix, superlu, mumps, spooles etc.
3439: PETSc must have been config/configure.py to use the external solver, using the option --download-package
3441: Level: intermediate
3443: .seealso: MatCopy(), MatDuplicate(), MatGetFactorAvailable()
3444: @*/
3445: PetscErrorCode MatGetFactor(Mat mat, const MatSolverPackage type,MatFactorType ftype,Mat *f)
3446: {
3447: PetscErrorCode ierr;
3448: char convname[256];
3449: PetscErrorCode (*conv)(Mat,MatFactorType,Mat*);
3455: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3456: MatPreallocated(mat);
3458: PetscStrcpy(convname,"MatGetFactor_");
3459: PetscStrcat(convname,((PetscObject)mat)->type_name);
3460: PetscStrcat(convname,"_");
3461: PetscStrcat(convname,type);
3462: PetscStrcat(convname,"_C");
3463: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3464: if (!conv) {
3465: PetscTruth flag;
3466: PetscStrcasecmp(MAT_SOLVER_PETSC,type,&flag);
3467: if (flag) {
3468: SETERRQ1(PETSC_ERR_SUP,"Matrix format %s does not have a built-in PETSc direct solver",((PetscObject)mat)->type_name);
3469: } else {
3470: 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);
3471: }
3472: }
3473: (*conv)(mat,ftype,f);
3474: return(0);
3475: }
3479: /*@C
3480: MatGetFactorAvailable - Returns a a flag if matrix supports particular package and factor type
3482: Collective on Mat
3484: Input Parameters:
3485: + mat - the matrix
3486: . type - name of solver type, for example, spooles, superlu, plapack, petsc (to use PETSc's default)
3487: - ftype - factor type, MAT_FACTOR_LU, MAT_FACTOR_CHOLESKY, MAT_FACTOR_ICC, MAT_FACTOR_ILU,
3489: Output Parameter:
3490: . flg - PETSC_TRUE if the factorization is available
3492: Notes:
3493: Some PETSc matrix formats have alternative solvers available that are contained in alternative packages
3494: such as pastix, superlu, mumps, spooles etc.
3496: PETSc must have been config/configure.py to use the external solver, using the option --download-package
3498: Level: intermediate
3500: .seealso: MatCopy(), MatDuplicate(), MatGetFactor()
3501: @*/
3502: PetscErrorCode MatGetFactorAvailable(Mat mat, const MatSolverPackage type,MatFactorType ftype,PetscTruth *flg)
3503: {
3504: PetscErrorCode ierr;
3505: char convname[256];
3506: PetscErrorCode (*conv)(Mat,MatFactorType,PetscTruth*);
3512: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3513: MatPreallocated(mat);
3515: PetscStrcpy(convname,"MatGetFactorAvailable_");
3516: PetscStrcat(convname,((PetscObject)mat)->type_name);
3517: PetscStrcat(convname,"_");
3518: PetscStrcat(convname,type);
3519: PetscStrcat(convname,"_C");
3520: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3521: if (!conv) {
3522: *flg = PETSC_FALSE;
3523: } else {
3524: (*conv)(mat,ftype,flg);
3525: }
3526: return(0);
3527: }
3532: /*@
3533: MatDuplicate - Duplicates a matrix including the non-zero structure.
3535: Collective on Mat
3537: Input Parameters:
3538: + mat - the matrix
3539: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
3540: values as well or not
3542: Output Parameter:
3543: . M - pointer to place new matrix
3545: Level: intermediate
3547: Concepts: matrices^duplicating
3549: .seealso: MatCopy(), MatConvert()
3550: @*/
3551: PetscErrorCode MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
3552: {
3554: Mat B;
3555: PetscInt i;
3561: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3562: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3563: MatPreallocated(mat);
3565: *M = 0;
3566: if (!mat->ops->duplicate) {
3567: SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
3568: }
3569: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
3570: (*mat->ops->duplicate)(mat,op,M);
3571: B = *M;
3572: if (mat->mapping) {
3573: MatSetLocalToGlobalMapping(B,mat->mapping);
3574: }
3575: if (mat->bmapping) {
3576: MatSetLocalToGlobalMappingBlock(B,mat->bmapping);
3577: }
3578: PetscMapCopy(((PetscObject)mat)->comm,mat->rmap,B->rmap);
3579: PetscMapCopy(((PetscObject)mat)->comm,mat->cmap,B->cmap);
3580:
3581: B->stencil.dim = mat->stencil.dim;
3582: B->stencil.noc = mat->stencil.noc;
3583: for (i=0; i<=mat->stencil.dim; i++) {
3584: B->stencil.dims[i] = mat->stencil.dims[i];
3585: B->stencil.starts[i] = mat->stencil.starts[i];
3586: }
3588: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
3589: PetscObjectStateIncrease((PetscObject)B);
3590: return(0);
3591: }
3595: /*@
3596: MatGetDiagonal - Gets the diagonal of a matrix.
3598: Collective on Mat and Vec
3600: Input Parameters:
3601: + mat - the matrix
3602: - v - the vector for storing the diagonal
3604: Output Parameter:
3605: . v - the diagonal of the matrix
3607: Level: intermediate
3609: Concepts: matrices^accessing diagonals
3611: .seealso: MatGetRow(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs()
3612: @*/
3613: PetscErrorCode MatGetDiagonal(Mat mat,Vec v)
3614: {
3621: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3622: if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3623: MatPreallocated(mat);
3625: (*mat->ops->getdiagonal)(mat,v);
3626: PetscObjectStateIncrease((PetscObject)v);
3627: return(0);
3628: }
3632: /*@
3633: MatGetRowMin - Gets the minimum value (of the real part) of each
3634: row of the matrix
3636: Collective on Mat and Vec
3638: Input Parameters:
3639: . mat - the matrix
3641: Output Parameter:
3642: + v - the vector for storing the maximums
3643: - idx - the indices of the column found for each row (optional)
3645: Level: intermediate
3647: Notes: The result of this call are the same as if one converted the matrix to dense format
3648: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
3650: This code is only implemented for a couple of matrix formats.
3652: Concepts: matrices^getting row maximums
3654: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(),
3655: MatGetRowMax()
3656: @*/
3657: PetscErrorCode MatGetRowMin(Mat mat,Vec v,PetscInt idx[])
3658: {
3665: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3666: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3667: MatPreallocated(mat);
3669: (*mat->ops->getrowmin)(mat,v,idx);
3670: PetscObjectStateIncrease((PetscObject)v);
3671: return(0);
3672: }
3676: /*@
3677: MatGetRowMinAbs - Gets the minimum value (in absolute value) of each
3678: row of the matrix
3680: Collective on Mat and Vec
3682: Input Parameters:
3683: . mat - the matrix
3685: Output Parameter:
3686: + v - the vector for storing the minimums
3687: - idx - the indices of the column found for each row (optional)
3689: Level: intermediate
3691: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
3692: row is 0 (the first column).
3694: This code is only implemented for a couple of matrix formats.
3696: Concepts: matrices^getting row maximums
3698: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMaxAbs(), MatGetRowMin()
3699: @*/
3700: PetscErrorCode MatGetRowMinAbs(Mat mat,Vec v,PetscInt idx[])
3701: {
3708: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3709: if (!mat->ops->getrowminabs) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3710: MatPreallocated(mat);
3711: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
3713: (*mat->ops->getrowminabs)(mat,v,idx);
3714: PetscObjectStateIncrease((PetscObject)v);
3715: return(0);
3716: }
3720: /*@
3721: MatGetRowMax - Gets the maximum value (of the real part) of each
3722: row of the matrix
3724: Collective on Mat and Vec
3726: Input Parameters:
3727: . mat - the matrix
3729: Output Parameter:
3730: + v - the vector for storing the maximums
3731: - idx - the indices of the column found for each row (optional)
3733: Level: intermediate
3735: Notes: The result of this call are the same as if one converted the matrix to dense format
3736: and found the minimum value in each row (i.e. the implicit zeros are counted as zeros).
3738: This code is only implemented for a couple of matrix formats.
3740: Concepts: matrices^getting row maximums
3742: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMaxAbs(), MatGetRowMin()
3743: @*/
3744: PetscErrorCode MatGetRowMax(Mat mat,Vec v,PetscInt idx[])
3745: {
3752: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3753: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3754: MatPreallocated(mat);
3756: (*mat->ops->getrowmax)(mat,v,idx);
3757: PetscObjectStateIncrease((PetscObject)v);
3758: return(0);
3759: }
3763: /*@
3764: MatGetRowMaxAbs - Gets the maximum value (in absolute value) of each
3765: row of the matrix
3767: Collective on Mat and Vec
3769: Input Parameters:
3770: . mat - the matrix
3772: Output Parameter:
3773: + v - the vector for storing the maximums
3774: - idx - the indices of the column found for each row (optional)
3776: Level: intermediate
3778: Notes: if a row is completely empty or has only 0.0 values then the idx[] value for that
3779: row is 0 (the first column).
3781: This code is only implemented for a couple of matrix formats.
3783: Concepts: matrices^getting row maximums
3785: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
3786: @*/
3787: PetscErrorCode MatGetRowMaxAbs(Mat mat,Vec v,PetscInt idx[])
3788: {
3795: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3796: if (!mat->ops->getrowmaxabs) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3797: MatPreallocated(mat);
3798: if (idx) {PetscMemzero(idx,mat->rmap->n*sizeof(PetscInt));}
3800: (*mat->ops->getrowmaxabs)(mat,v,idx);
3801: PetscObjectStateIncrease((PetscObject)v);
3802: return(0);
3803: }
3807: /*@
3808: MatGetRowSum - Gets the sum of each row of the matrix
3810: Collective on Mat and Vec
3812: Input Parameters:
3813: . mat - the matrix
3815: Output Parameter:
3816: . v - the vector for storing the maximums
3818: Level: intermediate
3820: Notes: This code is slow since it is not currently specialized for different formats
3822: Concepts: matrices^getting row sums
3824: .seealso: MatGetDiagonal(), MatGetSubMatrices(), MatGetSubmatrix(), MatGetRowMax(), MatGetRowMin()
3825: @*/
3826: PetscErrorCode MatGetRowSum(Mat mat, Vec v)
3827: {
3828: PetscInt start, end, row;
3829: PetscScalar *array;
3836: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3837: MatPreallocated(mat);
3838: MatGetOwnershipRange(mat, &start, &end);
3839: VecGetArray(v, &array);
3840: for(row = start; row < end; ++row) {
3841: PetscInt ncols, col;
3842: const PetscInt *cols;
3843: const PetscScalar *vals;
3845: array[row - start] = 0.0;
3846: MatGetRow(mat, row, &ncols, &cols, &vals);
3847: for(col = 0; col < ncols; col++) {
3848: array[row - start] += vals[col];
3849: }
3850: MatRestoreRow(mat, row, &ncols, &cols, &vals);
3851: }
3852: VecRestoreArray(v, &array);
3853: PetscObjectStateIncrease((PetscObject) v);
3854: return(0);
3855: }
3859: /*@
3860: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
3862: Collective on Mat
3864: Input Parameter:
3865: + mat - the matrix to transpose
3866: - reuse - store the transpose matrix in the provided B
3868: Output Parameters:
3869: . B - the transpose
3871: Notes:
3872: If you pass in &mat for B the transpose will be done in place
3874: Level: intermediate
3876: Concepts: matrices^transposing
3878: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose()
3879: @*/
3880: PetscErrorCode MatTranspose(Mat mat,MatReuse reuse,Mat *B)
3881: {
3887: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3888: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3889: if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
3890: MatPreallocated(mat);
3892: PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
3893: (*mat->ops->transpose)(mat,reuse,B);
3894: PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
3895: if (B) {PetscObjectStateIncrease((PetscObject)*B);}
3896: return(0);
3897: }
3901: /*@
3902: MatIsTranspose - Test whether a matrix is another one's transpose,
3903: or its own, in which case it tests symmetry.
3905: Collective on Mat
3907: Input Parameter:
3908: + A - the matrix to test
3909: - B - the matrix to test against, this can equal the first parameter
3911: Output Parameters:
3912: . flg - the result
3914: Notes:
3915: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
3916: has a running time of the order of the number of nonzeros; the parallel
3917: test involves parallel copies of the block-offdiagonal parts of the matrix.
3919: Level: intermediate
3921: Concepts: matrices^transposing, matrix^symmetry
3923: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
3924: @*/
3925: PetscErrorCode MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
3926: {
3927: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);
3933: PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
3934: PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
3935: if (f && g) {
3936: if (f==g) {
3937: (*f)(A,B,tol,flg);
3938: } else {
3939: SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
3940: }
3941: }
3942: return(0);
3943: }
3947: /*@
3948: MatIsHermitianTranspose - Test whether a matrix is another one's Hermitian transpose,
3950: Collective on Mat
3952: Input Parameter:
3953: + A - the matrix to test
3954: - B - the matrix to test against, this can equal the first parameter
3956: Output Parameters:
3957: . flg - the result
3959: Notes:
3960: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
3961: has a running time of the order of the number of nonzeros; the parallel
3962: test involves parallel copies of the block-offdiagonal parts of the matrix.
3964: Level: intermediate
3966: Concepts: matrices^transposing, matrix^symmetry
3968: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian(), MatIsTranspose()
3969: @*/
3970: PetscErrorCode MatIsHermitianTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
3971: {
3972: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);
3978: PetscObjectQueryFunction((PetscObject)A,"MatIsHermitianTranspose_C",(void (**)(void))&f);
3979: PetscObjectQueryFunction((PetscObject)B,"MatIsHermitianTranspose_C",(void (**)(void))&g);
3980: if (f && g) {
3981: if (f==g) {
3982: (*f)(A,B,tol,flg);
3983: } else {
3984: SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for Hermitian test");
3985: }
3986: }
3987: return(0);
3988: }
3992: /*@
3993: MatPermute - Creates a new matrix with rows and columns permuted from the
3994: original.
3996: Collective on Mat
3998: Input Parameters:
3999: + mat - the matrix to permute
4000: . row - row permutation, each processor supplies only the permutation for its rows
4001: - col - column permutation, each processor needs the entire column permutation, that is
4002: this is the same size as the total number of columns in the matrix
4004: Output Parameters:
4005: . B - the permuted matrix
4007: Level: advanced
4009: Concepts: matrices^permuting
4011: .seealso: MatGetOrdering()
4012: @*/
4013: PetscErrorCode MatPermute(Mat mat,IS row,IS col,Mat *B)
4014: {
4023: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4024: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4025: if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"MatPermute not available for Mat type %s",((PetscObject)mat)->type_name);
4026: MatPreallocated(mat);
4028: (*mat->ops->permute)(mat,row,col,B);
4029: PetscObjectStateIncrease((PetscObject)*B);
4030: return(0);
4031: }
4035: /*@
4036: MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the
4037: original and sparsified to the prescribed tolerance.
4039: Collective on Mat
4041: Input Parameters:
4042: + A - The matrix to permute
4043: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
4044: . frac - The half-bandwidth as a fraction of the total size, or 0.0
4045: . tol - The drop tolerance
4046: . rowp - The row permutation
4047: - colp - The column permutation
4049: Output Parameter:
4050: . B - The permuted, sparsified matrix
4052: Level: advanced
4054: Note:
4055: The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
4056: restrict the half-bandwidth of the resulting matrix to 5% of the
4057: total matrix size.
4059: .keywords: matrix, permute, sparsify
4061: .seealso: MatGetOrdering(), MatPermute()
4062: @*/
4063: PetscErrorCode MatPermuteSparsify(Mat A, PetscInt band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
4064: {
4065: IS irowp, icolp;
4066: const PetscInt *rows, *cols;
4067: PetscInt M, N, locRowStart, locRowEnd;
4068: PetscInt nz, newNz;
4069: const PetscInt *cwork;
4070: PetscInt *cnew;
4071: const PetscScalar *vwork;
4072: PetscScalar *vnew;
4073: PetscInt bw, issize;
4074: PetscInt row, locRow, newRow, col, newCol;
4075: PetscErrorCode ierr;
4082: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
4083: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
4084: if (!A->ops->permutesparsify) {
4085: MatGetSize(A, &M, &N);
4086: MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
4087: ISGetSize(rowp, &issize);
4088: if (issize != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for row permutation, should be %D", issize, M);
4089: ISGetSize(colp, &issize);
4090: if (issize != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for column permutation, should be %D", issize, N);
4091: ISInvertPermutation(rowp, 0, &irowp);
4092: ISGetIndices(irowp, &rows);
4093: ISInvertPermutation(colp, 0, &icolp);
4094: ISGetIndices(icolp, &cols);
4095: PetscMalloc(N * sizeof(PetscInt), &cnew);
4096: PetscMalloc(N * sizeof(PetscScalar), &vnew);
4098: /* Setup bandwidth to include */
4099: if (band == PETSC_DECIDE) {
4100: if (frac <= 0.0)
4101: bw = (PetscInt) (M * 0.05);
4102: else
4103: bw = (PetscInt) (M * frac);
4104: } else {
4105: if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
4106: bw = band;
4107: }
4109: /* Put values into new matrix */
4110: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
4111: for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
4112: MatGetRow(A, row, &nz, &cwork, &vwork);
4113: newRow = rows[locRow]+locRowStart;
4114: for(col = 0, newNz = 0; col < nz; col++) {
4115: newCol = cols[cwork[col]];
4116: if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
4117: cnew[newNz] = newCol;
4118: vnew[newNz] = vwork[col];
4119: newNz++;
4120: }
4121: }
4122: MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
4123: MatRestoreRow(A, row, &nz, &cwork, &vwork);
4124: }
4125: PetscFree(cnew);
4126: PetscFree(vnew);
4127: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
4128: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
4129: ISRestoreIndices(irowp, &rows);
4130: ISRestoreIndices(icolp, &cols);
4131: ISDestroy(irowp);
4132: ISDestroy(icolp);
4133: } else {
4134: (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
4135: }
4136: PetscObjectStateIncrease((PetscObject)*B);
4137: return(0);
4138: }
4142: /*@
4143: MatEqual - Compares two matrices.
4145: Collective on Mat
4147: Input Parameters:
4148: + A - the first matrix
4149: - B - the second matrix
4151: Output Parameter:
4152: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
4154: Level: intermediate
4156: Concepts: matrices^equality between
4157: @*/
4158: PetscErrorCode MatEqual(Mat A,Mat B,PetscTruth *flg)
4159: {
4169: MatPreallocated(B);
4170: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4171: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4172: 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);
4173: if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)A)->type_name);
4174: if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)B)->type_name);
4175: 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);
4176: MatPreallocated(A);
4178: (*A->ops->equal)(A,B,flg);
4179: return(0);
4180: }
4184: /*@
4185: MatDiagonalScale - Scales a matrix on the left and right by diagonal
4186: matrices that are stored as vectors. Either of the two scaling
4187: matrices can be PETSC_NULL.
4189: Collective on Mat
4191: Input Parameters:
4192: + mat - the matrix to be scaled
4193: . l - the left scaling vector (or PETSC_NULL)
4194: - r - the right scaling vector (or PETSC_NULL)
4196: Notes:
4197: MatDiagonalScale() computes A = LAR, where
4198: L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)
4200: Level: intermediate
4202: Concepts: matrices^diagonal scaling
4203: Concepts: diagonal scaling of matrices
4205: .seealso: MatScale()
4206: @*/
4207: PetscErrorCode MatDiagonalScale(Mat mat,Vec l,Vec r)
4208: {
4214: if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4217: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4218: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4219: MatPreallocated(mat);
4221: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4222: (*mat->ops->diagonalscale)(mat,l,r);
4223: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4224: PetscObjectStateIncrease((PetscObject)mat);
4225: return(0);
4226: }
4230: /*@
4231: MatScale - Scales all elements of a matrix by a given number.
4233: Collective on Mat
4235: Input Parameters:
4236: + mat - the matrix to be scaled
4237: - a - the scaling value
4239: Output Parameter:
4240: . mat - the scaled matrix
4242: Level: intermediate
4244: Concepts: matrices^scaling all entries
4246: .seealso: MatDiagonalScale()
4247: @*/
4248: PetscErrorCode MatScale(Mat mat,PetscScalar a)
4249: {
4255: if (a != 1.0 && !mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4256: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4257: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4258: MatPreallocated(mat);
4260: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4261: if (a != 1.0) {
4262: (*mat->ops->scale)(mat,a);
4263: PetscObjectStateIncrease((PetscObject)mat);
4264: }
4265: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4266: return(0);
4267: }
4271: /*@
4272: MatNorm - Calculates various norms of a matrix.
4274: Collective on Mat
4276: Input Parameters:
4277: + mat - the matrix
4278: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
4280: Output Parameters:
4281: . nrm - the resulting norm
4283: Level: intermediate
4285: Concepts: matrices^norm
4286: Concepts: norm^of matrix
4287: @*/
4288: PetscErrorCode MatNorm(Mat mat,NormType type,PetscReal *nrm)
4289: {
4297: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4298: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4299: if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4300: MatPreallocated(mat);
4302: (*mat->ops->norm)(mat,type,nrm);
4303: return(0);
4304: }
4306: /*
4307: This variable is used to prevent counting of MatAssemblyBegin() that
4308: are called from within a MatAssemblyEnd().
4309: */
4310: static PetscInt MatAssemblyEnd_InUse = 0;
4313: /*@
4314: MatAssemblyBegin - Begins assembling the matrix. This routine should
4315: be called after completing all calls to MatSetValues().
4317: Collective on Mat
4319: Input Parameters:
4320: + mat - the matrix
4321: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4322:
4323: Notes:
4324: MatSetValues() generally caches the values. The matrix is ready to
4325: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4326: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4327: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4328: using the matrix.
4330: Level: beginner
4332: Concepts: matrices^assembling
4334: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
4335: @*/
4336: PetscErrorCode MatAssemblyBegin(Mat mat,MatAssemblyType type)
4337: {
4343: MatPreallocated(mat);
4344: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
4345: if (mat->assembled) {
4346: mat->was_assembled = PETSC_TRUE;
4347: mat->assembled = PETSC_FALSE;
4348: }
4349: if (!MatAssemblyEnd_InUse) {
4350: PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
4351: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4352: PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
4353: } else {
4354: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
4355: }
4356: return(0);
4357: }
4361: /*@
4362: MatAssembled - Indicates if a matrix has been assembled and is ready for
4363: use; for example, in matrix-vector product.
4365: Collective on Mat
4367: Input Parameter:
4368: . mat - the matrix
4370: Output Parameter:
4371: . assembled - PETSC_TRUE or PETSC_FALSE
4373: Level: advanced
4375: Concepts: matrices^assembled?
4377: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
4378: @*/
4379: PetscErrorCode MatAssembled(Mat mat,PetscTruth *assembled)
4380: {
4385: *assembled = mat->assembled;
4386: return(0);
4387: }
4391: /*
4392: Processes command line options to determine if/how a matrix
4393: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
4394: */
4395: PetscErrorCode MatView_Private(Mat mat)
4396: {
4397: PetscErrorCode ierr;
4398: PetscTruth flg1,flg2,flg3,flg4,flg6,flg7,flg8;
4399: static PetscTruth incall = PETSC_FALSE;
4400: #if defined(PETSC_USE_SOCKET_VIEWER)
4401: PetscTruth flg5;
4402: #endif
4405: if (incall) return(0);
4406: incall = PETSC_TRUE;
4407: PetscOptionsBegin(((PetscObject)mat)->comm,((PetscObject)mat)->prefix,"Matrix Options","Mat");
4408: PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg1);
4409: PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg2);
4410: PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg3);
4411: PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg4);
4412: #if defined(PETSC_USE_SOCKET_VIEWER)
4413: PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg5);
4414: #endif
4415: PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg6);
4416: PetscOptionsName("-mat_view_draw","Draw the matrix nonzero structure","MatView",&flg7);
4417: PetscOptionsEnd();
4419: if (flg1) {
4420: PetscViewer viewer;
4422: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4423: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
4424: MatView(mat,viewer);
4425: PetscViewerPopFormat(viewer);
4426: }
4427: if (flg2) {
4428: PetscViewer viewer;
4430: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4431: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
4432: MatView(mat,viewer);
4433: PetscViewerPopFormat(viewer);
4434: }
4435: if (flg3) {
4436: PetscViewer viewer;
4438: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4439: MatView(mat,viewer);
4440: }
4441: if (flg4) {
4442: PetscViewer viewer;
4444: PetscViewerASCIIGetStdout(((PetscObject)mat)->comm,&viewer);
4445: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
4446: MatView(mat,viewer);
4447: PetscViewerPopFormat(viewer);
4448: }
4449: #if defined(PETSC_USE_SOCKET_VIEWER)
4450: if (flg5) {
4451: MatView(mat,PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4452: PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)mat)->comm));
4453: }
4454: #endif
4455: if (flg6) {
4456: MatView(mat,PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4457: PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)mat)->comm));
4458: }
4459: if (flg7) {
4460: PetscOptionsHasName(((PetscObject)mat)->prefix,"-mat_view_contour",&flg8);
4461: if (flg8) {
4462: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm),PETSC_VIEWER_DRAW_CONTOUR);
4463: }
4464: MatView(mat,PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4465: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4466: if (flg8) {
4467: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
4468: }
4469: }
4470: incall = PETSC_FALSE;
4471: return(0);
4472: }
4476: /*@
4477: MatAssemblyEnd - Completes assembling the matrix. This routine should
4478: be called after MatAssemblyBegin().
4480: Collective on Mat
4482: Input Parameters:
4483: + mat - the matrix
4484: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
4486: Options Database Keys:
4487: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
4488: . -mat_view_info_detailed - Prints more detailed info
4489: . -mat_view - Prints matrix in ASCII format
4490: . -mat_view_matlab - Prints matrix in Matlab format
4491: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
4492: . -display <name> - Sets display name (default is host)
4493: . -draw_pause <sec> - Sets number of seconds to pause after display
4494: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
4495: . -viewer_socket_machine <machine>
4496: . -viewer_socket_port <port>
4497: . -mat_view_binary - save matrix to file in binary format
4498: - -viewer_binary_filename <name>
4500: Notes:
4501: MatSetValues() generally caches the values. The matrix is ready to
4502: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
4503: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
4504: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
4505: using the matrix.
4507: Level: beginner
4509: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
4510: @*/
4511: PetscErrorCode MatAssemblyEnd(Mat mat,MatAssemblyType type)
4512: {
4513: PetscErrorCode ierr;
4514: static PetscInt inassm = 0;
4515: PetscTruth flg;
4521: inassm++;
4522: MatAssemblyEnd_InUse++;
4523: if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
4524: PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
4525: if (mat->ops->assemblyend) {
4526: (*mat->ops->assemblyend)(mat,type);
4527: }
4528: PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
4529: } else {
4530: if (mat->ops->assemblyend) {
4531: (*mat->ops->assemblyend)(mat,type);
4532: }
4533: }
4535: /* Flush assembly is not a true assembly */
4536: if (type != MAT_FLUSH_ASSEMBLY) {
4537: mat->assembled = PETSC_TRUE; mat->num_ass++;
4538: }
4539: mat->insertmode = NOT_SET_VALUES;
4540: MatAssemblyEnd_InUse--;
4541: PetscObjectStateIncrease((PetscObject)mat);
4542: if (!mat->symmetric_eternal) {
4543: mat->symmetric_set = PETSC_FALSE;
4544: mat->hermitian_set = PETSC_FALSE;
4545: mat->structurally_symmetric_set = PETSC_FALSE;
4546: }
4547: if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
4548: MatView_Private(mat);
4549: PetscOptionsHasName(((PetscObject)mat)->prefix,"-mat_is_symmetric",&flg);
4550: if (flg) {
4551: PetscReal tol = 0.0;
4552: PetscOptionsGetReal(((PetscObject)mat)->prefix,"-mat_is_symmetric",&tol,PETSC_NULL);
4553: MatIsSymmetric(mat,tol,&flg);
4554: if (flg) {
4555: PetscPrintf(((PetscObject)mat)->comm,"Matrix is symmetric (tolerance %G)\n",tol);
4556: } else {
4557: PetscPrintf(((PetscObject)mat)->comm,"Matrix is not symmetric (tolerance %G)\n",tol);
4558: }
4559: }
4560: }
4561: inassm--;
4562: return(0);
4563: }
4568: /*@
4569: MatCompress - Tries to store the matrix in as little space as
4570: possible. May fail if memory is already fully used, since it
4571: tries to allocate new space.
4573: Collective on Mat
4575: Input Parameters:
4576: . mat - the matrix
4578: Level: advanced
4580: @*/
4581: PetscErrorCode MatCompress(Mat mat)
4582: {
4588: MatPreallocated(mat);
4589: if (mat->ops->compress) {(*mat->ops->compress)(mat);}
4590: return(0);
4591: }
4595: /*@
4596: MatSetOption - Sets a parameter option for a matrix. Some options
4597: may be specific to certain storage formats. Some options
4598: determine how values will be inserted (or added). Sorted,
4599: row-oriented input will generally assemble the fastest. The default
4600: is row-oriented, nonsorted input.
4602: Collective on Mat
4604: Input Parameters:
4605: + mat - the matrix
4606: . option - the option, one of those listed below (and possibly others),
4607: - flg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE)
4609: Options Describing Matrix Structure:
4610: + MAT_SYMMETRIC - symmetric in terms of both structure and value
4611: . MAT_HERMITIAN - transpose is the complex conjugation
4612: . MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
4613: - MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag
4614: you set to be kept with all future use of the matrix
4615: including after MatAssemblyBegin/End() which could
4616: potentially change the symmetry structure, i.e. you
4617: KNOW the matrix will ALWAYS have the property you set.
4620: Options For Use with MatSetValues():
4621: Insert a logically dense subblock, which can be
4622: . MAT_ROW_ORIENTED - row-oriented (default)
4624: Note these options reflect the data you pass in with MatSetValues(); it has
4625: nothing to do with how the data is stored internally in the matrix
4626: data structure.
4628: When (re)assembling a matrix, we can restrict the input for
4629: efficiency/debugging purposes. These options include
4630: + MAT_NEW_NONZERO_LOCATIONS - additional insertions will be
4631: allowed if they generate a new nonzero
4632: . MAT_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
4633: . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
4634: . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
4635: - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
4637: Notes:
4638: Some options are relevant only for particular matrix types and
4639: are thus ignored by others. Other options are not supported by
4640: certain matrix types and will generate an error message if set.
4642: If using a Fortran 77 module to compute a matrix, one may need to
4643: use the column-oriented option (or convert to the row-oriented
4644: format).
4646: MAT_NEW_NONZERO_LOCATIONS set to PETSC_FALSE indicates that any add or insertion
4647: that would generate a new entry in the nonzero structure is instead
4648: ignored. Thus, if memory has not alredy been allocated for this particular
4649: data, then the insertion is ignored. For dense matrices, in which
4650: the entire array is allocated, no entries are ever ignored.
4651: Set after the first MatAssemblyEnd()
4653: MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
4654: that would generate a new entry in the nonzero structure instead produces
4655: an error. (Currently supported for AIJ and BAIJ formats only.)
4656: This is a useful flag when using SAME_NONZERO_PATTERN in calling
4657: KSPSetOperators() to ensure that the nonzero pattern truely does
4658: remain unchanged. Set after the first MatAssemblyEnd()
4660: MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
4661: that would generate a new entry that has not been preallocated will
4662: instead produce an error. (Currently supported for AIJ and BAIJ formats
4663: only.) This is a useful flag when debugging matrix memory preallocation.
4665: MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
4666: other processors should be dropped, rather than stashed.
4667: This is useful if you know that the "owning" processor is also
4668: always generating the correct matrix entries, so that PETSc need
4669: not transfer duplicate entries generated on another processor.
4670:
4671: MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
4672: searches during matrix assembly. When this flag is set, the hash table
4673: is created during the first Matrix Assembly. This hash table is
4674: used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
4675: to improve the searching of indices. MAT_NEW_NONZERO_LOCATIONS flag
4676: should be used with MAT_USE_HASH_TABLE flag. This option is currently
4677: supported by MATMPIBAIJ format only.
4679: MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
4680: are kept in the nonzero structure
4682: MAT_IGNORE_ZERO_ENTRIES - for AIJ/IS matrices this will stop zero values from creating
4683: a zero location in the matrix
4685: MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
4686: ROWBS matrix types
4688: Level: intermediate
4690: Concepts: matrices^setting options
4692: @*/
4693: PetscErrorCode MatSetOption(Mat mat,MatOption op,PetscTruth flg)
4694: {
4700: if (((int) op) < 0 || ((int) op) >= NUM_MAT_OPTIONS) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Options %d is out of range",(int)op);
4701: MatPreallocated(mat);
4702: switch (op) {
4703: case MAT_SYMMETRIC:
4704: mat->symmetric = flg;
4705: if (flg) mat->structurally_symmetric = PETSC_TRUE;
4706: mat->symmetric_set = PETSC_TRUE;
4707: mat->structurally_symmetric_set = flg;
4708: break;
4709: case MAT_HERMITIAN:
4710: mat->hermitian = flg;
4711: if (flg) mat->structurally_symmetric = PETSC_TRUE;
4712: mat->hermitian_set = PETSC_TRUE;
4713: mat->structurally_symmetric_set = flg;
4714: break;
4715: case MAT_STRUCTURALLY_SYMMETRIC:
4716: mat->structurally_symmetric = flg;
4717: mat->structurally_symmetric_set = PETSC_TRUE;
4718: break;
4719: case MAT_SYMMETRY_ETERNAL:
4720: mat->symmetric_eternal = flg;
4721: break;
4722: default:
4723: break;
4724: }
4725: if (mat->ops->setoption) {
4726: (*mat->ops->setoption)(mat,op,flg);
4727: }
4728: return(0);
4729: }
4733: /*@
4734: MatZeroEntries - Zeros all entries of a matrix. For sparse matrices
4735: this routine retains the old nonzero structure.
4737: Collective on Mat
4739: Input Parameters:
4740: . mat - the matrix
4742: Level: intermediate
4744: Concepts: matrices^zeroing
4746: .seealso: MatZeroRows()
4747: @*/
4748: PetscErrorCode MatZeroEntries(Mat mat)
4749: {
4755: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4756: if (mat->insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for matrices where you have set values but not yet assembled");
4757: if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4758: MatPreallocated(mat);
4760: PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
4761: (*mat->ops->zeroentries)(mat);
4762: PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
4763: PetscObjectStateIncrease((PetscObject)mat);
4764: return(0);
4765: }
4769: /*@C
4770: MatZeroRows - Zeros all entries (except possibly the main diagonal)
4771: of a set of rows of a matrix.
4773: Collective on Mat
4775: Input Parameters:
4776: + mat - the matrix
4777: . numRows - the number of rows to remove
4778: . rows - the global row indices
4779: - diag - value put in all diagonals of eliminated rows (0.0 will even eliminate diagonal entry)
4781: Notes:
4782: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4783: but does not release memory. For the dense and block diagonal
4784: formats this does not alter the nonzero structure.
4786: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4787: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4788: merely zeroed.
4790: The user can set a value in the diagonal entry (or for the AIJ and
4791: row formats can optionally remove the main diagonal entry from the
4792: nonzero structure as well, by passing 0.0 as the final argument).
4794: For the parallel case, all processes that share the matrix (i.e.,
4795: those in the communicator used for matrix creation) MUST call this
4796: routine, regardless of whether any rows being zeroed are owned by
4797: them.
4799: Each processor can indicate any rows in the entire matrix to be zeroed (i.e. each process does NOT have to
4800: list only rows local to itself).
4802: Level: intermediate
4804: Concepts: matrices^zeroing rows
4806: .seealso: MatZeroRowsIS(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4807: @*/
4808: PetscErrorCode MatZeroRows(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4809: {
4816: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4817: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4818: if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",((PetscObject)mat)->type_name);
4819: MatPreallocated(mat);
4821: (*mat->ops->zerorows)(mat,numRows,rows,diag);
4822: MatView_Private(mat);
4823: PetscObjectStateIncrease((PetscObject)mat);
4824: return(0);
4825: }
4829: /*@C
4830: MatZeroRowsIS - Zeros all entries (except possibly the main diagonal)
4831: of a set of rows of a matrix.
4833: Collective on Mat
4835: Input Parameters:
4836: + mat - the matrix
4837: . is - index set of rows to remove
4838: - diag - value put in all diagonals of eliminated rows
4840: Notes:
4841: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4842: but does not release memory. For the dense and block diagonal
4843: formats this does not alter the nonzero structure.
4845: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4846: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4847: merely zeroed.
4849: The user can set a value in the diagonal entry (or for the AIJ and
4850: row formats can optionally remove the main diagonal entry from the
4851: nonzero structure as well, by passing 0.0 as the final argument).
4853: For the parallel case, all processes that share the matrix (i.e.,
4854: those in the communicator used for matrix creation) MUST call this
4855: routine, regardless of whether any rows being zeroed are owned by
4856: them.
4858: Each processor should list the rows that IT wants zeroed
4860: Level: intermediate
4862: Concepts: matrices^zeroing rows
4864: .seealso: MatZeroRows(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4865: @*/
4866: PetscErrorCode MatZeroRowsIS(Mat mat,IS is,PetscScalar diag)
4867: {
4868: PetscInt numRows;
4869: const PetscInt *rows;
4876: ISGetLocalSize(is,&numRows);
4877: ISGetIndices(is,&rows);
4878: MatZeroRows(mat,numRows,rows,diag);
4879: ISRestoreIndices(is,&rows);
4880: return(0);
4881: }
4885: /*@C
4886: MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
4887: of a set of rows of a matrix; using local numbering of rows.
4889: Collective on Mat
4891: Input Parameters:
4892: + mat - the matrix
4893: . numRows - the number of rows to remove
4894: . rows - the global row indices
4895: - diag - value put in all diagonals of eliminated rows
4897: Notes:
4898: Before calling MatZeroRowsLocal(), the user must first set the
4899: local-to-global mapping by calling MatSetLocalToGlobalMapping().
4901: For the AIJ matrix formats this removes the old nonzero structure,
4902: but does not release memory. For the dense and block diagonal
4903: formats this does not alter the nonzero structure.
4905: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4906: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4907: merely zeroed.
4909: The user can set a value in the diagonal entry (or for the AIJ and
4910: row formats can optionally remove the main diagonal entry from the
4911: nonzero structure as well, by passing 0.0 as the final argument).
4913: Level: intermediate
4915: Concepts: matrices^zeroing
4917: .seealso: MatZeroRows(), MatZeroRowsLocalIS(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
4918: @*/
4919: PetscErrorCode MatZeroRowsLocal(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4920: {
4927: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4928: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4929: MatPreallocated(mat);
4931: if (mat->ops->zerorowslocal) {
4932: (*mat->ops->zerorowslocal)(mat,numRows,rows,diag);
4933: } else {
4934: IS is, newis;
4935: const PetscInt *newRows;
4937: if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
4938: ISCreateGeneral(PETSC_COMM_SELF,numRows,rows,&is);
4939: ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
4940: ISGetIndices(newis,&newRows);
4941: (*mat->ops->zerorows)(mat,numRows,newRows,diag);
4942: ISRestoreIndices(newis,&newRows);
4943: ISDestroy(newis);
4944: ISDestroy(is);
4945: }
4946: PetscObjectStateIncrease((PetscObject)mat);
4947: return(0);
4948: }
4952: /*@C
4953: MatZeroRowsLocalIS - Zeros all entries (except possibly the main diagonal)
4954: of a set of rows of a matrix; using local numbering of rows.
4956: Collective on Mat
4958: Input Parameters:
4959: + mat - the matrix
4960: . is - index set of rows to remove
4961: - diag - value put in all diagonals of eliminated rows
4963: Notes:
4964: Before calling MatZeroRowsLocalIS(), the user must first set the
4965: local-to-global mapping by calling MatSetLocalToGlobalMapping().
4967: For the AIJ matrix formats this removes the old nonzero structure,
4968: but does not release memory. For the dense and block diagonal
4969: formats this does not alter the nonzero structure.
4971: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS,PETSC_TRUE) the nonzero structure
4972: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4973: merely zeroed.
4975: The user can set a value in the diagonal entry (or for the AIJ and
4976: row formats can optionally remove the main diagonal entry from the
4977: nonzero structure as well, by passing 0.0 as the final argument).
4979: Level: intermediate
4981: Concepts: matrices^zeroing
4983: .seealso: MatZeroRows(), MatZeroRowsLocal(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
4984: @*/
4985: PetscErrorCode MatZeroRowsLocalIS(Mat mat,IS is,PetscScalar diag)
4986: {
4988: PetscInt numRows;
4989: const PetscInt *rows;
4995: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4996: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4997: MatPreallocated(mat);
4999: ISGetLocalSize(is,&numRows);
5000: ISGetIndices(is,&rows);
5001: MatZeroRowsLocal(mat,numRows,rows,diag);
5002: ISRestoreIndices(is,&rows);
5003: return(0);
5004: }
5008: /*@
5009: MatGetSize - Returns the numbers of rows and columns in a matrix.
5011: Not Collective
5013: Input Parameter:
5014: . mat - the matrix
5016: Output Parameters:
5017: + m - the number of global rows
5018: - n - the number of global columns
5020: Note: both output parameters can be PETSC_NULL on input.
5022: Level: beginner
5024: Concepts: matrices^size
5026: .seealso: MatGetLocalSize()
5027: @*/
5028: PetscErrorCode MatGetSize(Mat mat,PetscInt *m,PetscInt* n)
5029: {
5032: if (m) *m = mat->rmap->N;
5033: if (n) *n = mat->cmap->N;
5034: return(0);
5035: }
5039: /*@
5040: MatGetLocalSize - Returns the number of rows and columns in a matrix
5041: stored locally. This information may be implementation dependent, so
5042: use with care.
5044: Not Collective
5046: Input Parameters:
5047: . mat - the matrix
5049: Output Parameters:
5050: + m - the number of local rows
5051: - n - the number of local columns
5053: Note: both output parameters can be PETSC_NULL on input.
5055: Level: beginner
5057: Concepts: matrices^local size
5059: .seealso: MatGetSize()
5060: @*/
5061: PetscErrorCode MatGetLocalSize(Mat mat,PetscInt *m,PetscInt* n)
5062: {
5067: if (m) *m = mat->rmap->n;
5068: if (n) *n = mat->cmap->n;
5069: return(0);
5070: }
5074: /*@
5075: MatGetOwnershipRangeColumn - Returns the range of matrix columns owned by
5076: this processor.
5078: Not Collective, unless matrix has not been allocated, then collective on Mat
5080: Input Parameters:
5081: . mat - the matrix
5083: Output Parameters:
5084: + m - the global index of the first local column
5085: - n - one more than the global index of the last local column
5087: Notes: both output parameters can be PETSC_NULL on input.
5089: Level: developer
5091: Concepts: matrices^column ownership
5093: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), MatGetOwnershipRangesColumn()
5095: @*/
5096: PetscErrorCode MatGetOwnershipRangeColumn(Mat mat,PetscInt *m,PetscInt* n)
5097: {
5105: MatPreallocated(mat);
5106: if (m) *m = mat->cmap->rstart;
5107: if (n) *n = mat->cmap->rend;
5108: return(0);
5109: }
5113: /*@
5114: MatGetOwnershipRange - Returns the range of matrix rows owned by
5115: this processor, assuming that the matrix is laid out with the first
5116: n1 rows on the first processor, the next n2 rows on the second, etc.
5117: For certain parallel layouts this range may not be well defined.
5119: Not Collective, unless matrix has not been allocated, then collective on Mat
5121: Input Parameters:
5122: . mat - the matrix
5124: Output Parameters:
5125: + m - the global index of the first local row
5126: - n - one more than the global index of the last local row
5128: Note: both output parameters can be PETSC_NULL on input.
5130: Level: beginner
5132: Concepts: matrices^row ownership
5134: .seealso: MatGetOwnershipRanges(), MatGetOwnershipRangeColumn(), MatGetOwnershipRangesColumn()
5136: @*/
5137: PetscErrorCode MatGetOwnershipRange(Mat mat,PetscInt *m,PetscInt* n)
5138: {
5146: MatPreallocated(mat);
5147: if (m) *m = mat->rmap->rstart;
5148: if (n) *n = mat->rmap->rend;
5149: return(0);
5150: }
5154: /*@C
5155: MatGetOwnershipRanges - Returns the range of matrix rows owned by
5156: each process
5158: Not Collective, unless matrix has not been allocated, then collective on Mat
5160: Input Parameters:
5161: . mat - the matrix
5163: Output Parameters:
5164: . ranges - start of each processors portion plus one more then the total length at the end
5166: Level: beginner
5168: Concepts: matrices^row ownership
5170: .seealso: MatGetOwnershipRange(), MatGetOwnershipRangeColumn(), MatGetOwnershipRangesColumn()
5172: @*/
5173: PetscErrorCode MatGetOwnershipRanges(Mat mat,const PetscInt **ranges)
5174: {
5180: MatPreallocated(mat);
5181: PetscMapGetRanges(mat->rmap,ranges);
5182: return(0);
5183: }
5187: /*@C
5188: MatGetOwnershipRangesColumn - Returns the range of local columns for each process
5190: Not Collective, unless matrix has not been allocated, then collective on Mat
5192: Input Parameters:
5193: . mat - the matrix
5195: Output Parameters:
5196: . ranges - start of each processors portion plus one more then the total length at the end
5198: Level: beginner
5200: Concepts: matrices^column ownership
5202: .seealso: MatGetOwnershipRange(), MatGetOwnershipRangeColumn(), MatGetOwnershipRanges()
5204: @*/
5205: PetscErrorCode MatGetOwnershipRangesColumn(Mat mat,const PetscInt **ranges)
5206: {
5212: MatPreallocated(mat);
5213: PetscMapGetRanges(mat->cmap,ranges);
5214: return(0);
5215: }
5219: /*@
5220: MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
5221: Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
5222: to complete the factorization.
5224: Collective on Mat
5226: Input Parameters:
5227: + mat - the matrix
5228: . row - row permutation
5229: . column - column permutation
5230: - info - structure containing
5231: $ levels - number of levels of fill.
5232: $ expected fill - as ratio of original fill.
5233: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
5234: missing diagonal entries)
5236: Output Parameters:
5237: . fact - new matrix that has been symbolically factored
5239: Notes:
5240: See the users manual for additional information about
5241: choosing the fill factor for better efficiency.
5243: Most users should employ the simplified KSP interface for linear solvers
5244: instead of working directly with matrix algebra routines such as this.
5245: See, e.g., KSPCreate().
5247: Level: developer
5249: Concepts: matrices^symbolic LU factorization
5250: Concepts: matrices^factorization
5251: Concepts: LU^symbolic factorization
5253: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
5254: MatGetOrdering(), MatFactorInfo
5256: @*/
5257: PetscErrorCode MatILUFactorSymbolic(Mat fact,Mat mat,IS row,IS col,const MatFactorInfo *info)
5258: {
5268: if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %D",(PetscInt)info->levels);
5269: if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %G",info->fill);
5270: if (!(fact)->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",((PetscObject)mat)->type_name);
5271: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
5272: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
5273: MatPreallocated(mat);